如何从GEOMETRYCOLLECTION返回MULTIPOLYGON?

14

我有一个全球国家的数据集,希望能够在本初子午线上进行分割,并重新调整数据以关注太平洋。

我尝试使用简单要素(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")

但这会导致以下错误:“Error in m[1, ] : incorrect number of dimensions”

然后我尝试了:

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

但是这样只会得到一个 POLYGON 对象,我不知道如何正确地将其分组为一个 MULTIPOLYGON

有没有人知道更好的拆分和移动方法,或者从 GEOMETRYCOLLECTIONMULTIPOLYGON 的简单方法?

这是我想要的结果:

enter image description here


2
sf 还有一个函数 st_collection_extract() 可能会有帮助。 - RobertMyles
@RobertMyles st_collection_extract 只会提取 "POLYGON"、"POINT"、"LINESTRING" 中的一种。令人气愤的是,它不会提取任何其他类型的几何图形。 - Earlien
1个回答

7
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
    )

enter image description here

你剩下的绘图代码无法重现,所以我将让你自行解决。


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接