2022-11-01_30daymapchallenge_points.R

library(getdata)
library(ggplot2)
library(mapmaryland)
library(maplayer)
library(dplyr)

dead_zones <-
  get_esri_data(
    url = "https://geodata.md.gov/imap/rest/services/Environment/MD_ChesapeakeBayDeadZones/FeatureServer/0",
    crs = 3857
  )

streams <-
  get_water_data(
    location = mapbaltimore::baltimore_msa_counties,
    crs = 3857
  )

ponds <-
  get_esri_data(
    url = "https://geodata.md.gov/imap/rest/services/Agriculture/MD_NutrientManagementSetbacksFromWaterways/MapServer/6",
    location = mapbaltimore::baltimore_msa_counties,
    crs = 3857
  )

msa_area <-
  mapbaltimore::baltimore_msa_counties %>% 
  sfext::st_bbox_ext(dist = 3, unit = "mi", crs = 3857, class = "sf")

area_circle <-
  dead_zones %>% 
  sf::st_crop(msa_area) %>% 
  sf::st_union() %>% 
  sf::st_convex_hull() %>% 
  sf::st_inscribed_circle(dTolerance = 1000)

streams_trim <-
  sfext::st_trim(streams, area_circle)

ponds_trim <-
  sfext::st_trim(ponds, area_circle)

ggplot() +
  layer_location_data(
    data = mapmaryland::md_counties_detailed,
    size = 0.1, alpha = 0.2,
    color = "gray85"
  ) +
  layer_location_data(
    data = streams_trim,
    color = "steelblue2", size = 0.4, alpha = 0.2
  ) +
  layer_location_data(
    data = ponds_trim,
    fill = "steelblue",
    color = NA,
    alpha = 0.25
  ) +
  layer_location_data(
    data = dead_zones,
    fn = ~ .x %>% filter(minimum < 5),
    aes(color = minimum),
    size = 0.29, alpha = 0.8, shape = 20
  ) +
  theme_void() +
  hrbrthemes::theme_ipsum_rc() +
  layer_neatline(
    data = msa_circle,
    crs = 3857, dist = 1, unit = "mi",
    color = "gray60", size = 0.75,
    asp = sfext::get_paper("Instagram story")$asp
  ) +
  theme_legend("bottomleft", bgcolor = NA) +
  labs(
    title = "Dissolved oxygen levels of 3-5 mg/l are stressful to\nsome organisms and 2-3 mg/l are borderline lethal",
    color = "Dissolved\noxygen (mg/l)",
    caption = "Data: Chesapeake Bay Dead Zones (2019), Maryland DNR/MD iMap"
  ) +
  scale_color_distiller(type = "seq", palette = "YlOrRd", direction = -1) +
  theme(
    legend.background = element_blank(),
    plot.title = element_text(size = 12),
    panel.grid = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
    )

ggsave_ext(
  filename = "chesapeake_bay_dead_zones.png",
  paper = "Instagram story",
  scale = 1.25
)

2022-11-02_30daymapchallenge_lines.R

library(dplyr)
library(tidyverse)

ports <-
  tibble::tribble(
             ~Port,        ~lat,         ~lon,
       "Baltimore", 39.25659062, -76.54659451,
         "Norfolk", 36.91749489, -76.29947623,
        "New York", 40.68478577, -74.00619792,
    "Philadelphia", 39.98048588, -75.08957318,
        "Delaware", 39.71528893, -75.52992711
    )


data <-
  tibble::tribble(
          ~State, ~Baltimore, ~Norfolk, ~`New York`, ~Philadelphia, ~Delaware,  ~Total,
      "Maryland",     20868L,     580L,     771L,          392L,       25L,  22636L,
  "Pennsylvania",      3224L,     219L,   10727L,        11963L,       71L,  26204L,
    "New Jersey",      3152L,     182L,   47529L,         5025L,       46L,  55934L,
      "Michigan",      1847L,     526L,    1727L,          228L,        3L,   4331L,
          "Ohio",      1276L,     893L,    3061L,          785L,        3L,   6019L,
      "New York",       980L,     262L,   23941L,          936L,       75L,  26194L,
      "Kentucky",       826L,    1759L,     572L,           25L,        0L,   3183L,
      "Illinois",       757L,     443L,    1664L,          954L,        2L,   3819L,
      "Virginia",       631L,   14475L,     525L,          146L,        6L,  15783L,
       "Indiana",       289L,     296L,     999L,          108L,        1L,   1692L,
      "Delaware",       122L,      11L,     406L,          145L,     1211L,   1895L,
     "Tennessee",       106L,     132L,     175L,            2L,        0L,    415L,
         "Total",     34078L,   19778L,   92097L,        20710L,     1443L, 168107L
  )


data <- data %>% 
  pivot_longer(
    cols = -State,
    names_to = "Port"
  )

port_totals <-
  data %>% 
  filter(State == "Total", Port != "Total") %>% 
  left_join(ports, by = "Port") %>% 
  sf::st_as_sf(
    coords = c("lon", "lat"),
    crs = 4326
  )

state_totals <-
  data %>% 
  filter(Port == "Total", State != "Total")

us_states <- tigris::states(cb = TRUE)

us_states_pt <- sf::st_centroid(us_states) %>% 
  rename(State = NAME)

state_data <-
  data %>% 
  filter(Port != "Total", State != "Total") %>% 
  left_join(us_states_pt, by = "State") %>% 
  sf::st_as_sf()

port_data <-
  data %>% 
  filter(Port != "Total", State != "Total") %>% 
  left_join(ports, by = "Port") %>% 
  sf::st_as_sf(
    coords = c("lon", "lat"),
    crs = 4326
  )

lines <- 
  map_dfr(
   seq_along(port_data$geometry),
    ~ sfext::as_sf(
       sf::st_cast(
         c(
           state_data$geometry[[.x]],
           port_data$geometry[[.x]]
           # )
         ),
         to = "LINESTRING"
       ),
       crs = 4326
     )   
  )

data <- data %>% 
  filter(Port != "Total", State != "Total")
  
sf::st_geometry(data) <- as_sfc(lines)

ggplot() +
  layer_mapbox(
    data = data,
    style_url = "mapbox://styles/elipousson-baltimorecity/cla0jnoxl000414mjexsdhgwf",
    diag_ratio = 0.05,
    asp = 6 / 4
    ) +
  geom_sf(
    data = data,
    aes(lwd = value, color = Port),
    alpha = 0.2
    ) +
  geom_sf(
    data = port_totals,
    aes(size = value, color = Port),
    alpha = 0.8
  ) +
  # geomtextpath::geom_labelsf(
  #   data = data %>% filter(Port == "Baltimore"),
  #   aes(label = value,
  #       lwd = value, color = Port),
  # ) +
  scale_color_brewer(type = "qual", palette = "Dark2") +
  scale_size_continuous(labels = scales::label_dollar()) +
 # gghighlight::gghighlight(Port == "Baltimore") +
  theme_legend("topleft") +
  guides(color = "none") +
  labs(
    lwd = "Millions ($)"
  )