从状态组创建区域经纬度数据框答案

【问题标题】：Create regional lat, long dataframe from group of states从状态组创建区域经纬度数据框
【发布时间】：2020-12-09 20:17:23
【问题描述】：

我一直在使用urbnmapr 进行所有地图制作，它包含用于绘制美国各州的漂亮边界。我现在需要几个州的外边界来制作区域地图，并且正在寻找一种方法来过滤数据框以仅包含外边界。

我将使用这个外边框来遮盖光栅图像。我正在合作的地区是

filter(urbnmapr::states, state_name %in% c("South Dakota", "Nebraska", "Iowa", "Minnesota",
                                              "Missouri", "Michigan", "Indiana", "Illinois",
                                              "Wisconsin", "Kansas", "Ohio", "North Dakota"))

欢迎任何可以给我一个对象的想法，我可以用它来掩盖这个区域的光栅。

【问题讨论】：

标签： r sf

【解决方案1】：

使用urbnmapr 和sf。

在你想要的状态的 sf 对象上使用st_union() 会给你外部边界。

我认为sf 对象可用于屏蔽/裁剪光栅对象。

library(urbnmapr)
library(tidyverse)
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

states_all <- get_urbn_map(map = 'states', sf = TRUE)
my_states_vec <-  c("South Dakota", "Nebraska", "Iowa", "Minnesota",
                                           "Missouri", "Michigan", "Indiana", "Illinois",
                                           "Wisconsin", "Kansas", "Ohio", "North Dakota")
my_states <- states_all %>% 
  filter(state_name %in% my_states_vec) %>%
  st_union()

head(my_states)
#> Geometry set for 1 feature 
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -340177.8 ymin: -950695.4 xmax: 1627432 ymax: 498098.3
#> projected CRS:  US National Atlas Equal Area
#> MULTIPOLYGON (((1422979 -227154.2, 1423979 -225...

ggplot(my_states) + 
  geom_sf(fill = NA)

^{由reprex package (v0.3.0) 于 2020 年 12 月 10 日创建}

【讨论】：

【解决方案2】：

这是一些示例数据，首先是 raster/sp

library(raster)
us <- getData("GADM", level=1, country="USA")
us <- us[!(us$NAME_1 %in% c("Alaska", "Hawaii")), ]
r <- raster(us, res=1)
values(r) <- 1:ncell(r)
states <- c("South Dakota", "Nebraska", "Iowa", "Minnesota", "Missouri", "Michigan", "Indiana", "Illinois", "Wisconsin", "Kansas", "Ohio", "North Dakota")

s <- us[us$NAME_1 %in% states, ]
m <- mask(r, s)
plot(m)
lines(us)
lines(s, col="red", lwd=2)

如您所见，您不需要外边框来遮盖，但如果您想要它们，您可以使用

d <- aggregate(s)
lines(d, col="blue", lwd=3)

现在使用 urbanmapr / sf

#remotes::install_github("UrbanInstitute/urbnmapr")
library(urbnmapr)
library(sf) 
us2 <- get_urbn_map(map="states", sf = TRUE)
us2 <- sf::st_transform(us2, crs(r))
s2 <- subset(us2, state_name %in% states)
r <- raster(us, res=1)
values(r) <- 1:ncell(r)

m2 <- mask(r, s2)
plot(m2)
lines(as(us2, "Spatial"))
lines(as(s2, "Spatial"), col="red")

【讨论】：