[R-sig-Geo] Exporting Multipolygon geojson with rgdal::writeOGR
Andy Teucher
andy.teucher at gmail.com
Tue Oct 27 00:10:25 CET 2015
Fantastic, thanks Roger. I will follow along in the r-forge repository.
Andy
On Mon, Oct 26, 2015 at 3:59 PM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> 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