R st_centroid GEOS error Desconocido WKB tipo 12

Question

Estoy tratando de calcular los centroides de los polígonos de las áreas de los códigos postales (archivo shape wfs), pero estoy recibiendo el siguiente mensaje de error:

Error in CPL_geos_op(“centroid”, x, numeric(0), integer(0), numeric(0), : Evaluation error: ParseException: Unknown WKB type 12.

¿Hay alguna forma de arreglar esto? No tengo ni idea de qué significa esto. Supongo que hay algo inusual con el formato de la columna de geometría (?).

> sessionInfo()
R version 4.0.4 (2021-02-15)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)

Matrix products: default

locale:
[1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United Kingdom.1252   
[3] LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C                           
[5] LC_TIME=English_United Kingdom.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ows4R_0.1-5     geometa_0.6-3   forcats_0.5.1   stringr_1.4.0   dplyr_1.0.5    
 [6] purrr_0.3.4     readr_1.4.0     tidyr_1.1.3     tibble_3.1.0    ggplot2_3.3.3  
[11] tidyverse_1.3.0 httr_1.4.2      sf_0.9-7        sp_1.4-5        RPostgres_1.3.1

Ejemplo reproducible (obtener datos del servidor WFS)

library(tidyverse)
library(sf)
library(httr)
library(ows4R)

# get postcode data ####

# following info from
# https://inbo.github.io/tutorials/tutorials/spatial_wfs_services/

# connect to Statistics Finland geoserver ####
wfs_regions <- "http://geo.stat.fi/geoserver/postialue/wfs"
regions_client <- WFSClient$new(wfs_regions, 
                                serviceVersion = "2.0.0")
regions_client$getFeatureTypes(pretty = TRUE)

url <- parse_url(wfs_regions)
url$query <- list(service = "wfs",
                  request = "GetFeature",
                  typename = "postialue:pno_2018",
                  srsName = "EPSG:3067"
                  )

request <- build_url(url)

# grab data ####
postal_codes <- read_sf(request)
# st_geometry(postal_codes)

ggplot(postal_codes) +
  geom_sf()

# subset to smaller area, e.g. #### 
pc.helsinki <- postal_codes[postal_codes$kunta == "091",]

ggplot(pc.helsinki) +
  geom_sf()

# centroids of postcode polygons ####
pc.helsinki$centroid <- st_centroid(pc.helsinki$geom)

Answer 1

El tipo WKB 12 significa MultiSurface https://en.wikipedia.org/wiki/Well-known_text_representation_of_geometry .

Y efectivamente ese servidor devuelve geometrías MultiSurface

http://geo.stat.fi/geoserver/postialue/wfs?service=WFS&version=2.0.0&request=GetFeature&typenames=postialue:pno_2018&srsName=EPSG:3067&count=2

Las geometrías no son realmente curvas porque dentro de MultiSurface sólo hay anillos lineales

<gml:MultiSurface srsName="http://www.opengis.net/gml/srs/epsg.xml#3067" srsDimension="2">
<gml:surfaceMember>
<gml:Polygon>
<gml:exterior>
<gml:LinearRing>
...

De todos modos, su software parece pensar que no sabe manejar las geometrías curvas y por lo tanto no tiene ni siquiera un intento con estos datos.

Sugiero guardar los datos en el disco como XML utilizando la URL anterior pero sin &count=2 límite. Entonces puede utilizar ogr2ogr https://gdal.org/programs/ogr2ogr.html y convertir los datos en algún formato que su software pueda leer con una opción -nlt CONVERT_TO_LINEAR . Creo que su software puede abrir ese conjunto de datos convertido.

Answer 2

Para manejar la conversión mencionada por @user30184 dentro de R, se puede utilizar la función ogr2ogr() proporcionada por mi gdalUtilities paquete.

Para demostrar su uso, he aquí una pequeña función de utilidad adaptada de una que he utilizado en el pasado cuando necesitaba convertir MULTISURFACE geometrías a MULTIPOLYGON geometrías.

library(gdalUtilities)

ensure_multipolygons <- function(X) {
    tmp1 <- tempfile(fileext = ".gpkg")
    tmp2 <- tempfile(fileext = ".gpkg")
    st_write(X, tmp1)
    ogr2ogr(tmp1, tmp2, f = "GPKG", nlt = "MULTIPOLYGON")
    Y <- st_read(tmp2)
    st_sf(st_drop_geometry(X), geom = st_geometry(Y))
}

## Try it on your data
pc.helsinki_2 <- ensure_multipolygons(pc.helsinki)
pc.helsinki_2_centroids <- st_centroid(pc.helsinki_2) 

Si desea asegurarse de que cada uno de los "centroides" se encuentre en algún lugar dentro del MULTIPOLYGON del que se deriva, puede utilizar en su lugar st_point_on_surface() así:

ggplot(pc.helsinki_2) +
    geom_sf() +
    geom_sf(data = st_point_on_surface(pc.helsinki_2))