[R-sig-Geo] rgdal::writeOGR with driver='ESRI Shapefile' converts Polygon object into a hole
Roger Bivand
Roger@B|v@nd @end|ng |rom nhh@no
Sat Aug 17 11:06:02 CEST 2019
Please follow up both here and on:
https://github.com/r-spatial/sf/issues/1130
as the problem is also seen in the sf package using the same GDAL ESRI
Shapefile driver.
Roger
On Fri, 16 Aug 2019, Roger Bivand wrote:
> On Fri, 16 Aug 2019, Roger Bivand wrote:
>
>> 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.
>
> Without the comment set, OGRGeometryFactory::organizePolygons() is used; with
> it set, OGRGeometryFactory::organizePolygons() is not used, because the
> object is declared as two exterior rings. In both cases, we have the output
> object written out and read back in incorrectly with the ESRI shapefile
> driver, but SF-compliant drivers round-trip (in the test GeoJSON), correctly.
>
> It is likely that the changes made in 2015 to accommodate GeoJSON led to this
> possible regression for the ESRI Shapefile driver. I'm adding this geometry
> to tests/tripup.R and data files; without code changes the hole slot is wrong
> and the ring direction is changes to match, so the coordinates change too.
>
> Reading using the deprecated maptools::readShapeSpatial() also gets a hole in
> the second external ring. However, writing with deprecated
> maptools:: writeSpatialShape() yields a shapefile that when read with
> maptools:: readShapeSpatial() gives the correct two exterior ring, no hole
> object. When read with sf::st_read() and rgdal::readOGR(), the object is also
> correct. So the problem definitely lies in rgdal::writeOGR(), and
> sf::st_write() - roundtripping with sf::st_write() and sf::st_read() degrades
> from MULTIPOLYGON to POLYGON with the ESRI Shapefile driver.
>
> Roger
>
>>
>> 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
More information about the R-sig-Geo
mailing list