3

In my python code I run the following command:

command = 'gdalwarp rgb.tif rgb_output_cut.tif -cutline extent.kml -dstnodata 0 -q'
os.system(command)

This cuts rgb.tif based on the extents specified in extent.kml and outputs a file called rgb_output_cut.tif. I would like to run this command without gdalwarp but instead with Python bindings. Is anyone familiar with this? Also, for the purposes of my task, having -dstnodata 0 is optional, it is just nice to have.

I found similar examples, however, they are just using gdalwarp for reprojection. For my case I'd like to cut out an extent.

Here, in my extent.kml file it is just a series of LinearRing Polygons, for example:

  <Placemark>
<Style><LineStyle><color>ff0000ff</color></LineStyle><PolyStyle><fill>0</fill></PolyStyle></Style>
  <Polygon><altitudeMode>clampToGround</altitudeMode><outerBoundaryIs><LinearRing><altitudeMode>clampToGround</altitudeMode><coordinates>-104.959386510378,39.7075208980859 -104.959394642483,39.7020277743643 -104.968738431214,39.7020465443243 -104.968779091739,39.7075083857725 -104.959386510378,39.7075208980859</coordinates></LinearRing></outerBoundaryIs></Polygon>

Not, I must use gdal version 1.10.1.

jlcv
  • 345
  • 5
  • 13
  • Extent, as in envelope? Consider GDAL_Translate -projwin Xmin Ymax Xmax Ymin (I think the parameters are in that order).. doesn't warp unless you tell it to adjust the pixel size so is much faster, you would just need a way to get the extent out of the KML. Instead of os.system have a look at subprocess.Popen, I find it works better. – Michael Stimson Nov 16 '17 at 04:53
  • My extent.kml file is just a series of LinearRing Polygons, I updated my post with an example. – jlcv Nov 16 '17 at 05:11
  • https://stackoverflow.com/questions/44975952/get-feature-extent-using-gdal-ogr use the KML or LIBKML driver, whichever works on your install. inDriver = ogr.GetDriverByName("KML") then the rest of it should work as normal.. that is if the layer from the kml doesn't support GetExtent http://gdal.org/python/osgeo.ogr.Layer-class.html#GetExtent – Michael Stimson Nov 16 '17 at 05:17

1 Answers1

9

Using the GDAL Python bindings (GDAL >= 2.1), it should be:

from osgeo import gdal

input_raster = "path/to/rgb.tif"
output_raster = "path/to/rgb_output_cut.tif"
input_kml = "path/to/extent.kml"

ds = gdal.Warp(output_raster,
              input_raster,
              format = 'GTiff',
              cutlineDSName = input_kml,
              cutlineLayer = 'extent',
              dstNodata = 0)
ds = None
Antonio Falciano
  • 14,333
  • 2
  • 36
  • 66