8

EDIT: PARTIAL SOLUTION

Using gdal_translate in the command line seems to do the trick, even though the Python binding doesn't work.

This worked to crop a GeoTiff removing 300 pixel padding from top and left and keeping the next 2000x2000 pixels. gdal_translate -srcwin 300 300 2000 2000 input.tif output.tif

ORIGINAL QUESTION

I've spent an embarrassingly long time trying to figure this out.

I have a series of satellite images in GEOTiff format. Each image has a 300 pixel buffer where it overlaps with the image next to it.

Goal: I am trying to crop off the 300 pixel buffer off of each image and then use as a raster with GDAL.

Constraints: 1) I need to keep all of the metadata and coordinate system information associated with the files 2) I want to do this entirely in Python (no command line)

What I've tried unsuccessfully:

1) using gdal.translate's srcWin function:

raster_data = gdal.Open('image.tif')
x_size = raster_data.RasterXSize
y_size = raster_data.RasterYSize
raster_data_unpadded = gdal.Translate('temp.tif', raster_data,
                                       srcWin = [padding,padding, 
                                       x_size - padding*2, 
                                       y_size-padding*2])

The problem: This produces a black image with no data

2) Cropping image using PIL and then saving back as TIF

 from PIL import Image
 img = Image.open(image.tif)
 x_size, y_size = img.size
 box = (padding, padding, x_size-padding, y_size - padding)
 img_unpadded = img.crop(box)
 img_unpadded.save('unpadded_image.tif')

The problem: PIL fails to save the file. It just hangs. Trying to save as a ".tiff" file produces the error "encoder error -9"

3) Using Rasterio

with rasterio.open("image.tif") as src:

    out_image, out_transform = mask(src, geoms, crop=True)

out_meta = src.meta.copy()

The problem: Rasterio doesn't seem to accept masks in Pixels format (e.g. 300 pixels). Only takes geometries, like a polygon, in the coordinate system of the file. I don't know how to translate between Pixels and coordinates.

PAL
  • 81
  • 1
  • 1
  • 4
  • Can you elaborate on what you tried thusfar? If its a a few lines of code you should modify your question by editing this into it below bullet 3. – ZF007 Jan 30 '18 at 22:05
  • Any reason for the contempt towards command line? Contrary to the python gdal, it seems to be working – Val Jan 31 '18 at 09:22
  • Confirmed it works using gdal command line. `gdal_translate -srcwin 300 300 2000 2000 input.tif output.tif` [creates a cropped GeoTIFF 300 pixels cropped off top and left side and 2000 pixels high and long]. Weird the same thing doesn't work with the Python bindings.... – PAL Jan 31 '18 at 17:58

3 Answers3

5

The following should work:

import rasterio
from rasterio.windows import Window

with rasterio.open('input.tif') as src:
    window = Window(padding, padding, src.width - 2 * padding, src.height - 2 * padding)

    kwargs = src.meta.copy()
    kwargs.update({
        'height': window.height,
        'width': window.width,
        'transform': rasterio.windows.transform(window, src.transform)})

    with rasterio.open('cropped.tif', 'w', **kwargs) as dst:
        dst.write(src.read(window=window))
Sam
  • 171
  • 1
  • 4
3

This is how you use gdal_translate in python via gdal.Translate.

Best option: projwin

The easiest way is with the projwin flag, which takes 4 values:

window = (upper_left_x, upper_left_y, lower_right_x, lower_right_y)

These values are in map coordinates. The bounds of the input file can be obtained via gdalinfo input_raster.tif from the command line.

NOTE: for many coordinate systems, ymax is actually less than ymin, so it's important to use "upper_left" and "lower_right" to identify the coordinates instead of "max" and "min." Akin's answer didn't work for me because of this difference.

The complete solution, then is:

from osgeo import gdal

upper_left_x = 699934.584491
upper_left_y = 6169364.0093
lower_right_x = 700160.946739
lower_right_y = 6168703.00544
window = (upper_left_x,upper_left_y,lower_right_x,lower_right_y)

gdal.Translate('output_crop_raster.tif', 'input_raster.tif', projWin = window)

Additional option: srcwin

srcwin is another gdal_translate flag similar to projwin, but takes in the pixel and line window via an offset and size, instead of using the map coordinate bounds. You would use it like this, but OP seemed to have problems with this method, so better to just use projwin:

window = (offset_x, offset_y, size_x, size_y)
gdal.Translate('output_crop_raster.tif', 'input_raster.tif', srcWin = window)
Daniel C Jacobs
  • 691
  • 8
  • 18
0

Here is a script to remove the padding using Daniel's solution.

#!/usr/bin/env python3
'''
This script clips all files in a directory to remove the overlapping pixels.
It is based on https://stackoverflow.com/a/70777036/13001393.
'''

import os
import re
from osgeo import gdal

gdal.UseExceptions()
dirname = '/home/username/gisdata/ned/1m/'
filenames = os.listdir(dirname)
output = dirname + 'temp.tiff'

def flatten(l):
    return [item for sublist in l for item in sublist]

def GeoCoord(dataset, x, y):
    GT = dataset.GetGeoTransform()
    if GT:
        return [GT[0] + x*GT[1] + y*GT[2], GT[3] + x*GT[4] + y*GT[5]]

for file in filenames:
    match  = re.search(r'\.tif$', file)

    if match:
        path = dirname + file
        dataset = gdal.Open(path, gdal.GA_ReadOnly)
        size = [dataset.RasterXSize, dataset.RasterYSize]
        pad = [(size[0] % 100) / 2, (size[1] % 100) / 2]
        window = flatten([GeoCoord(dataset, pad[0], pad[1]),
                          GeoCoord(dataset, size[0] - pad[0], size[1] - pad[1])])
        gdal.Translate(output, path, projWin = window)
        os.remove(path)
        os.rename(output, path)
Cal Abel
  • 173
  • 1
  • 4