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=($)#($)(#$'])
_______________________________________________
gdal-dev mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
