2

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

polygon on the raster

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.

enter image description here

  • Now, you can try to subset the raster to each individual feature of the shapefile and then access them separately. This code can be perfectly used in this case, but it is a little bit more complicate. If you have some issues ask for advice. – xunilk Sep 9 at 11:51
3

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):

enter image description here

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

enter image description here

Editing Note:

I added two more features to my original layer:

enter image description here

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

enter image description here

  • This way works, thanks. What if I have several features in my shape file? – Cz Gong Sep 5 at 7:04
  • It also works. Please, see my Editing Note. – xunilk Sep 5 at 9:15
  • Thanks, I got that. But I meant to subset the raster to each individual feature of the shape file, and then access them separately. – Cz Gong Sep 5 at 9:30
  • I know how to do that but your initial question was answered. You now need to loop all features for producing different clipped raster with above script, however, this is another question. – xunilk Sep 5 at 13:03

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.