Exporting Multipolygon geojson with rgdal::writeOGR

On Thu, 22 Oct 2015, Andy Teucher wrote:

I?m finding that writeOGR isn?t exporting multipolygons properly using
the GeoJSON driver. I have a simple test case (borrowed from here:

http://gis.stackexchange.com/questions/137977/writeogr-alters-multipolygon-holes)
with a geojson string with one multipolygon containing two polygons. I
use readOGR to create a SpatialPolygonsDataFrame out of it, then write
it with writeOGR:

library(rgdal)

js <- '{
"type": "MultiPolygon",
"coordinates": [[[[102.0, 2.0], [103.0, 2.0], [103.0, 3.0], [102.0, 3.0],
[102.0, 2.0]]], [[[100.0, 0.0], [101.0, 0.0], [101.0, 1.0], [100.0, 1.0],
[100.0, 0.0]]]]
} '

spdf <- readOGR(js, layer='OGRGeoJSON', verbose=FALSE) # Create a
SpatialPoygonsDataFrame

temp <- tempfile()
writeOGR(spdf, dsn = temp, layer = "", driver="GeoJSON")
cat(readLines(temp))

# Output:
{ "type": "FeatureCollection", "crs": { "type": "name", "properties":
{ "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } }, "features": [ { "type":
"Feature", "id": 0, "properties": { "FID": 0 }, "geometry": { "type":
"MultiPolygon", "coordinates": [ [ [ [ 102.0, 2.0 ], [ 102.0, 3.0 ], [
103.0, 3.0 ], [ 103.0, 2.0 ], [ 102.0, 2.0 ] ], [ [ 100.0, 0.0 ], [
100.0, 1.0 ], [ 101.0, 1.0 ], [ 101.0, 0.0 ], [ 100.0, 0.0 ] ] ] ] } }
] }

I can replicate this with GDAL 2.0.1 and rgdal 1.0-7 - reading temp in gives
the orphaned hole. I'm surprised that passing the js object to readOGR()
doesn't fail, but that isn't the source of the problem. spdf appears to be
structured correctly. Please provide your GDAL and rgdal versions, and OS
details from sessionInfo().

We are completely relying on the GDAL drivers here - we can't cherry pick
for particular drivers (historically excepting ESRI Shapefile), so your
debugging will need to examine writeOGR(), the C code it calls, and the
interactions between the rgdal C code and called OGR functions. Please check
which other drivers are affected (I think KML is, is GML?). Can you place
debugging Rprintf(...); in the rgdal C code to see where this is coming
from?

I can start looking at this sometime next week, so please make a start
yourself; contributions from others are very welcome.

Roger

If you look closely at the output, you can see that the 'coordinates'
array now contains a single polygon array with two coordinate arrays:
the boundary, and a second one which is now treated as a hole of the
first (orphaned as it is outside the bounds of the polygon). The
original 'coordinates' array consists of two polygon arrays, each
consisting of a single coordinate array which defining a polygon (with
no holes), which is correct according the GeoJSON spec:
http://geojson.org/geojson-spec.html#polygon.

I'm always hesitant to call things a bug, but this doesn't appear to
happen using ogr2ogr on the command line:

writeOGR(spdf, ".", "test", driver = "ESRI Shapefile") # Write a
shapefile to convert using ogr2ogr
system("ogr2ogr -f GeoJSON test_from_shp.geojson test.shp")
cat(readLines("test_from_shp.geojson"))

# Output:
{ "type": "FeatureCollection",  "features": [ { "type": "Feature",
"properties": { "FID": 0 }, "geometry": { "type": "MultiPolygon",
"coordinates": [ [ [ [ 102.0, 2.0 ], [ 102.0, 3.0 ], [ 103.0, 3.0 ], [
103.0, 2.0 ], [ 102.0, 2.0 ] ] ], [ [ [ 100.0, 0.0 ], [ 100.0, 1.0 ],
[ 101.0, 1.0 ], [ 101.0, 0.0 ], [ 100.0, 0.0 ] ] ] ] } } ] }

The resulting output is correct.


Thanks in advance for any help on this.

Andy

--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
e-mail: Roger.Bivand at nhh.no