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.
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
