Sf: Reading `gdb` file now fails after move from 0.4.3 to 0.5.2

Created on 17 Jul 2017 · 14Comments · Source: r-spatial/sf

The new CRAN release has broken some previously operational code with regards to reading a .gdb file (zipped gdb used in example attached).

packageVersion("sf")
#> [1] '0.4.3'
gdb <- "./test.gdb"
x <- sf::read_sf(gdb, layer = "example")
class(x)
#> [1] "sf"         "data.frame"
sf::st_bbox(x)
#>        xmin        ymin        xmax        ymax 
#> -179.999989   -3.730833  179.999989   31.797868

After updating to current CRAN version:

packageVersion("sf")
#> [1] '0.5.2'
gdb <- "./test.gdb"
x <- sf::read_sf(gdb, layer = "example")
#> Error in UseMethod("st_bbox"): no applicable method for 'st_bbox' applied to an object of class "c('XY', 'MULTISURFACE', 'sfg')"

test.gdb.zip

Source

adamdsmith

All 14 comments

I have no trouble reading this:

> x <- sf::st_read(gdb, layer = "example")
Reading layer `example' from data source `/tmp/test.gdb' using driver `OpenFileGDB'
Simple feature collection with 90 features and 30 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -180 ymin: -3.730833 xmax: 180 ymax: 31.79787
epsg (SRID):    4269
proj4string:    +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs

What do you get when you do

> sf::sf_extSoftVersion()
          GEOS           GDAL         proj.4         lwgeom 
       "3.5.1"        "2.1.3"        "4.9.2" "2.3.2 r15302"

? What is your sessionInfo() output? Maybe this is gdal 2.2 now supporting MULTISURFACE?

edzer on 17 Jul 2017

Could be, as I'm running gdal 2.2

sf::sf_extSoftVersion()
#>    GEOS    GDAL  proj.4  lwgeom 
#> "3.6.1" "2.2.0" "4.9.3"      NA

devtools::session_info()
#> Session info -------------------------------------------------------------
#>  setting  value                       
#>  version  R version 3.4.0 (2017-04-21)
#>  system   x86_64, mingw32             
#>  ui       RTerm                       
#>  language (EN)                        
#>  collate  English_United States.1252  
#>  tz       GMT                         
#>  date     2017-07-17
#> Packages -----------------------------------------------------------------
#>  package   * version date       source        
#>  backports   1.1.0   2017-05-22 CRAN (R 3.4.0)
#>  base      * 3.4.0   2017-04-21 local         
#>  compiler    3.4.0   2017-04-21 local         
#>  datasets  * 3.4.0   2017-04-21 local         
#>  DBI         0.7     2017-06-18 CRAN (R 3.4.0)
#>  devtools    1.13.2  2017-06-02 CRAN (R 3.4.0)
#>  digest      0.6.12  2017-01-27 CRAN (R 3.4.0)
#>  evaluate    0.10.1  2017-06-24 CRAN (R 3.4.0)
#>  graphics  * 3.4.0   2017-04-21 local         
#>  grDevices * 3.4.0   2017-04-21 local         
#>  grid        3.4.0   2017-04-21 local         
#>  htmltools   0.3.6   2017-04-28 CRAN (R 3.4.0)
#>  knitr       1.16    2017-05-18 CRAN (R 3.4.0)
#>  magrittr    1.5     2014-11-22 CRAN (R 3.4.0)
#>  memoise     1.1.0   2017-04-21 CRAN (R 3.4.0)
#>  methods   * 3.4.0   2017-04-21 local         
#>  Rcpp        0.12.11 2017-05-22 CRAN (R 3.4.0)
#>  rmarkdown   1.6     2017-06-15 CRAN (R 3.4.0)
#>  rprojroot   1.2     2017-01-16 CRAN (R 3.4.0)
#>  sf          0.5-2   2017-07-12 CRAN (R 3.4.1)
#>  stats     * 3.4.0   2017-04-21 local         
#>  stringi     1.1.5   2017-04-07 CRAN (R 3.4.0)
#>  stringr     1.2.0   2017-02-18 CRAN (R 3.4.0)
#>  tools       3.4.0   2017-04-21 local         
#>  udunits2    0.13    2016-11-17 CRAN (R 3.4.0)
#>  units       0.4-5   2017-06-15 CRAN (R 3.4.0)
#>  utils     * 3.4.0   2017-04-21 local         
#>  withr       1.0.2   2016-06-20 CRAN (R 3.4.0)
#>  yaml        2.1.14  2016-11-12 CRAN (R 3.4.0)

adamdsmith on 17 Jul 2017

Should work, now.

edzer on 19 Jul 2017

👍1

Thanks @edzer. Seems ready to close.

adamdsmith on 19 Jul 2017

Well, maybe not. Generating possible downstream complication with st_intersects. Working on a reprex.

adamdsmith on 19 Jul 2017

library(dplyr)
library(sf)
gdb <- "./test.gdb"
x <- sf::read_sf(gdb, layer = "example")

# Create some overlapping points
x_pts <- data.frame(id = 1:3,
                       lon = c(-169, -164, -160),
                       lat = c(16, 25, 0)) %>%
  st_as_sf(., coords = c("lon", "lat"), crs = st_crs(x))

st_intersects(x_pts, x)
#> although coordinates are longitude/latitude, it is assumed that they are planar
#> Error in CPL_geos_binop(st_geometry(x), st_geometry(y), op, par, pattern, : Evaluation error: ParseException: Unknown WKB type 12.

adamdsmith on 19 Jul 2017

That is a limitation of the underlying GEOS library, not of sf: GEOS supports types 1-7 AFAIK.

The question is what to do next? You can convert CIRCULARSTRING into LINESTRING with st_cast, but your data contains things like

> st_as_text(st_geometry(x)[47])
[1] "MULTISURFACE(CURVEPOLYGON(COMPOUNDCURVE(LINESTRING(-159.399779123 22.226016471, -159.399699153 22.226276431, -159.398736217 22.226118372, -159.398260872 22.226095318, -159.398140792 22.2260564590001, -159.398163058 22.2257268010001, -159.397882642 22.225394244, -159.397397157 22.225057335, -159.397318825 22.2251780230001, -159.396993115 22.225177984, -159.396748242 22.2248808800001, -159.396901679 22.224770398, -159.396876329 22.224673093, -159.399167008 22.224731392, -159.399502204 22.225551382), CIRCULARSTRING(-159.399502204 22.225551382, -159.399622762077 22.2257930044972, -159.399779123 22.226016471))))"

which is a bit more tricky.

edzer on 19 Jul 2017

Tricky indeed. Path of lease resistance seems to be to export the layer of interest from the gdb into shp and work on that...

Thanks again for your help.

x <- sf::read_sf("./example.shp")
x_pts <- data.frame(id = 1:3,
                    lon = c(-169, -164, -160),
                    lat = c(16, 25, 0)) %>%
  st_as_sf(., coords = c("lon", "lat"), crs = st_crs(x))

st_intersects(x_pts, x)
#> although coordinates are longitude/latitude, it is assumed that they are planar
#> [[1]]
#> [1] 79
#> 
#> [[2]]
#> [1] 78
#> 
#> [[3]]
#> [1] 90

adamdsmith on 19 Jul 2017

Interestingly, PostGIS can do this kind of thing:

postgis=# select * from st_intersects('POINT(-159.399699153 22.226276431)'::geometry, 'MULTISURFACE(CURVEPOLYGON(COMPOUNDCURVE(LINESTRING(-159.399779123 22.226016471, -159.399699153 22.226276431, -159.398736217 22.226118372, -159.398260872 22.226095318, -159.398140792 22.2260564590001, -159.398163058 22.2257268010001, -159.397882642 22.225394244, -159.397397157 22.225057335, -159.397318825 22.2251780230001, -159.396993115 22.225177984, -159.396748242 22.2248808800001, -159.396901679 22.224770398, -159.396876329 22.224673093, -159.399167008 22.224731392, -159.399502204 22.225551382), CIRCULARSTRING(-159.399502204 22.225551382, -159.399622762077 22.2257930044972, -159.399779123 22.226016471))))'::geometry) ;
 st_intersects 
---------------
 t
(1 row)

but I have no clue how it does it.

edzer on 19 Jul 2017

I think I found it: here and here - looks quite doable.

edzer on 19 Jul 2017

👍1

With your original data I now see

> xx = st_cast(x,  "MULTIPOLYGON")
> x_pts <- data.frame(id = 1:3,
+                     lon = c(-169, -164, -160),
+                     lat = c(16, 25, 0)) %>%
+   st_as_sf(., coords = c("lon", "lat"), crs = st_crs(x))
> 
> st_intersects(x_pts, xx)
[[1]]
[1] 79

[[2]]
[1] 78

[[3]]
[1] 90

edzer on 19 Jul 2017

👍1

Works perfectly on this end as well. Many thanks!

adamdsmith on 19 Jul 2017

Great, thanks!

edzer on 19 Jul 2017

I had a similar issue solved by casting to MULTIPOLYGON. You can check if your layer contains problematic MULTISURFACE geometries with: unique(st_geometry_type(st_geometry(foo)))