6

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

Improve this question
4
  • 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
    Commented Jul 26, 2016 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
    Commented Jul 26, 2016 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
    Commented Jul 27, 2016 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
    Commented Jul 27, 2016 at 5:19

2 Answers 2

Reset to default
17

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
Improve this answer
2
  • 1
    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
    Commented Jul 9, 2018 at 11:50
  • 4
    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
    Commented Jul 9, 2018 at 15:05
5

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.")
Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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