[gdal-dev] Python RasterizeLayer and gdal_rasterize

[Ivan Lucena]

Thanks Mike,

I found the problem.

When I call CreateLayer() I am passing the SRS but if I try to print it just before calling RasterizeLayer(), it is not there.

I tried to call FlushCache() and SyncToDisk() but the only thing that solve the problem is to close the dataset and open it again.

        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)​

        geom_ds.FlushCache() # That doesn't
        geom_ds.SyncToDisk()​ # help

        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()​
   ​
        print(geom_lyr.GetSpatialRef()) # None!!

        # That's the solution: close and open the dataset again
        geom_ds = None
        geom_ds = ogr.Open('E:\TMP\MEMORY_%s_GEOM.JSON' % name)​
        geom_lyr = geom_ds.GetLayerByName('MEMORY_LAYER')
        # todo: how can I do that in memory?

        print(geom_lyr.GetSpatialRef()) # Good!!

        err = gdal.RasterizeLayer(mem_ds, [1], geom_lyr, burn_values=[1])​
        if err != 0:​
            print(err)

Even, does it look like a bug to you? A GEOJSON driver bug maybe. I probably did something like that once in a raster driver, wait for the closing of the file to update the SRS.

Thanks,

Ivan



________________________________
From: Mike Taves <mwtoews at gmail.com>
Sent: Wednesday, July 24, 2019 4:40 PM
To: Ivan Lucena <ivan.lucena at outlook.com>
Cc: gdal-dev at lists.osgeo.org <gdal-dev at lists.osgeo.org>
Subject: Re: [gdal-dev] Python RasterizeLayer and gdal_rasterize

On Thu, 25 Jul 2019 at 01:19, Ivan Lucena <ivan.lucena at outlook.com> wrote:
>         mem_ds.SetGeoTransform((xmin, xResolut, 0, ymax, 0, yResolut,))

Check your value for yResolut, as it should be negative. In my
experience, a raster full of zeros is when the geometry is outside the
raster window, so be sure to manually check the raster and vector
geometries in (e.g.) QGIS. Lastly, to properly close the file, you
only need to dereference the object:

mem_ds = None
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20190724/9b7d6027/attachment-0001.html>