EMPTYs with mixed dimension cause error

[paleolimbot]

The error message reported in paleolimbot/ggspatial#100 seems to be a result of mixed-dimension sfc objects, where calculating the Z range fails. I think this kind of thing can get generated from EMTPYs that weren't properly dimensioned:

sf::st_as_sfc(c("LINESTRING Z (1 2 3, 4 5 6)", "LINESTRING EMPTY"))
#> Error in CPL_get_z_range(obj, 1): z error - expecting three columns;

I can't seem to replicate the error reported in ggspatial, which occurs when running this:

sf::read_sf(system.file("longlake/LongLakeMarshWetlands.shp", package = "ggspatial"))
#> Simple feature collection with 18 features and 21 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XYZ
#> Bounding box:  xmin: 406919 ymin: 5081208 xmax: 414747.9 ymax: 5089171
#> z_range:       zmin: 4.5 zmax: 18.4
#> Projected CRS: NAD83 / UTM zone 20N
#> # A tibble: 18 × 22
#>    OBJECTID_1 OBJECTID FEAT_CODE HID   PROVKEY ZVALUE STARTDATE  ENDDATE   
#>         <dbl>    <int> <chr>     <chr> <chr>    <dbl> <date>     <date>    
#>  1     116918    55607 WASW40    <NA>  <NA>      9999 2009-10-01 NA        
#>  2     109511    48171 WASW40    <NA>  <NA>      9999 2009-10-01 NA        
#>  3     112880    51656 WASW40    <NA>  <NA>      9999 2009-10-01 NA        
#>  4     109515    48175 WASW40    <NA>  <NA>      9999 2009-10-01 NA        
#>  5     115493    55110 WASW40    <NA>  <NA>      9999 2009-10-01 NA        
#>  6     109511    48171 WASW40    <NA>  <NA>      9999 2009-10-01 NA        
#>  7          0        0 <NA>      <NA>  <NA>         0 1899-12-30 1899-12-30
#>  8     109484    48144 WASW40    <NA>  <NA>      9999 2009-10-01 NA        
#>  9     116913    55602 WASW40    <NA>  <NA>      9999 2009-10-01 NA        
#> 10     116828    55517 WASW40    <NA>  <NA>      9999 2009-10-01 NA        
#> 11     116788    55477 WASW40    <NA>  <NA>         9 2009-10-01 NA        
#> 12      13843    13673 WALK40    <NA>  <NA>         0 2009-10-09 NA        
#> 13     109481    48141 WASW40    <NA>  <NA>        11 2009-10-01 NA        
#> 14     116582    55271 WASW40    <NA>  <NA>         9 2009-10-01 NA        
#> 15     116914    55603 WASW40    <NA>  <NA>      9999 2009-10-01 NA        
#> 16     116904    55593 WASW40    <NA>  <NA>      9999 2009-10-01 NA        
#> 17     116913    55602 WASW40    <NA>  <NA>      9999 2009-10-01 NA        
#> 18          0        0 <NA>      <NA>  <NA>         0 1899-12-30 1899-12-30
#> # … with 14 more variables: PRODUCT <chr>, SCALE <chr>, COLLECTOR <chr>,
#> #   CAPTURE <chr>, PRODYEAR <chr>, PRODMONTH <chr>, X_Y_ACC <chr>, Z_ACC <chr>,
#> #   MINZ <dbl>, MAXZ <dbl>, POLY_CLASS <int>, SHAPE_LENG <dbl>,
#> #   SHAPE_LE_1 <dbl>, geometry <MULTIPOLYGON [m]>

The files have been the same for a long time and I can't spot anything funky with the shapefile (written using QGIS ages ago). Let me know if I can do any digging to help!

[edzer]

I can reproduce this; polygon 9 is "XY": when debugging in st_read:

Browse[3]> sapply(x$geometry, class)
     [,1]           [,2]           [,3]           [,4]           [,5]          
[1,] "XYZ"          "XYZ"          "XYZ"          "XYZ"          "XYZ"         
[2,] "MULTIPOLYGON" "MULTIPOLYGON" "MULTIPOLYGON" "MULTIPOLYGON" "MULTIPOLYGON"
[3,] "sfg"          "sfg"          "sfg"          "sfg"          "sfg"         
     [,6]           [,7]           [,8]           [,9]           [,10]         
[1,] "XYZ"          "XYZ"          "XYZ"          "XY"           "XYZ"         
[2,] "MULTIPOLYGON" "MULTIPOLYGON" "MULTIPOLYGON" "MULTIPOLYGON" "MULTIPOLYGON"
[3,] "sfg"          "sfg"          "sfg"          "sfg"          "sfg"         
     [,11]          [,12]          [,13]          [,14]          [,15]         
[1,] "XYZ"          "XYZ"          "XYZ"          "XYZ"          "XYZ"         
[2,] "MULTIPOLYGON" "MULTIPOLYGON" "MULTIPOLYGON" "MULTIPOLYGON" "MULTIPOLYGON"
[3,] "sfg"          "sfg"          "sfg"          "sfg"          "sfg"         
     [,16]          [,17]          [,18]         
[1,] "XYZ"          "XYZ"          "XYZ"         
[2,] "MULTIPOLYGON" "MULTIPOLYGON" "MULTIPOLYGON"
[3,] "sfg"          "sfg"          "sfg"         

funny enough, when looking at the output of

ogrinfo /usr/local/lib/R/site-library/ggspatial/longlake/LongLakeMarshWetlands.shp LongLakeMarshWetlands

and searching for feature 8 (counting starts at 0), it gives a MULTIPOLYGON Z with the proper (varying) Z coordinates ... I'm a bit at a loss here.

[rsbivand]

Please see my follow-up on paleolimbot/ggspatial#100 ; it looks as though the driver changed at some point.

[edzer]

Handling this more gracefully (and usefully!) on the sf side is an open issue, see #1592 and #1682, where I'm still hoping that @dcooley will come with a solution!

[rsbivand]

And paleolimbot/ggspatial#101 - converting to GPKG resolves problem.

[paleolimbot]

I'll leave the fix to this for those other two issues! Thanks!