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"?

  • 1
    A polygon object as for example the Shapely.Polygon will not work, but an opened vector dataset (ds = ogr.Open(vectorfn)) should work. – Johan Dec 4 '19 at 15:36
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)
|improve this answer|||||

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.