[GIS] How to clip raster inside of circle (python gdal)

clipgdalmaskingpolygonpython

I have a source raster and an ERSI shape file, want to crop the raster into the circle.

raster = gdal.Open("raster.tif", gdal.GA_ReadOnly) #read raster
projection=raster.GetProjectionRef()
VectorDriver = ogr.GetDriverByName('ESRI Shapefile') #intialize vector
VectorDataset = VectorDriver.Open("boxes.shp", 1) 
layer = VectorDataset.GetLayer()
feature=layer[0] #select the first polygon (the circle shown in image)
geom = feature.GetGeometryRef() 
minX, maxX, minY, maxY = geom.GetEnvelope()
OutTile = gdal.Warp("cut", raster, format='GTiff', outputBounds=[minX, minY, maxX, maxY],xRes=0.05,yRes=0.05,dstSRS=projection)

But when I plot the OutTile, it shows only the corners of the polygon. I also tried the "cutline" argument, but then it croped the whole raster.

Best Answer

If you create a clipping layer with only the circular feature to use, your code can be simplified in this way (I use my own paths):

from osgeo import gdal, ogr

OutTile = gdal.Warp("/home/zeito/pyqgis_data/cut.tif", 
                    "/home/zeito/pyqgis_data/utah_demUTM12.tif", 
                    cutlineDSName='/home/zeito/pyqgis_data/boxes.shp',
                    cropToCutline=True,
                    dstNodata = 0)

OutTile = None

On the other hand, you don't need to explicitly read raster or get drivers, datasets or ogr features. You can use creation parameters of resulting raster inside gdal.Warp function.

In my case, involved layers look in Map View of QGIS as follow (but it is not necessary to load them for running script):

After running above script, resulting layer loaded in QGIS looks as expected:

Editing Note:

I added two more features to my original layer:

and after running same script, resulting layer loaded in QGIS is also looking as expected:

Related Solutions

[GIS] How to create mask bands in GDAL and Python

for one particular band, here's how I create a nodata mask.

import gdal
import numpy as np

nodata = -9999#or set it based on the raster nodata value.
r = gdal.Open(rastername)
raster_arr = np.array(r.GetRasterBand(band).ReadAsArray())

nodatamask = raster_arr == nodata

#do your thing and end up with 
#a result_raster that you'd like to save

result_raster[nodatamask] = nodata
#then save.

Of course, if you're doing logarithms and you have negative values, you'd need to make a temporary nodata value.

nodatamask = raster_arr == nodata

raster_arr[nodatamask] = 100000000

#do your thing

result_raster[nodatamask] = nodata
#then save.

[GIS] How to modify raster pixel inside a polygon

Cutline does the exact opposite of what you are trying to do. If you use cutline with gdalwarp, it should crop your raster to the area INSIDE the polygon and set everything outside to no data (see below):

before after

I don't know why nothing happened in your case (perhaps you forgot to specify the cutline layer option -cl - who knows).

To do what you want there are several ways of going about it. One would be to create a polygon with the same extent as your raster and then calculate the difference of that polygon with your desired polygons. This will effectively give you a new polygon which is the area you want to keep. You can then use that as your cutline.

An alternative is to rasterize your polygons with the same resolution as your existing raster using gdal_rasterize, but set the -burn value to be the same as your NoData value (which you can also specify) and your -init value to be 1. You can then simply multiply your new raster with the existing one.

Best Answer

Related Solutions

[GIS] How to create mask bands in GDAL and Python

[GIS] How to modify raster pixel inside a polygon

Related Question