Sf: Problem with MULTISURFACE sf object retrieved from geonode

[lbusett]

Hi,

I am currently trying to retrieve a POLYGON spatial object served by a geonode server, using something like:

vect <- sf::st_read(dsn = "http://XXXXX/geoserver/ows?service=wfs&version=2.0.0&request=GetCapabilities", layer = "geonode:bf_jolanda_appezzamenti_2018_1r")

(sorry, but I can not share the "real " address).

The problem is that the object is apparently retrieved correctly, but with a "MULTISURFACE" geometry type (in origin, before upload to geonode, it is a POLYGON):

Simple feature collection with 165 features and 8 fields
geometry type: MULTISURFACE
dimension: XY
bbox: xmin: 729592.6 ymin: 4968100 xmax: 737675.3 ymax: 4974834
epsg (SRID): 32632
proj4string: +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs
First 10 features:
gml_id id_geom azienda tenimento id_appez nome_appez anno raccolto the_geom
1 bf_jolanda_appezzamenti_2018_1r.1 0003 BONIFICHE FERRARESI JOLANDA 0003 NA 2018 primo MULTISURFACE (POLYGON ((731...
2 bf_jolanda_appezzamenti_2018_1r.2 0007 BONIFICHE FERRARESI JOLANDA 0007 NA 2018 primo MULTISURFACE (POLYGON ((731...

, and many sf functions do not appear to work properly on MULTISURFACE (e.g., plot, conversion to sp).

I tried to re-cast to POLYGON, but I got the following error:

```

sf::st_cast(vect, "MULTIPOLYGON")
Error in which_sfc_col(from_cls) : st_cast for MULTISURFACE not supported
```

, although the help of the function seems to suggest the existence of a method for that, and indeed it appears that the method exists here:

https://github.com/r-spatial/sf/blob/cd9acec231891c10ad44080081f89f38729067b7/R/cast_sfg.R#L202

Is there an alternative way to easily convert the MULTISURFACE object to a geometry type more easily managed by sf ?

[mdsumner]

Could you share at least a single sfg? I've never seen this type and don't know how to find examples. I've had good results with multipatch though, and guess this is similar but with triangulations rather than just flat patches?

[edzer]

See https://github.com/r-spatial/sf/commit/8162d66091b43b35dd05007a534816b93a0f2e62 ; example code:

library(sf)
demo(nc)
pol = st_geometry(nc)[1:2]
class(pol) = c("XY", "MULTISURFACE", "sfg")
pol
try(st_cast(pol, "MULTIPOLYGON")) # fails!
st_cast(st_sfc(pol), "MULTIPOLYGON")
st_cast(st_sf(a = 1, st_sfc(pol)), "MULTIPOLYGON")

library(sf)
# Linking to GEOS 3.6.2, GDAL 2.2.3, proj.4 4.9.3
demo(nc)


#   demo(nc)
#   ---- ~~

## this object was created as follows:
library(sf)

# nc = st_read(system.file("shapes/", package="maptools"), "sids")
# st_crs(nc) = 4267 # "+proj=longlat +ellps=clrk66" or "+proj=longlat +datum=NAD27"
# print(nc, n = 3)
# st_write(nc, "nc.gpkg", "nc.gpkg", driver = "GPKG")

# description of the dataset, see vignette in package spdep:
# https://cran.r-project.org/web/packages/spdep/vignettes/sids.pdf

datasource = { if ("GPKG" %in% st_drivers()$name)
    system.file("gpkg/nc.gpkg", package="sf")
else
    system.file("shape/nc.shp", package="sf")
}

agr = c(AREA = "aggregate", PERIMETER = "aggregate", CNTY_ = "identity",
        CNTY_ID = "identity", NAME = "identity", FIPS = "identity", FIPSNO = "identity",
        CRESS_ID = "identity", BIR74 = "aggregate", SID74 = "aggregate", NWBIR74 = "aggregate",
        BIR79 = "aggregate", SID79 = "aggregate", NWBIR79  = "aggregate")

nc = st_read(datasource, agr = agr)
# Reading layer `nc.gpkg' from data source `/home/edzer/R/x86_64-pc-linux-gnu-library/3.4/sf/gpkg/nc.gpkg' using driver `GPKG'
# Simple feature collection with 100 features and 14 fields
# Attribute-geometry relationship: 0 constant, 8 aggregate, 6 identity
# geometry type:  MULTIPOLYGON
# dimension:      XY
# bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
# epsg (SRID):    4267
# proj4string:    +proj=longlat +datum=NAD27 +no_defs
pol = st_geometry(nc)[1:2]
class(pol) = c("XY", "MULTISURFACE", "sfg")
pol
# MULTISURFACE (MULTIPOLYGON (((-81.47276 36.23436, -81.54084 36.27251, -81.56198 36.27359, -81.63306 36.34069, -81.74107 36.39178, -81.69828 36.47178, -81.7028 36.51934, -81.67 36.58965, -81.3453 36.57286, -81.34754 36.53791, -81.32478 36.51368, -81.31332 36.4807, -81.26624 36.43721, -81.26284 36.40504, -81.24069 36.37942, -81.23989 36.36536, -81.26424 36.35241, -81.32899 36.3635, -81.36137 36.35316, -81.36569 36.33905, -81.35413 36.29972, -81.36745 36.2787, -81.40639 36.28505, -81.41233 36.26729, -81.43104 36.26072, -81.45289 36.23959, -81.47276 36.23436))), MULTIPOLYGON (((-81.23989 36.36536, -81.24069 36.37942, -81.26284 36.40504, -81.26624 36.43721, -81.31332 36.4807, -81.32478 36.51368, -81.34754 36.53791, -81.3453 36.57286, -80.90344 36.56521, -80.93355 36.49831, -80.96577 36.46722, -80.94967 36.41473, -80.95639 36.4038, -80.97795 36.39138, -80.98284 36.37183, -81.00278 36.36668, -81.02464 36.37783, -81.0428 36.41034, -81.08425 36.42992, -81.09856 36.43115, -81.11331 36.42285, -81.12938 36.42633, -81.1384 36.41763, -81.15337 36.42474, -81.17667 36.41544, -81.23989 36.36536))))
try(st_cast(pol, "MULTIPOLYGON")) # fails!
# OGR: Corrupt data
# Error in CPL_multisurface_to_multipolygon(structure(list(x), crs = NA_crs_,  : 
#   OGR error
st_cast(st_sfc(pol), "MULTIPOLYGON")
# Geometry set for 2 features 
# geometry type:  MULTIPOLYGON
# dimension:      XY
# bbox:           xmin: -81.74107 ymin: 36.23436 xmax: -80.90344 ymax: 36.58965
# epsg (SRID):    NA
# proj4string:    NA
# MULTIPOLYGON (((-81.47276 36.23436, -81.54084 3...
# MULTIPOLYGON (((-81.23989 36.36536, -81.24069 3...
st_cast(st_sf(a = 1, st_sfc(pol)), "MULTIPOLYGON")
# Simple feature collection with 2 features and 1 field
# geometry type:  MULTIPOLYGON
# dimension:      XY
# bbox:           xmin: -81.74107 ymin: 36.23436 xmax: -80.90344 ymax: 36.58965
# epsg (SRID):    NA
# proj4string:    NA
#   a                       geometry
# 1 1 MULTIPOLYGON (((-81.47276 3...
# 2 1 MULTIPOLYGON (((-81.23989 3...
# Warning message:
# In st_cast.sf(st_sf(a = 1, st_sfc(pol)), "MULTIPOLYGON") :
#   repeating attributes for all sub-geometries for which they may not be constant

[lbusett]

I will have to ask if I can share it.... "privacy" reasons, sorry.

Do not know if this can help, but the WKT seems quite straightforward:

st_as_text(vect$the_geom[1]) [1] "MULTISURFACE (POLYGON ((731124.7 4974124, 730038.1 4974120, 730033.4 4974669, 730036.1 4974686, 730039.2 4974692, 730246.3 4974694, 730453.6 4974692, 730718.9 4974696, 730948.5 4974696, 731042.5 4974697, 731044.6 4974617, 731105.6 4974597, 731125.5 4974519, 731124.7 4974124)))"

seems like a "regular" POLYGON, just with an additional "MULTISURFACE" prefix...

[edzer]

Ah, this seems a bit more work: MULTISURFACE can contain both POLYGON or MULTIPOLYGON objects.

[lbusett]

@edzer

thanks for the reply. I will try that ASAP.

[mdsumner]

So, many of these "exotics" are just special case GCs, any plain explanation anywhere? I'll revisit the docs to see, this is good information. Was it apparent already to others?

[lbusett]

@edzer

unfortunately, the fix does not work, I figure due to the fact that my MULTISURFACE is a POLYGON.
However, changing this line :

https://github.com/r-spatial/sf/blob/8162d66091b43b35dd05007a534816b93a0f2e62/R/cast_sfc.R#L28

to:

MULTISURFACE = 3

and using:

st_cast(vect, "POLYGON")

seems to work, though I fear it may break if the MULTISURFACE is a MULTIPOLYGON .

BTW: here is a single geometry for testing:

    geom <- structure(list(structure(list(structure(list(structure(c(731124.744433883, 
730038.092612545, 730033.362977928, 730036.061967524, 730039.177946841, 
730246.308388667, 730453.557829005, 730718.85983322, 730948.518104378, 
731042.512398251, 731044.56433467, 731105.589862708, 731125.479666974, 
731124.744433883, 4974123.95500705, 4974119.67965301, 4974668.51650436, 
4974686.02837158, 4974691.94532496, 4974693.58219819, 4974691.91309645, 
4974695.65791511, 4974695.64078853, 4974696.54272954, 4974616.58933008, 
4974597.22944158, 4974518.74802878, 4974123.95500705), .Dim = c(14L, 
2L))), class = c("XY", "POLYGON", "sfg"))), class = c("XY", "MULTISURFACE", 
"sfg"))), class = c("sfc_MULTISURFACE", "sfc"), precision = 0, bbox = structure(c(730033.362977928, 
4974119.67965301, 731125.479666974, 4974696.54272954), .Names = c("xmin", 
"ymin", "xmax", "ymax"), class = "bbox"), crs = structure(list(
    epsg = 32632L, proj4string = "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs"), .Names = c("epsg", 
"proj4string"), class = "crs"), n_empty = 0L)

[edzer]

Using the fact that MULTISURFACE "is" a GEOMETRYCOLLECTION, this however seems to work:

st_cast(geom, "GEOMETRYCOLLECTION") %>% st_collection_extract("POLYGON")

[lbusett]

Indeed it works like a charm, both on sfc and sf objects !

Many thanks!

[edzer]

@mdsumner : other exotics CURVE, COMPOUNDCURVE, CIRCULARSTRING got gdal support recently (limited to converting to linear); this is now in sf::st_cast. These are not GCs.

[mdsumner]

Understood, I just didn't realize about the GC analogs. Thanks!

[rafaleo]

I've got something similar I guess:

library(rwfs)

dsn = "http://inspire.redcar-cleveland.gov.uk/geoserver/RCBC_INSPIRE_WFS/wfs?request=getfeature&version=2.0.0&typeName=RCBC_INSPIRE_WFS:RCBC-LANDSCAPE_CHARACTER_TRACT&outputformat=GML32"

fileName <- tempfile()
download.file(dsn, fileName)

request <- rwfs::GMLFile$new(fileName)
client <- rwfs::WFSCachingClient$new(request)

client$listLayers()
layer <- client$getLayer("RCBC-LANDSCAPE_CHARACTER_TRACT")

print(layer$the_geom)
plot(layer)
unlink(fileName)

st_cast(layer, "GEOMETRYCOLLECTION") %>% st_collection_extract("POLYGON")

This last line returns Error in which_sfc_col(from_cls) : st_cast for MULTISURFACE not supported. What I'm doing wrong? Should I change this line in sf/R/cast_sfc.R?

[rsbivand]

This is a closed issue. However, look at sf::gdal_utils() using the "vectortranslate" utility to access the freestanding ogr2ogr facilities, https://gdal.org/programs/ogr2ogr.html#ogr2ogr, including options for conversion from internally unsupported geometries. You'll need to put the downloaded features into a file I think (there was a case with CURVEPOLYGON objects but I can't find it now).

[rafaleo]

Thanks! You've lead me to the trails. I used:

> ogrinfo(fileName, so=TRUE)
[1] "Had to open data source read-only."                                               
[2] "INFO: Open of `C:\\Users\\cp\\AppData\\Local\\Temp\\RtmpojE2WU\\file15841aed7dab'"
[3] "      using driver `GML' successful."                                             
[4] "1: RCBC-LANDSCAPE_CHARACTER_TRACT (Multi Polygon)"                                

and
ogr2ogr(fileName, "sic.shp", "RCBC-LANDSCAPE_CHARACTER_TRACT", nlt = "POLYGON")
sic <- readOGR("sic.shp")
plot(sic)
```
and got a plot.