On Tue, 13 Aug 2019, Phil Haines wrote:
Dear list,
I have a single Polygons object containing multiple Polygon objects
that share a common border. When I output this using writeOGR() one of
the Polygon objects becomes a hole, as the following example shows.
Create a Polygons object containing two adjoining Polygon objects
library(rgdal)
r1 <- rbind(c(1,1),c(1,2),c(2,2),c(2,1),c(1,1))
r2 <- r1; r2[,1] <- r2[,1]+1
Ps1 = Polygons(list(Polygon(r1),Polygon(r2)),ID=1)
SPDF = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1)),
data.frame(Example=c("Minimal")))
Perform a write/readOGR() cycle
fn <- tempfile()
writeOGR(SPDF, fn, layer='test', driver='ESRI Shapefile')
SPDF2 <- readOGR(dsn=fn,layer='test')
Second Polygon object is now a hole
sapply(SPDF2@polygons[[1]]@Polygons,slot,"hole")
[1] FALSE TRUE
I see from the sp documentation that "Polygon objects belonging to a
Polygons object should either not overlap one-other, or should be
fully included" but I am not sure how this relates to bordering
Polygon objects. I would welcome any advice as to whether what I am
asking of writeOGR is reasonable?
The problem is with the 'ESRI Shapefile' representation and driver:
library(rgdal)
r1 <- rbind(c(1,1),c(1,2),c(2,2),c(2,1),c(1,1))
r2 <- r1; r2[,1] <- r2[,1]+1
Ps1 = Polygons(list(Polygon(r1),Polygon(r2)),ID=1)
SPDF = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1)),
data.frame(Example=c("Minimal")))
sapply(slot(slot(SPDF, "polygons")[[1]], "Polygons"), slot, "hole")
sapply(slot(slot(SPDF, "polygons")[[1]], "Polygons"), slot, "ringDir")
# which constructs a MULTIPOLYGON object:
rgeos::writeWKT(SPDF)
library(sf)
st_as_text(st_geometry(st_as_sf(SPDF)))
# The 'ESRI Shapefile' driver is not Simple-Feature compliant (it
# pre-dates it), so the failure occurs by seeing the second exterior ring
# as an interior ring
fn <- tempfile()
writeOGR(SPDF, fn, layer='test', driver='ESRI Shapefile')
SPDF2 <- readOGR(dsn=fn, layer='test')
sapply(slot(slot(SPDF2, "polygons")[[1]], "Polygons"), slot, "hole")
sapply(slot(slot(SPDF2, "polygons")[[1]], "Polygons"), slot, "ringDir")
rgeos::writeWKT(SPDF2)
st_as_text(st_geometry(st_as_sf(SPDF2)))
# This happens with sf too, using the same GDAL driver:
sf2 <- st_read(dsn=fn, layer='test')
st_geometry(sf2)
library(sf)
st_as_text(st_geometry(st_as_sf(sf2)))
rgeos::writeWKT(as(sf2, "Spatial"))
# Adding the comment fix doesn't help:
comment(slot(SPDF, "polygons")[[1]])
SPDF_c <- rgeos::createSPComment(SPDF)
comment(slot(SPDF_c, "polygons")[[1]])
writeOGR(SPDF_c, fn, layer='test_c', driver='ESRI Shapefile',
verbose=TRUE)
# reports
# Object initially classed as: wkbPolygon
# SFS comments in Polygons objects
# Object reclassed as: wkbMultiPolygon
SPDF2_c <- readOGR(dsn=fn, layer='test_c')
sapply(slot(slot(SPDF2_c, "polygons")[[1]], "Polygons"), slot, "hole")
sapply(slot(slot(SPDF2_c, "polygons")[[1]], "Polygons"), slot, "ringDir")
rgeos::writeWKT(SPDF2_c)
st_as_text(st_geometry(st_as_sf(SPDF2_c)))
# If the input object is written out with the GeoPackage driver:
fn1 <- tempfile(fileext=".gpkg")
writeOGR(SPDF, fn1, layer="test", driver='GPKG')
sf2a <- st_read(dsn=fn1,layer='test')
st_coordinates(st_geometry(sf2a))
SPDF2a <- readOGR(dsn=fn1)
sapply(slot(slot(SPDF2a, "polygons")[[1]], "Polygons"), slot, "hole")
sapply(slot(slot(SPDF2a, "polygons")[[1]], "Polygons"), slot, "ringDir")
rgeos::writeWKT(SPDF2a)
st_as_text(st_geometry(st_as_sf(SPDF2a)))
# the issue is resolved. If we separate the exterior rings:
r2a <- r1; r2a[,1] <- r2a[,1]+1.00001
Ps1a = Polygons(list(Polygon(r1),Polygon(r2a)),ID=1)
SPDFa = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1a)),
data.frame(Example=c("Minimal")))
fna <- tempfile()
writeOGR(SPDFa, fna, layer='test', driver='ESRI Shapefile')
SPDF2_a <- readOGR(dsn=fna, layer='test')
sapply(slot(slot(SPDF2_a, "polygons")[[1]], "Polygons"), slot, "hole")
sapply(slot(slot(SPDF2_a, "polygons")[[1]], "Polygons"), slot, "ringDir")
rgeos::writeWKT(SPDF2_a)
st_as_text(st_geometry(st_as_sf(SPDF2_a)))
# we are OK as the two exterior rings do not touch.
# does using sf make a difference?
fn_s <- tempfile(fileext=".shp")
st_write(st_as_sf(SPDF), dsn=fn_s)
sf_in <- st_read(fn_s)
st_as_text(st_geometry(st_as_sf(sf_in)))
# No
fn_s <- tempfile(fileext=".shp")
st_write(st_as_sf(SPDF_c), dsn=fn_s)
sf_in_c <- st_read(fn_s)
st_as_text(st_geometry(st_as_sf(sf_in_c)))
# nor with the pretend-SF-compliant comment set either.
So the weakness is in the "ESRI Shapefile" write driver, or possibly in
the OGRGeometryFactory::organizePolygons() function in GDAL used in
OGR_write() (a C++ function) called by writeOGR(). If sf::st_write() also
calls OGRGeometryFactory::organizePolygons(), we'd maybe consider that it
has a weakness for the "ESRI Shapefile" driver, but which does not affect
SF-compliant drivers.
