[gdal-dev] Modify DST_SRS in Python RasterizeLayer

[Starms, William (MU-Student)]

Howdy! I’m trying to rasterize a shape that is in absolute pixel/line coordinates, but it gets burned wrong if the destination dataset has an srs due to coordinate conversion. There aren’t any docs for that specific function (https://gdal.org/en/stable/api/python/utilities.html#osgeo.gdal.RasterizeLayer) so I dug into the C++/CLI options (https://gdal.org/en/latest/api/gdal_alg.html#_CPPv432GDALCreateGenImgProjTransformer212GDALDatasetH12GDALDatasetHPPc) but couldn’t get them working.

If my layer’s srs is configured as None, it generates a warning that it’s assuming the destination srs. If I set it to an unconfigured srs (whatever that is), the warning goes away, but I assume it’s using an identity transform for the layer. I couldn’t find an srs I could apply to the layer that would disable coordinate conversion to the dst raster.

I know I could get the geotransform and convert it, if it exists, but my target images may not have a reference system and may not have had their gcps applied yet. I’m also hesitant to try and remove it and then reapply after burning because I’m unsure of what other effects it could have, such as with gcps.

Also, is there any chance there's a secret way to enable inverse mode on RasterizeLayer? Would be nice to know if there is, I don't have the experience with rasterize to generate an inverse and merge, I just end up inverting my input, but there's always a chance to get the raster bounds wrong that way.

Appreciate the help,

Will


Sample code snippets, not sure if this supports attachments:

from osgeo import gdal, ogr, osr
gdal.UseExceptions()
def make_ogr_feature(spatial_ref, shape):
    # Cobbled together from hours of stackoverflow and I don’t fully understand it, but it works
    rast_ogr_ds = ogr.GetDriverByName('Memory').CreateDataSource('wrk')
    rast_mem_lyr = rast_ogr_ds.CreateLayer('poly', srs=spatial_ref)
    feat = ogr.Feature(rast_mem_lyr.GetLayerDefn())
    feat.SetGeometryDirectly(ogr.Geometry(wkt=shape))
    rast_mem_lyr.CreateFeature(feat)
    return rast_mem_lyr, rast_ogr_ds

# create a blank image, no/blank/identity SRS
driver = gdal.GetDriverByName('GTiff')
ds = driver.Create('unset.tif', xsize=100, ysize=100, bands=3, eType=gdal.GDT_Byte)
# Create a layer with no srs
layer, layer_ds = make_ogr_feature(None, 'Polygon ((6 29, 24 22, 3 5, 6 29))')
# Burn to image, close
gdal.RasterizeLayer(ds, [1,2,3], layer, burn_values=[255,255,255])
# New image, set GCPs and SRS
ds = driver.Create('set.tif', xsize=100, ysize=100, bands=3, eType=gdal.GDT_Byte)
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
ds.SetGCPs([gdal.GCP(30, 30, 0, 0.5, 0.5, '', 'UpperLeft'),
            gdal.GCP(40, 30, 0, 99.5, 0.5, '', 'UpperRight'),
            gdal.GCP(40,40,0,99.5,99.5,'','LowerRight'),
            gdal.GCP(30,40,0,0.5,99.5, '','LowerLeft')],srs)
# try to rasterize, but nothing will render because it assumes dest srs
gdal.RasterizeLayer(ds, [1,2,3], layer, burn_values=[255,255,255])
# found some docs about transform options
gdal.RasterizeLayer(ds, [1,2,3], layer, burn_values=[255,255,255], options=['SRC_METHOD=NO_GEOTRANSFORM'])
gdal.RasterizeLayer(ds, [1,2,3], layer, burn_values=[255,255,255], options=['SRC_METHOD=NO_GEOTRANSFORM','DST_METHOD=NO_GEOTRANSFORM'])
# but it doesn't seem to do anything or be validated?
gdal.RasterizeLayer(ds, [1,2,3], layer, burn_values=[255,255,255], options=[SRC_METHOD=($)#($)(#$'])
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20241218/8df839f8/attachment-0001.htm>