Exporting Multipolygon geojson with rgdal::writeOGR

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

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