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

Roger Bivand Roger.Bivand at nhh.no
Mon Oct 26 23:59:06 CET 2015


On Mon, 26 Oct 2015, Andy Teucher wrote:

> Hi Roger,
>
> Sorry I haven't been much help here - it's not for lack of trying,
> just lack of skills. I've been poring through OGR_write.cpp and trying
> to figure out language and the flow of things... a few thoughts
> occurred to me (which are probably obvious to you):

Thanks, Andy.

We're making some progress but are not yet there. I'm uncertain about 
whether a layer can have Polygon and MultiPolygon features, but we've 
found an OGR function that aims to organise rings for which their 
exterior/interior status is not well known. This passes back the 
appropriate class, bur we then force it to MultiPolygon if any object is a 
MultiPolygon. When we have something to test, we'll get back to you.

Best wishes,

Roger

>
> 1) Is this bug because holes are standalone polygons in sp objects,
> while they are members of polygon objects in geojson, wkt, etc?
>
> 2) On a related note, I notice that a single polygon with a hole is
> identified as a multipolygon containing a single polygon with a hole,
> when using writeOGR and one of the affected drivers.  Practically I
> don't think this matters too much, but it's not quite right.
>
> 3) In lines 116-160 the polygon/multipolygon and line/multiline
> determination is made for the entire layer. Can/should this be made
> for each feature in the layer, as you could have a layer with a
> mixture of single polygons and multipolygons?
>
> 4) I feel like there ought to be one more layer of looping through the
> polyons. e.g.:
>
> loop over each feature;
>     determine whether each feature is a single polygon (length 1 or
> all rings but one are holes), or a multipolygon:
>        if a polygon loop over and add ring(s);
>        if a multipolygon loop over polygons:
>            for each polygon add ring(s) and somehow assign holes to
> parent polygons
>
> Again, apologies if I'm being totally obvious (or way off base).
> Thanks for all your work on these packages, and on this issue.
>
> Andy
>
>
> On Thu, Oct 22, 2015 at 11:07 PM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
>> On Fri, 23 Oct 2015, Andy Teucher wrote:
>>
>>> On my Mac, I am running rgdal 1.0-7 and GDAL 1.11.3
>>>
>>>
>>> sessionInfo()
>>>
>>
>> OK, thanks.
>>
>> The problem also affects the SQLite driver (the CRAN binaries do not have
>> this driver). Could someone please check the PostGIS driver - my guess is
>> that all OGC SFS compliant drivers may be affected.
>>
>> Could somebody also please check whether this is a recently introduced issue
>> or not? For example somebody still on rgdal < 1.0? src/OGR_write.cpp was
>> changed 2015-08-21, but the last previous change was 2015-06-11, then
>> 2015-05-31, 2014-08-17 ... from the ChangeLog visible on the package CRAN
>> page.
>>
>> Roger
>>
>>
>>>
>>>
>>>
>>>
>>> R version 3.2.2 (2015-08-14)
>>>
>>> Platform: x86_64-apple-darwin13.4.0 (64-bit)
>>>
>>> Running under: OS X 10.11 (El Capitan)
>>>
>>>
>>>
>>>
>>> locale:
>>>
>>> [1] en_CA.UTF-8/en_CA.UTF-8/en_CA.UTF-8/C/en_CA.UTF-8/en_CA.UTF-8
>>>
>>>
>>>
>>>
>>> attached base packages:
>>>
>>> [1] stats     graphics  grDevices utils     datasets  methods   base
>>>
>>>
>>>
>>>
>>> other attached packages:
>>>
>>> [1] rgdal_1.0-7 sp_1.2-1
>>>
>>>
>>>
>>>
>>> loaded via a namespace (and not attached):
>>>
>>> [1] tools_3.2.2     grid_3.2.2      lattice_0.20-33
>>>
>>>
>>>
>>>
>>> Andy
>>>
>>> On Thu, Oct 22, 2015 at 11:28 AM, Andy Teucher <andy.teucher at gmail.com>
>>> wrote:
>>>
>>>> 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
>>>
>>>
>>
>> --
>> 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
>

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


More information about the R-sig-Geo mailing list