3

I would like to know how can I clip a raster image using a bounding box in python.

So far I can do it with gdal_translate:

gdal_translate -projwin -75.3 5.5 -73.5 3.7 -of GTiff original.tif new.tif

I tried:

import gdal
from gdalconst import *
src_filename = 'ori.tif' 
dst_filename = 'new.tif'
src_ds = gdal.Open(src_filename, GA_ReadOnly)
format = "GTiff"
driver = gdal.GetDriverByName( format )
dst_ds = driver.CreateCopy( dst_filename, src_ds, 0, [ '-projwin=-75.3 5.5 -73.5 3.7'] )
dst_ds = None
src_ds = None

I know that the -projwin option should not be defined this way. I get the message:

Driver GTiff does not support -projwin creation option

  • 1
    Hi, if you can provide some of your code, that might be helpful for someone trying to improve on it. E.g. is the bounding box you mention in a shapefile, in WKT format, or based off something else? And maybe you can post a screenshot on what that shift looks like. – weiji14 Jul 26 '16 at 22:37
  • I think you will have to calculate a multiple of the cell size offset from the current extents so that it will not change the current pixels. – caiohamamura Jul 26 '16 at 22:52
  • Use -projwin ulx uly lrx lry or -srcwin xoff yoff xsize ysize in translate. Projwin units are goreferenced coordinates, srcwin units are rows and columns. Neither should resample the raster; internally -projwin extent is converted to -srcwin offset/count in the translate so either works exactly the same. I would use the python subprocess module to execute the command but could easily be done using QGIS GUI. That is assuming you're not wanting to do all the hard work yourself in python. – Michael Stimson Jul 27 '16 at 0:33
  • Note that the -projwin command only works properly in GDAL 1.11.4 or earlier. Not sure if its a bug or a deliberate change in 2.1.0. – obrl_soil Jul 27 '16 at 5:19
12

Since GDAL 2.1 (more info here) GDAL and OGR utilities can be used as library functions, so this task is incredibly simple now:

from osgeo import gdal

ds = gdal.Open('original.tif')
ds = gdal.Translate('new.tif', ds, projWin = [-75.3, 5.5, -73.5, 3.7])
ds = None
  • Very elegant solution. I wonder whether there is any way of storing the cropped raster in memory rather than copying then into new files. I would like to crop rasters on the fly to create numpy arrays stacking these cropped images for further analysis. Is this possible? Thanks! – Irene Jul 9 '18 at 11:50
  • 2
    I found the solution in the following Python code: svn.osgeo.org/gdal/trunk/autotest/utilities/…. Apparently, to load it in memory, "new.tif" would be an empty string, and then you can pass to the gdal.Translate() function a format="MEM" parameter as seen in the link pasted before. – Irene Jul 9 '18 at 15:05
4

You can use gdal.Warp or gdal.ReprojectImage to do this. Below, I have illustrated how gdal.Warp can be used to clip a raster to a bounding area obtained from a shapefile feature.

A list of Python command options for gdal.Warp are available here.

from osgeo import ogr, gdal

InputImage = 'XXX.tif'
Shapefile = 'XXX.shp'
RasterFormat = 'GTiff'
PixelRes = 0.5
VectorFormat = 'ESRI Shapefile'

# Open datasets
Raster = gdal.Open(InputImage, gdal.GA_ReadOnly)
Projection = Raster.GetProjectionRef()

VectorDriver = ogr.GetDriverByName(VectorFormat)
VectorDataset = VectorDriver.Open(Shapefile, 0) # 0=Read-only, 1=Read-Write
layer = VectorDataset.GetLayer()
FeatureCount = layer.GetFeatureCount()
print("Feature Count:",FeatureCount)

# Iterate through the shapefile features
Count = 0
for feature in layer:
    Count += 1
    print("Processing feature "+str(Count)+" of "+str(FeatureCount)+"...")

    geom = feature.GetGeometryRef() 
    minX, maxX, minY, maxY = geom.GetEnvelope() # Get bounding box of the shapefile feature

    # Create raster
    OutTileName = str(Count)+'.SomeTileName.tif'
    OutTile = gdal.Warp(OutTileName, Raster, format=RasterFormat, outputBounds=[minX, minY, maxX, maxY], xRes=PixelRes, yRes=PixelRes, dstSRS=Projection, resampleAlg=gdal.GRA_NearestNeighbour, options=['COMPRESS=DEFLATE'])
    OutTile = None # Close dataset

# Close datasets
Raster = None
VectorDataset.Destroy()
print("Done.")

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.