sf:如何从 GEOMETRYCOLLECTION 返回 MULTIPOLYGON?
sf: How to get back to MULTIPOLYGON from GEOMETRYCOLLECTION?
我有一个世界国家数据集,我想在本初子午线上拆分它,然后将数据重新居中以关注太平洋。
我正在尝试使用简单功能 (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
的简单方法?
这是我想要的结果:
GEOMETRYCOLLECTION 是一个几何元素列表,因此我们可以提取各个几何元素。
幸运的是,您的每个 GEOMETRYCOLLECTION 几何图形都是 POLYGONS,因此我们可以很好地将它们包装成 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
)
您的其余绘图代码不可重现,因此我将其留给您整理。
我有一个世界国家数据集,我想在本初子午线上拆分它,然后将数据重新居中以关注太平洋。
我正在尝试使用简单功能 (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
的简单方法?
这是我想要的结果:
GEOMETRYCOLLECTION 是一个几何元素列表,因此我们可以提取各个几何元素。
幸运的是,您的每个 GEOMETRYCOLLECTION 几何图形都是 POLYGONS,因此我们可以很好地将它们包装成 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
)
您的其余绘图代码不可重现,因此我将其留给您整理。