6

I am trying to get tiles of raster contained within a bounding polygon in a shapefile. My output raster should have values for all pixels within the polygon feature, and NA or 0 everywhere outside the bounding polygon. I am trying to use gdal.warp for this. This is the code which I tried, and while it works to trim the original raster to the extent of the polygon, I am not able to make the pixel values zero outside the polygon.

InputImage = 'inputs.tif'
Shapefile = 'outcrops.shp'
RasterFormat = 'GTiff'
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()

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

OutTileName = 'SomeTileName.tif'
OutTile = gdal.Warp(OutTileName, Raster, format=RasterFormat, outputBounds=[minX, minY, maxX, maxY], dstSRS=Projection, cutlineLayer=layer, cropToCutline=True)
OutTile = None # Close dataset

Raster = None
VectorDataset.Destroy()
print("Done.")

Can someone point me towards where I am going wrong?

2
  • With gdalwarp you would set nodata with -dstnodata or -dstalpha. Do something similar with your python.
    – user30184
    Oct 2, 2017 at 20:39
  • OutTile = gdal.Warp(OutTileName, Raster, format=RasterFormat, dstSRS=Projection, dstNodata=999, cutlineLayer=Shapefile, options=['COMPRESS=DEFLATE']) Tried this, but I am not the getting the clipped raster. Instead I am getting the original raster in the output file.
    – RRC
    Oct 2, 2017 at 21:06

2 Answers 2

6

The simplest solution I've been able to come up with makes use gdal.Warp in combination with OGR SQL. I've been tackling some days with this problem and just decided to give this a try, and it worked like a prayer:

import gdal

# Allow for verbose exception reporting
gdal.UseExceptions()

# Build the OGR SQL
sql = "SELECT * FROM your_shapefile"

# Clip the input Raster
result = gdal.Warp('path/to/clipped_raster',
                   '/path/to/raster',
                   cutlineDSName='/path/to/Vector',
                   cutlineSQL=sql,
                   cropToCutline=True)

# Initiate writing operations
result = None

As you can see, we don't read any datasets ourselves. We only provide the paths to gdal.Warp and the rest is taken care of. Keep in mind though, that I'm not aware of you shapefile dataset and the SQL-clause is only an example. There are some gotchas also, which are not directly related to the gdal.Warp, such as matching coordinate systems.

In case your shapefile would contain multiple instances, such as multiple crop fields, you could try implementing a WHERE-clause in to the SQL. Example of this would look like

SELECT * FROM your_shapefile WHERE a_field = 'a_value'

I truly hope this helps you!

5

This is easier to do with gdal_translate and gdal_rasterize (http://www.gdal.org/gdal_translate.html, http://www.gdal.org/gdal_rasterize.html)

#Read in datasets
Raster = gdal.Open('/path/to/raster')
VectorDriver = ogr.GetDriverByName(VectorFormat)
Vector = VectorDriver.Open('/path/to/.shp', 0)

#Get shapefile bounding box:
layer = Vector.GetLayer()
feature = layer.GetFeature(0)
geom = feature.GetGeometryRef()
minX, maxX, minY, maxY = geom.GetEnvelope()
bbox = (minX,maxY,maxX,minY) #Reorder bbox to use with gdal_translate

#Clip raster to shapefile bounding box
clip = gdal.Translate('path/to/output', Raster, projWin=bbox)

#Create new values
output = gdal.Rasterize(clip, [1], Vector, burn_values = [1], init_value = [0])

Translate is a better method of cutting rasters by a bounding box than warp because it's a little faster and uses slightly less memory to perform the operation.

The rasterize function is used to burn values into any raster data format according to a vector geometry. We initialize the raster with a value of 0 and assign a new value of 1 everywhere inside of the geometry.

I don't have data to test this code on, so I can't guarantee that it will work, but it should get you going on the right track.


It would also be relatively straightforward to throw the clipped image from gdal.Warp into numpy to add a value of 0 to all pixels outside of the geometry. Assuming we are starting with your output (OutTile):

band = OutTile.GetRasterBand(1) #Assuming we are working with one band output
array = band.ReadAsArray() #Numpy array
array[np.where(array!=1)] = 0 #Assuming all values inside polygon have value of 1
band.WriteArray(array)
band.FlushCache()
OutTile = None
1
  • Thank you for your response. I want the pixel values of the original raster to remain for pixels inside the polygon, and NoDATA outside the polygon. Your method makes it 1 inside, 0 outside. How do I keep the original raster values inside the polygon in the final output?
    – RRC
    Oct 3, 2017 at 15:41

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.