<!DOCTYPE html>
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
</head>
<body>
<p>William,<br>
</p>
<div class="moz-cite-prefix">Le 18/12/2024 à 13:05, Starms, William
(MU-Student) via gdal-dev a écrit :<br>
</div>
<blockquote type="cite"
cite="mid:FCB7608A-E1BC-4670-A070-37F08AEAD6D2@missouri.edu">
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<div dir="ltr">
<p class="MsoNormal"
style="font-size: 12pt; -webkit-text-size-adjust: auto; margin: 0in; font-family: Aptos, sans-serif;">
<span style="font-size: 11pt;">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 (<a
href="https://gdal.org/en/stable/api/python/utilities.html#osgeo.gdal.RasterizeLayer"
style="color: rgb(70, 120, 134);" moz-do-not-send="true"
class="moz-txt-link-freetext">https://gdal.org/en/stable/api/python/utilities.html#osgeo.gdal.RasterizeLayer</a>)
so I dug into the C++/CLI options (<a
href="https://gdal.org/en/latest/api/gdal_alg.html#_CPPv432GDALCreateGenImgProjTransformer212GDALDatasetH12GDALDatasetHPPc"
style="color: rgb(70, 120, 134);" moz-do-not-send="true"
class="moz-txt-link-freetext">https://gdal.org/en/latest/api/gdal_alg.html#_CPPv432GDALCreateGenImgProjTransformer212GDALDatasetH12GDALDatasetHPPc</a>)
but couldn’t get them working.</span></p>
</div>
</blockquote>
<p>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
<a class="moz-txt-link-freetext" href="https://github.com/OSGeo/gdal/blob/master/alg/gdalrasterize.cpp#L1659">https://github.com/OSGeo/gdal/blob/master/alg/gdalrasterize.cpp#L1659</a></p>
<p>There were 2 issues in your code:</p>
<p>- you didn't attach the SRS to the vector layer</p>
<p>- 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)</p>
<p>Fixed version below, with a modified geometry within the raster
extent and with an option to use a geotransform or GCPs:</p>
<p>"""<br>
</p>
<p>from osgeo import gdal, ogr, osr<br>
<br>
gdal.UseExceptions()<br>
<br>
def make_ogr_feature(spatial_ref, shape):<br>
rast_ogr_ds =
ogr.GetDriverByName('Memory').CreateDataSource('wrk')<br>
rast_mem_lyr = rast_ogr_ds.CreateLayer('poly',
srs=spatial_ref)<br>
feat = ogr.Feature(rast_mem_lyr.GetLayerDefn())<br>
feat.SetGeometryDirectly(ogr.Geometry(wkt=shape))<br>
rast_mem_lyr.CreateFeature(feat)<br>
return rast_mem_lyr, rast_ogr_ds<br>
<br>
driver = gdal.GetDriverByName('GTiff')<br>
ds = driver.Create('set.tif', xsize=100, ysize=100, bands=3,
eType=gdal.GDT_Byte)<br>
<br>
srs = osr.SpatialReference()<br>
srs.ImportFromEPSG(4326)<br>
srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) # ask
for longitude, latitude order<br>
<br>
layer, layer_ds = make_ogr_feature(srs, 'Polygon ((32 33,38 33,38
38,32 38,32 33))')<br>
<br>
use_geotransform = False<br>
if use_geotransform:<br>
ds.SetSpatialRef(srs)<br>
ds.SetGeoTransform([30,(40 - 30) / ds.RasterXSize,0,40,0, - (40
- 30) / ds.RasterYSize])<br>
else:<br>
ds.SetGCPs([gdal.GCP(30,40, 0, 0.5, 0.5, '', 'UpperLeft'),<br>
gdal.GCP(40,40, 0, 99.5, 0.5, '', 'UpperRight'),<br>
gdal.GCP(40,30,0,99.5,99.5,'','LowerRight'),<br>
gdal.GCP(30,30,0,0.5,99.5, '','LowerLeft')], srs)<br>
<br>
gdal.RasterizeLayer(ds, [1,2,3], layer, burn_values=[255,255,255])<br>
</p>
<p>"""<br>
</p>
<p>Even</p>
<span style="white-space: pre-wrap">
</span>
<pre class="moz-signature" cols="72">--
<a class="moz-txt-link-freetext" href="http://www.spatialys.com">http://www.spatialys.com</a>
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.</pre>
</body>
</html>