2

I am attempting to clip a raster using a polygon file (shapefile) and to save the result in memory. I succeed to do it using gdal.Warp and the option format=VRT using the following code. Please note that before using gdal.Warp, I use a function to fix conflict geometry of my shapefile.

 ## SAVING GDAL.WARP RESULT WITH VRT

from osgeo import ogr, gdal

# input data
raster = file.tif
raster = gdal.Open(raster)
shp = polygon.shp


# function to fix conflict geometries
def clean_shapefile(shp):
    shp = ogr.Open(shp, 1)
    layer = shp.GetLayer()

    for feature in layer:
        geom = feature.GetGeometryRef()
        if not geom.IsValid():
            feature.SetGeometry(geom.Buffer(0))
            layer.SetFeature(feature)
            assert feature.GetGeometryRef().IsValid()

    layer.ResetReading()
    assert all(feature.GetGeometryRef().IsValid() for feature in layer)


# function to clip raster by a polygon (shapefile)
def clip_raster(raster, shp):
    clean_shapefile(shp)  
    result = gdal.Warp("", raster, cutlineDSName=shp, cropToCutline=True, dstNodata=0, format="VRT")
        return result

# result
clipped_raster = clip_raster(raster, shp)

However, VRT is quiet slow and I want to save the clipped raster using the MEM driver of GDAL (apparently faster).

Problem: I get the following error when running gdal.Warp:

Warning 1: All options related to creation ignored in update mode.

Final result is 1 while I want a new gdal dataset. I believed it was due to the gdal.Open reading option (GA_ReadOnly = 0 VS GA_Update = 1), so I tried option 1 (GA_Update) but nothing changed. What is wrong with my code?

 ## SAVING GDAL.WARP RESULT WITH THE MEM DRIVER FROM GDAL

from osgeo import ogr, gdal

# input data
raster = file.tif
raster = gdal.Open(raster)
shp = polygon.shp


# function to fix conflict geometries
def clean_shapefile(shp):
    shp = ogr.Open(shp, 1)
    layer = shp.GetLayer()

    for feature in layer:
        geom = feature.GetGeometryRef()
        if not geom.IsValid():
            feature.SetGeometry(geom.Buffer(0))
            layer.SetFeature(feature)
            assert feature.GetGeometryRef().IsValid()

    layer.ResetReading()
    assert all(feature.GetGeometryRef().IsValid() for feature in layer)


# function to clip raster by a polygon (shapefile)
def clip_raster(raster, shp):
    clean_shapefile(shp)
    array = raster.GetRasterBand(1)
    array = array.ReadAsArray()
    x_pixels = array.shape[1]
    y_pixels = array.shape[0]
    driver = gdal.GetDriverByName("MEM")
    dataset = driver.Create("", x_pixels, y_pixels, 1)

    world_file = raster.GetGeoTransform()
    proj = raster.GetProjection()
    dataset.SetGeoTransform(world_file)
    dataset.SetProjection(proj)
    dataset.FlushCache()

    result = gdal.Warp(dataset, raster, cutlineDSName=shp, cropToCutline=True, dstNodata=0)
    return result

# result
clipped_raster = clip_raster(raster, shp)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.