[R-sig-Geo] rgdal::writeOGR with driver='ESRI Shapefile' converts Polygon object into a hole

[Roger Bivand]

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

This is probably related to a similar but inverse problem with the 
SF-compliant GeoJSON driver in 2015:

https://stat.ethz.ch/pipermail/r-sig-geo/2015-October/023609.html

continued the next month in:

https://stat.ethz.ch/pipermail/r-sig-geo/2015-November/023656.html

The details are in this SVN diff

https://r-forge.r-project.org/scm/viewvc.php/pkg/src/OGR_write.cpp?root=rgdal&r1=555&r2=571

up to the end og the list thread, and this one from then until now:

https://r-forge.r-project.org/scm/viewvc.php/pkg/src/OGR_write.cpp?root=rgdal&r1=571&r2=733

Summary: could you change drivers, or is it really necessary to fix an EOL 
problem? What is your use case?

>
> Thanks in advance for your time,

Thanks for a complete example,

Roger

> Phil
>
>> sessionInfo()
> R version 3.5.1 (2018-07-02)
> Platform: x86_64-w64-mingw32/x64 (64-bit)
> Running under: Windows >= 8 x64 (build 9200)
>
> Matrix products: default
>
> locale:
> [1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United
> Kingdom.1252    LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C
>                         LC_TIME=English_United Kingdom.1252
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] rgdal_1.4-4 sp_1.3-1
>
> loaded via a namespace (and not attached):
> [1] compiler_3.5.1  tools_3.5.1     yaml_2.2.0      grid_3.5.1
> lattice_0.20-35
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using 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; e-mail: Roger.Bivand using nhh.no
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en