[gdal-dev] Modify DST_SRS in Python RasterizeLayer

[Even Rouault]

William,

Le 18/12/2024 à 13:05, Starms, William (MU-Student) via gdal-dev a écrit :
>
> 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.
>
You can't directly pass transformer options that apply to 
GDALCreateGenImgProjTransformer() though GDALRasterizeLayer(). The 
former will use the later and compose transformer options with the logic 
at https://github.com/OSGeo/gdal/blob/master/alg/gdalrasterize.cpp#L1659

There were 2 issues in your code:

- you didn't attach the SRS to the vector layer

- the bounds of your raster minx,miny=(30,30)-maxx,maxy=(40,40) don't 
intersect the envelope of your geometry minx,miny=(3,5)-maxx,maxy=(24,29)

Fixed version below, with a modified geometry within the raster extent 
and with an option to use a geotransform or GCPs:

"""

from osgeo import gdal, ogr, osr

gdal.UseExceptions()

def make_ogr_feature(spatial_ref, shape):
     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

driver = gdal.GetDriverByName('GTiff')
ds = driver.Create('set.tif', xsize=100, ysize=100, bands=3, 
eType=gdal.GDT_Byte)

srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)  # ask for 
longitude, latitude order

layer, layer_ds = make_ogr_feature(srs, 'Polygon ((32 33,38 33,38 38,32 
38,32 33))')

use_geotransform = False
if use_geotransform:
   ds.SetSpatialRef(srs)
   ds.SetGeoTransform([30,(40 - 30) / ds.RasterXSize,0,40,0, - (40 - 30) 
/ ds.RasterYSize])
else:
   ds.SetGCPs([gdal.GCP(30,40, 0, 0.5, 0.5, '', 'UpperLeft'),
               gdal.GCP(40,40, 0, 99.5, 0.5, '', 'UpperRight'),
               gdal.GCP(40,30,0,99.5,99.5,'','LowerRight'),
               gdal.GCP(30,30,0,0.5,99.5, '','LowerLeft')], srs)

gdal.RasterizeLayer(ds, [1,2,3], layer, burn_values=[255,255,255])

"""

Even

-- 
http://www.spatialys.com
My software is free, but my time generally not.
Butcher of all kinds of standards, open or closed formats. At the end, this is just about bytes.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20241218/832dd916/attachment.htm>