[R-sig-Geo] Exporting Multipolygon geojson with rgdal::writeOGR

[Andy Teucher]

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 ] ] ] ] } }
] }

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