Handling PostGIS's geojson output with a Z dimensions

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

Handling PostGIS's geojson output with a Z dimensions

Rutger
Dear all,

I have a Python script rasterizing polygons which are imported with OGR from GeoJSON. The GeoJSON is created by PostGIS and send over the network (direct connection to the db is not possible).

This has worked very well until it started to crash the Python interpreter, a hard crash with no traceback/error whatsoever. Which made it hard to debug, but after isolating the specific case i noticed it crashed on an empty polygon with a Z dimension. Rasterizing an empty polygon with GDAL is useless, but normally works fine.

If you convert an empty 2D MultiPolgon in PostGIS to GeoJSON it looks like:
'{"type": "MultiPolygon", "coordinates": []}'

The above is handled very well, ogr.CreateGeometryFromJson() converts this into a proper geom, and its IsEmpty() method returns True.

Now an empty 25D MultiPolygon from PostGIS looks like:
'{"type": "MultiPolygon", "coordinates": [[[]]]}'

Which ogr converts to a geom, but the IsEmpty() method retuns False now, which is why it slipped through and ended up in my gdal.Rasterize call. Exporting this geom to WKT shows 'MULTIPOLYGON EMPTY' as you would expect.

Passing this to gdal.Rasterize causes the crash, whereas it properly rasterizes the 2D empty polygon.

A few resulting questions;
Is the output from PostGIS valid GeoJSON? I briefly skimmed the geojson spec, but couldnt find something about nested empty arrays. If its valid, should gdal.Rasterize handle this without crashing? Or if I as an end user have to handle this, what is the appropriate way? Since IsEmpty() doesnt catch it, a workaround might be to use an additional Export/Import step running it through wkt for example:
geom = ogr.CreateGeometryFromJson('{"type": "MultiPolygon", "coordinates": [[]]}')
geom = ogr.CreateGeometryFromWkt(geom.ExportToWkt())

I don't have full control over the PostGIS queries on the back-end, so even though ST_FORCE2D() would also work, i would rather catch this on the front-end as well.

Here is a notebook replicating the case:
http://nbviewer.jupyter.org/gist/RutgerK/47412685882ed3f518c70fd09f57691c


Regards,
Rutger



Reply | Threaded
Open this post in threaded view
|

Re: Handling PostGIS's geojson output with a Z dimensions

Even Rouault-2

On mardi 9 mai 2017 01:46:00 CEST Rutger wrote:

> Dear all,

>

> I have a Python script rasterizing polygons which are imported with OGR from

> GeoJSON. The GeoJSON is created by PostGIS and send over the network

> (direct connection to the db is not possible).

>

> This has worked very well until it started to crash the Python interpreter,

> a hard crash with no traceback/error whatsoever. Which made it hard to

> debug, but after isolating the specific case i noticed it crashed on an

> empty polygon with a Z dimension. Rasterizing an empty polygon with GDAL is

> useless, but normally works fine.

 

This sounds pretty much like

https://trac.osgeo.org/gdal/ticket/6844

Fixed in 2.2.0 and in 2.1 branch (post 2.1.3)

 

Your workaround of exporting to WKT and re-importing from it is a good one.

 

Even

 

--

Spatialys - Geospatial professional services

http://www.spatialys.com


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Handling PostGIS's geojson output with a Z dimensions

Rutger
Even,

It does indeed, thanks for pointing that out. I guess i should have read the release notes before posting this. :)

I'll test it as soon as "conda-forge" has a GDAL build of 2.2 to confirm.

Congratulations on the release btw, lots of good features!


Regards,
Rutger