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

[Andy Teucher]

Thanks very much Roger.

I've replicated this on my work machine (Windows) and on my Mac at home.

KML and GML are both affected. ESRI shapefile works fine and, though
I'm not very familiar with the format, 'MapInfo File' works too
(loading the file in qgis shows a multipart polygon with no geometry
errors).

sessionInfo() for my Windows machine (using rgdal 1.0-7 with included
gdal drivers):

R version 3.2.2 (2015-08-14)
Platform: i386-w64-mingw32/i386 (32-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

locale:
[1] LC_COLLATE=English_Canada.1252  LC_CTYPE=English_Canada.1252
LC_MONETARY=English_Canada.1252
[4] LC_NUMERIC=C                    LC_TIME=English_Canada.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] rgdal_1.0-7 sp_1.2-0

loaded via a namespace (and not attached):
[1] tools_3.2.2     grid_3.2.2      lattice_0.20-33


I'll send the details for my Mac environment this evening.

I'll dig into the source code as soon as I can - C is a completely
foreign language to me, so I can't promise anything fast!

Thanks again,
Andy

On Thu, Oct 22, 2015 at 12:27 AM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> 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
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
> --
> 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