Cropping a raster file using GDAL w/ Python
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.
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))
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)