1

I am new to ogr and vector data manipulation.

I am interested to cut an image (in this case a DEM file) using gdal.Warp() with the '-cutline' option to used the borders of a rectangular polygon to cut the underlying image. I saw good examples here

I wonder if it is possible to provide a "polygon object" instead of a "vector file"(.kml, .shp etc.) as the input to "cutlineDSName"?

user1771806
  • 61
  • 1
  • 1
  • 9

1 Answers1

1

I did that a few days ago. I am not an expert in Python so perhaps there are other ways to achieve that. Don't forget to import all the modules you need.

  • First, you can convert your DEM into an array of value 1 representing the extent of your DEM:
# input (dem)
dem = dem.tif

output (dem extent)


extent = extent.tif


convert the dem into an array of value 1


dem = gdal.Open(dem)
band = dem.GetRasterBand(1)
array = band.ReadAsArray()
array[:,:] = 1
x_pixels = array.shape[1]
y_pixels = array.shape[0]
driver = gdal.GetDriverByName("GTIFF")
dataset = driver.Create(extent, x_pixels, y_pixels, 1)
dataset.GetRasterBand(1).WriteArray(array)


georeferencing using the spatial reference of the input DEM


world_file = dem.GetGeoTransform()
proj = dem.GetProjection()
dataset.SetGeoTransform(world_file)
dataset.SetProjection(proj)
dataset.FlushCache()
  • Secondly, you can convert the resulted raster (extent.tif) into a polygon (.shp):
# input (dem extent)
extent = extent.tif

output (polygon in .shp)


shp = polygon.shp


generating the shapefile


extent = gdal.Open(extent)
band = extent.GetRasterBand(1)
driver = ogr.GetDriverByName("ESRI Shapefile")
dataset = driver.CreateDataSource(shp)


georeferencing and layer generation


srs = osr.SpatialReference()
srs.ImportFromWkt(extent.GetProjectionRef())
layer = dataset.CreateLayer("", srs)
field = ogr.FieldDefn("ID", ogr.OFTInteger)
layer.CreateField(field)
gdal.Polygonize(band, None, layer, 0, (), None)
dataset.Destroy()
  • Finaly, you can use the previous polygon to clip the DEM:
# input 1 (dem)
dem = dem.tif

input 2 (polygon in .shp)


shp = polygon.shp


output (clipped dem)


clipped_dem = clipped_dem.tif


clipping the DEM


with fiona.open(shp, "r") as shapefile:
    geometry = [feature["geometry"] for feature in shapefile]


with rasterio.open(dem) as src:
    out_image, out_transform = mask(src, geometry, crop=True)
    out_meta = src.meta.copy()
    out_meta.update({"driver": "GTiff",
                     "height": out_image.shape[1],
                     "width": out_image.shape[2],
                     "transform": out_transform})


with rasterio.open(clipped_dem, "w", **out_meta) as output:
    output.write(out_image)
rdmato33
  • 159
  • 7