sf：如何从 GEOMETRYCOLLECTION 回到 MULTIPOLYGON？答案

【问题标题】：sf: How to get back to MULTIPOLYGON from GEOMETRYCOLLECTION?sf：如何从 GEOMETRYCOLLECTION 回到 MULTIPOLYGON？
【发布时间】：2019-12-31 19:06:49
【问题描述】：

我有一个世界国家数据集，并希望将其拆分为本初子午线，并将数据重新集中在太平洋上。

我正在尝试使用简单功能 (sf) 来执行此操作，但遇到了我无法解决的对象类型问题。

为了拆分数据，我尝试了以下方法：


   st_wg84 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

   # world country layer
   sfpolys <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf") 
   %>% st_sfc(crs = st_wg84 )

   # shift central/prime meridian towards west 
   shift <- 152 

   # create "split line" to split worldmap (split at Prime Meridian)
   split.line <- st_linestring(
     x = cbind(matrix(shift-180, 181, 1), matrix(-90:90,181,1))
    ) %>% 
     st_sfc(crs=st_wg84)

   # split country polygons along prime meridian
   sfpolys.split <- lwgeom::st_split(sfpolys, split.line)

这会产生一个GEOMETRYCOLLECTION 对象，沿着所需的线分割，包含与传入MULTIPOLYGON 相同数量的特征。

接下来，我需要移动坐标以使地图重新居中，为此我必须将多边形坐标转换为数据框。

    countries <- data.table(map_data(as(sfpolys.split, "Spatial")))

    # Shift coordinates to fall correctly on shifted map
    countries$long.shift <- countries$long + shift
    countries$long.shift <- ifelse(countries$long.shift > 180, 
    countries$long.shift - 360, countries$long.shift)

    # plot shifted map
    ggplot() + 
      geom_polygon(data=countries, 
        aes(x=long.shift, y=lat, group=group), 
        colour="black", fill="gray80", size = 0.25) +
      coord_equal()

但是，此转换不适用于 GEOMETRYCOLLECTION，但适用于 MULTIPOLYGON。

所以为了回到MULTIPOLYGON，我首先尝试了以下操作：

sfpolys.split <- sfpolys.split %>% st_cast("MULTIPOLYGON")

但这会导致以下错误：“m[1, ] 中的错误：维数不正确”

然后我尝试了：

sfpolys.split <- sfpolys.split %>% st_collection_extract(type="POLYGON")

但这给出了一个POLYGON 对象，我不知道如何正确分组为MULTIPOLYGON。

有没有人知道进行这种拆分和转移的更好方法，或者从GEOMETRYCOLLECTION 到MULTIPOLYGON 的简单方法？

这是我想要的结果：

【问题讨论】：

sf 还有一个函数st_collection_extract() 可能会有所帮助。

标签： r polygon spatial sf

【解决方案1】：

GEOMETRYCOLLECTION 是一个几何体列表，因此我们可以提取单个几何体。

幸运的是，您的每个 GEOMETRYCOLLECTION 几何都是多边形，所以我们可以很好地将它们包装成 MULTIPOLYGONS

geoms <- lapply( sfpolys.split$geometry, `[` )
mp <- lapply( geoms, function(x) sf::st_multipolygon( x = x ) )

然后创建一个sfc

sfc_mp <- sf::st_sfc( mp )

并将其附加到您的对象上

sfpolys.split$mp <- sfc_mp
sfpolys.split <- sf::st_set_geometry( sfpolys.split, sfc_mp )

这是一个检查格陵兰岛已经分裂的阴谋。我在每个单独的多边形周围添加了一个白色边框

library(mapdeck)

sf_line <- sf::st_sf( geometry = split.line )

mapdeck() %>%
    add_path(
        data = sf_line
    ) %>%
    add_polygon(
        data = sfpolys.split
        , fill_colour = "name_pl"
        , stroke_colour = "#FFFFFF"
        , stroke_width = 50000
    )

您的绘图代码的其余部分不可重现，因此我将其留给您解决。

【讨论】：