如何在R制作的Map上加粗/加粗州边界,同时包括县边界?

au9on6nz  于 2023-06-03  发布在  其他
关注(0)|答案(1)|浏览(124)

我正在试图绘制美国所有50个州及其各自的县的Map。但我特别想加厚/大胆的美国国家边界,同时保持县边界薄。
这是我试图Map我的数据的代码。请注意,我不得不使用不同的代码来修改阿拉斯加和夏威夷的位置。
首先,这是绘制县边界和结果Map的代码。

map_sf <- tigris::counties(cb = T, class = 'sf')
# removed US territories
map_sf <- map_sf %>% filter(!STATEFP %in% c('60', '66', '69', '72', '78'))
# CRS code to modify Alaska and Hawaii location
crs_lambert <- "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"
map_sf <- map_sf %>%
  st_transform(crs = crs_lambert)
alaska <- map_sf %>% filter(STATE_NAME %in% 'Alaska')
alaska_g <- st_geometry(alaska)
alaska_centroid <- st_centroid(st_union(alaska_g))
rot <- function(a) matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)
alaska_trans <- (alaska_g - alaska_centroid) * rot(-39 * pi/180) / 2.3 + alaska_centroid + c(1000000, -5000000)
alaska <- alaska %>% st_set_geometry(alaska_trans) %>% st_set_crs(st_crs(df))
hawaii <- map_sf %>% filter(STATE_NAME %in% 'Hawaii')
hawaii_g <- st_geometry(hawaii)
hawaii_centroid <- st_centroid(st_union(hawaii_g))
hawaii_trans <- (hawaii_g - hawaii_centroid) * rot(-35 * pi/180) + hawaii_centroid + c(5200000, -1400000)
hawaii <- hawaii %>% st_set_geometry(hawaii_trans) %>% st_set_crs(st_crs(map_sf))
map_sf <- map_sf %>%
  filter(!STATE_NAME %in% c('Alaska', 'Hawaii')) %>%
  rbind(alaska) %>%
  rbind(hawaii)
map_sf <- map_sf %>% rename(county = "NAMELSAD", state = "STUSPS") %>% select(county, state, geometry)
# changed projection to longlat
crs_lambert <- "+proj=longlat +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"
map_sf <- map_sf %>%
  st_transform(crs = crs_lambert)

county map of US
这是我尝试为美国各州制作的Map,再次修改了阿拉斯加和夏威夷的位置。

map_us <- tigris::states(class = 'sf')
map_us <- map_us %>% filter(!STATEFP %in% c('60', '66', '69', '72', '78')) 
crs_lambert <- "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"
map_us <- map_us %>%
  st_transform(crs = crs_lambert)
alaska <- map_us %>% filter(NAME %in% 'Alaska')
alaska_g <- st_geometry(alaska)
alaska_centroid <- st_centroid(st_union(alaska_g))
rot <- function(a) matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)
alaska_trans <- (alaska_g - alaska_centroid) * rot(-39 * pi/180) / 2.3 + alaska_centroid + c(1000000, -5000000)
alaska <- alaska %>% st_set_geometry(alaska_trans) %>% st_set_crs(st_crs(df))
hawaii <- map_us %>% filter(NAME %in% 'Hawaii')
hawaii_g <- st_geometry(hawaii)
hawaii_centroid <- st_centroid(st_union(hawaii_g))
hawaii_trans <- (hawaii_g - hawaii_centroid) * rot(-35 * pi/180) + hawaii_centroid + c(5200000, -1400000)
hawaii <- hawaii %>% st_set_geometry(hawaii_trans) %>% st_set_crs(st_crs(map_us))
map_us <- map_us %>%
  filter(!NAME %in% c('Alaska', 'Hawaii')) %>%
  rbind(alaska) %>%
  rbind(hawaii)
map_us <- map_us %>% rename(state = "STUSPS") %>% select(state, geometry)
crs_lambert <- "+proj=longlat +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"
map_us <- map_us %>%
  st_transform(crs = crs_lambert)
ggplot() +
  geom_sf(data = map_us, linewidth = 1)

US state map borders
这是当我尝试叠加两张Map时的结果。

ggplot() +
  geom_sf(data = map_us, linewidth = 1) +
  geom_sf(data = map_sf, fill = NA)

overlay of county and state map
正如您所看到的,州边界与县Map上的州边界并不完全一致,尤其是阿拉斯加和夏威夷。但我确实对两张Map使用了完全相同的CRS修改。
有没有其他方法可以增加州边界的厚度,使其与县Map上的边界完全一致?

kgsdhlau

kgsdhlau1#

我认为你是过于复杂的代码(虽然这是真的,练习本身并不容易实现)。
首先,让我们提取数据(我清理了一点代码):

library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE
library(ggplot2)

## County map
map_sf <- tigris::counties(cb = T, class = "sf")

# removed US territories
map_sf <- map_sf %>% filter(!STATEFP %in% c("60", "66", "69", "72", "78"))

# CRS code to modify Alaska and Hawaii location
crs_lambert <- "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"

map_sf <- map_sf %>%
  st_transform(crs_lambert)

# Function to rotate
rot <- function(a) matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)

# Alaska
alaska <- map_sf %>% filter(STATE_NAME %in% "Alaska")
alaska_g <- st_geometry(alaska)
alaska_centroid <- st_centroid(st_union(alaska_g))
alaska_trans <- (alaska_g - alaska_centroid) * rot(-39 * pi / 180) / 2.3 + alaska_centroid + c(1000000, -5000000)

alaska <- alaska %>%
  st_set_geometry(alaska_trans) %>%
  st_set_crs(st_crs(map_sf))

# Hawaii
hawaii <- map_sf %>% filter(STATE_NAME %in% "Hawaii")
hawaii_g <- st_geometry(hawaii)
hawaii_centroid <- st_centroid(st_union(hawaii_g))
hawaii_trans <- (hawaii_g - hawaii_centroid) * rot(-35 * pi / 180) + hawaii_centroid + c(5200000, -1400000)
hawaii <- hawaii %>%
  st_set_geometry(hawaii_trans) %>%
  st_set_crs(st_crs(map_sf))

在这一点上,你已经拥有了你所需要的一切。

# Union all

map_sf <- map_sf %>%
  filter(!STATE_NAME %in% c("Alaska", "Hawaii")) %>%
  rbind(alaska) %>%
  rbind(hawaii) %>%
  rename(county = "NAMELSAD", state = "STUSPS") %>%
  select(county, state)

对象同时具有县和州。您可以在此处生成(即聚合)属于同一个州的所有县的形状,实际上,您将生成一个州的shapefile**,它与县**完全对齐:

# Create here the state map

map_us <- map_sf %>%
  group_by(state) %>%
  summarise()

现在最后的操纵和阴谋。边界应该完全对齐(因为州是从县的边界生成的),并且您不需要重复所有州级别信息的代码。

# changed projection to longlat
crs_lonlat <- "+proj=longlat +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"

map_sf <- map_sf %>%
  st_transform(crs_lonlat)

map_us <- map_us %>%
  st_transform(crs_lonlat)

ggplot() +
  geom_sf(data = map_us, linewidth = 1) +
  geom_sf(data = map_sf, fill = NA)

创建于2023-06-01使用reprex v2.0.2
希望这个能帮上忙

相关问题