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

[Andy Teucher]

One more test: WKT has the same issue:

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)

writeOGR(spdf, dsn = "test.csv", layer = "test", driver="CSV",
layer_options = "GEOMETRY=AS_WKT")

cat(readLines("test.csv"), sep = "\n")

WKT,FID,
"MULTIPOLYGON (((102 2,102 3,103 3,103 2,102 2),(100 0,100 1,101 1,101
0,100 0)))",0


It should be:

WKT,FID,
"MULTIPOLYGON (((102 2,102 3,103 3,103 2,102 2)),((100 0,100 1,101
1,101 0,100 0)))",0

Andy

On Thu, Oct 22, 2015 at 10:36 AM, Andy Teucher <andy.teucher at gmail.com> wrote:
> 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