[gdal-dev] Python RasterizeLayer and gdal_rasterize

[Ivan Lucena]

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





-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20190724/583db5b2/attachment.html>