Hi there,
I am having a problem with gdal.RasterizeLayer() in Python. I look at the
documentation and code example in forums and tutorial but I can't find a
solution to my problem.
My code is supposed to loops thru a shape files, with countries for example,
and create a boolean image of each country but the images always come as all
black pixels. It should use GDAL and OGR memory drivers but in order to debug
the issue, I am saving one GEOJSON and one GTIFF file for each country.
What is most intriguing is that I can run gdal_rasterize with any of the
countries pairs of file and it works as expected. I open it on QGIS and it all
looks perfectly georeferenced too. The resolution is coming from another remote
sensing image that I am using as reference. And again, everything works fine
with the gdal_rasterize command, so it should be right.
Here is my code:
geom_ds =
ogr.GetDriverByName('GEOJSON').CreateDataSource('E:\TMP\MEMORY_%s.JSON' % name)
geom_lyr = geom_ds.CreateLayer('MEMORY_LAYER',
geom_type=ogr.wkbPolygon, srs=ref_srs)
geom_feat = ogr.Feature(geom_lyr.GetLayerDefn())
geom_feat.SetGeometry(geom)
geom_lyr.CreateFeature(geom_feat)
mem_ds = gdal.GetDriverByName('GTIFF').Create('E:\TMP\MEMORY_%s.TIF' %
name, xcount, ycount, 1, gdal.GDT_Byte)
mem_ds.SetGeoTransform((xmin, xResolut, 0, ymax, 0, yResolut,))
mem_ds.SetProjection(ref_ds.GetProjectionRef())
mem_ds.FlushCache()
err = gdal.RasterizeLayer(mem_ds, [1], geom_lyr, burn_values=[1])
if err != 0:
print(err)
Here is the command line:
(base) E:\TMP>gdal_rasterize MEMORY_Peru.JSON" "MEMORY_Peru.TIF" -burn 1
I am using GDAL 2.3.3 and Python 3.7.3 installed by Anaconda3 (64-bit). The
function gdal.UseExceptions() was called but there is nothing been reported.
I might be asking too much, but if someone had the same experience and have
found a solution and want to share, that would be awesome. Or maybe there is an
obvious error on my code that you can point out. I tried several things like
flushing and closing and re-opening the files, to make sure that the data is
ready for the functions, but nothing really helped.
One other thing is, once I got it to work, I am thinking that I don't need to
prepare the little geometry file if I could call gdal.rasterizelayer() the
original shape with a *query* for each country at the time. That should be more
efficient, I guess.
Thank you for your attention and have a nice day everybody,
Ivan
_______________________________________________
gdal-dev mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
