如何绘制不同日期的地图变化

Question

我正在使用美国的地图形状文件。 根据 1900 年的地图，我想看看哪些部分在 1920 年改变了状态？它应该是一张由大白色和蓝色凹坑组成的地图，表明变化。

我正在使用 sf 包和 ggplot。我尝试过 st_difference 和 st_sym_difference 。但它们似乎不起作用。

形状文件下载自 https://digital.newberry.org/ahcb/downloads/united_states.html

我所做的是

states <- st_read("/US_HistStateTerr_Gen01_Shapefile/US_HistStateTerr_Gen01.shp")

states <- states %>%
  rename_all(tolower) 

state_map_1900 <- states %>% filter(ymd(19000101) %within% interval(start_date, end_date))
state_map_1920 <- states %>% filter(ymd(19201231) %within% interval(start_date, end_date))
map_change <- st_difference(state_map_1900, state_map_1920)

ggplot() +
  # Baseline map (1900)
  geom_sf(data = state_map_1900, fill = "white", color = "black") +  
  # Highlight changes
  geom_sf(data = map_change, fill = 'blue', color = NA) +  
  labs(title = "Areas that Changed States", fill = "") +
  theme_minimal()

但是整个地图必须是蓝色的。这不可能是正确的。

你能帮我一下吗？谢谢你。

Answer 1

如果您首先加入这些子集以生成具有多个几何列的

sf

对象，您可以获得成对的

st_sym_difference()

或

st_difference()

与

purrr::map2()

或

mapply()

:

library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
library(dplyr, warn.conflicts = FALSE)
library(lubridate, warn.conflicts = FALSE)
library(stringr)
library(tidyr)
library(purrr)
library(ggplot2)

# use GDAL virtual filesystems to fetch and extract shapes from online archive
vsi <- paste0(
  "/vsizip/vsicurl/",
  "https://digital.newberry.org/ahcb/downloads/gis/US_AtlasHCB_StateTerr_Gen01.zip",
  "/US_AtlasHCB_StateTerr_Gen01/US_HistStateTerr_Gen01_Shapefile"
)

# id-s can change, extract id prefix from id for joining
states <- 
  read_sf(vsi) |>  
  rename_with(tolower) |> 
  mutate(id_pref = str_split_i(id, "_", 1), .before = id) |> 
  select(id_num:end_date, area_sqmi, geometry)
print(states, n = 4)
#> Simple feature collection with 220 features and 8 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -179.1474 ymin: 18.91228 xmax: -66.94993 ymax: 71.38961
#> Geodetic CRS:  WGS 84
#> # A tibble: 220 × 9
#>   id_num name              id_pref id    version start_date end_date   area_sqmi
#>    <int> <chr>             <chr>   <chr>   <int> <date>     <date>         <dbl>
#> 1      1 Alaska Department ak      ak_d…       1 1867-03-30 1884-05-16    575301
#> 2      2 Alaska District   ak      ak_d…       1 1884-05-17 1912-08-23    575301
#> 3      3 Alaska            ak      ak_s…       1 1959-01-03 2000-12-31    575301
#> 4      4 Alaska Territory  ak      ak_t…       1 1912-08-24 1959-01-02    575301
#> # ℹ 216 more rows
#> # ℹ 1 more variable: geometry <MULTIPOLYGON [°]>

# full join subsets by id prefix, as_tibble(y) as full_join of 2 sf objects will fail,
# generates 2 geometry columns (geometry.00, geometry.20);
# iterate though geometry pairs to get pair-wise differences as 3rd geometry column, 
# add area from (sym) diff
state_map_1900_1920 <- 
  full_join(
    states |>             filter(ymd(19000101) %within% interval(start_date, end_date)),
    as_tibble(states) |>  filter(ymd(19201231) %within% interval(start_date, end_date)),
    by = "id_pref", suffix = c(".00", ".20")) |> 
  mutate(geom_diff = map2(geometry.00, geometry.20, st_sym_difference) |> st_sfc(crs = st_crs(geometry.00)),
         diff_area = st_area(geom_diff) |> units::set_units(mi^2))

glimpse(state_map_1900_1920)
#> Rows: 52
#> Columns: 19
#> $ id_num.00     <int> 2, 7, 11, 17, 18, 20, 25, 27, 28, 39, 44, 46, 49, 51, 55…
#> $ name.00       <chr> "Alaska District", "Alabama", "Arkansas", "Arizona Terri…
#> $ id_pref       <chr> "ak", "al", "ar", "az", "ca", "co", "ct", "dc", "de", "f…
#> $ id.00         <chr> "ak_district", "al_state", "ar_state", "az_terr", "ca_st…
#> $ version.00    <int> 1, 3, 2, 2, 1, 1, 4, 2, 1, 1, 3, 1, 1, 1, 1, 1, 4, 1, 2,…
#> $ start_date.00 <date> 1884-05-17, 1820-12-19, 1840-05-21, 1866-05-05, 1850-09…
#> $ end_date.00   <date> 1912-08-23, 2000-12-31, 2000-12-31, 1912-02-13, 1959-12…
#> $ area_sqmi.00  <dbl> 575301, 51656, 53179, 113999, 158097, 104093, 4975, 68, …
#> $ geometry.00   <MULTIPOLYGON [°]> MULTIPOLYGON (((-179.0575 5..., MULTIPOLYGO…
#> $ id_num.20     <int> 4, 7, 11, 15, 18, 20, 25, 27, 28, 39, 44, 48, 49, 51, 55…
#> $ name.20       <chr> "Alaska Territory", "Alabama", "Arkansas", "Arizona", "C…
#> $ id.20         <chr> "ak_terr", "al_state", "ar_state", "az_state", "ca_state…
#> $ version.20    <int> 1, 3, 2, 1, 1, 1, 4, 2, 1, 1, 3, 1, 1, 1, 1, 1, NA, 1, 2…
#> $ start_date.20 <date> 1912-08-24, 1820-12-19, 1840-05-21, 1912-02-14, 1850-09…
#> $ end_date.20   <date> 1959-01-02, 2000-12-31, 2000-12-31, 2000-12-31, 1959-12…
#> $ area_sqmi.20  <dbl> 575301, 51656, 53179, 113999, 158097, 104093, 4975, 68, …
#> $ geometry.20   <MULTIPOLYGON [°]> MULTIPOLYGON (((-179.0575 5..., MULTIPOLYGO…
#> $ geom_diff     <GEOMETRY [°]> GEOMETRYCOLLECTION EMPTY, GEOMETRYCOLLECTION EM…
#> $ diff_area     [mi^2] 0.0000 [mi^2], 0.0000 [mi^2], 0.0000 [mi^2], 0.0000 [mi…

变化：

state_map_1900_1920 |> 
  filter(!st_is_empty(geom_diff)) |> 
  st_drop_geometry() |>
  select(id.00, name.00, end_date.00, id.20, name.20, area_sqmi.00, area_sqmi.20, diff_area)
#> # A tibble: 5 × 8
#>   id.00    name.00 end_date.00 id.20 name.20 area_sqmi.00 area_sqmi.20 diff_area
#>   <chr>    <chr>   <date>      <chr> <chr>          <dbl>        <dbl>    [mi^2]
#> 1 it_indi… Indian… 1907-11-15  <NA>  <NA>           31075           NA  31074.  
#> 2 mi_state Michig… 1908-12-31  mi_s… Michig…        58121        58371    246.  
#> 3 ne_state Nebras… 1902-11-03  ne_s… Nebras…        77361        77353      7.56
#> 4 ok_terr  Oklaho… 1907-11-15  ok_s… Oklaho…        38824        69899  31074.  
#> 5 sd_state South … 1902-11-03  sd_s… South …        77107        77115      7.56

ggplot(state_map_1900_1920) +
  geom_sf(fill = "white", color = "black", linewidth = .05) +
  geom_sf(aes(geometry = geom_diff), fill = 'blue', color = "blue") +
  coord_sf(crs = 5070) +
  theme_minimal()

^{创建于 2024-09-23，使用 reprex v2.1.1}