global map

gage_timeseries_gif/_targets.R

@@ -47,7 +47,10 @@ usgs_logo <- magick::image_read('usgs_logo_grey.png') |>
   magick::image_resize('x80') |>
   magick::image_colorize(100, grey_light)
 
+# projected maps
 proj.string <- "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"
+# global map
+ortho_crs <- "+proj=ortho +lat_0=40 +lon_0=-100 +x_0=0 +y_0=0 +ellps=WGS84 +no_defs"
 
 p1_fetch_targets <- list(
   # Output of `national-flow-observations`
@@ -83,6 +86,12 @@ p1_fetch_targets <- list(
   tar_target(
     state_map_PR,
     extract_states(area_name = "PR")
+  ),
+
+  # prep globe for global layout
+  tar_target(
+    world_map_sf,
+    rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")
   )
 )
 
@@ -114,6 +123,16 @@ p2_process_targets <- list(
     site_map_PR,
     extract_sites(area_name = "PR", 
                   gage_info = gage_info)
+  ),
+
+  # merge all sites together for global view map
+  tar_target(
+    sites_all_global,
+    harmonize_sites(in_CONUS = site_map_CONUS,
+                    in_HI = site_map_HI,
+                    in_PR = site_map_PR,
+                    in_AK = site_map_AK,
+                    crs = ortho_crs)
   )
 )
 
@@ -224,6 +243,27 @@ downstream_targets <- list(
   )
 )
 
+# # # # # # # # # # # # # # # # # # # # # # #  
+#
+# 
+#       GLOBAL MAP of STREAMGAGES
+# 
+#
+global_targets <- list(
+  tar_target(
+    global_png,
+    compose_global_map(
+      gage_melt = gage_melt,
+      sites_in_sf = sites_all_global,
+      globe_in_sf = world_map_sf,
+      focal_year = year_for_plot,
+      crs = ortho_crs,
+      png_out = "out/global_gages_map.png"
+    ),
+    format = "file"
+  )
+)
+
 # # # # # # # # # # # # # # # # # # # # # # #  
 #
 # 
@@ -349,6 +389,7 @@ c(p1_fetch_targets,
   p3_viz_split_targets, 
   p3_viz_combine_targets, 
   downstream_targets,
+  global_targets,
   gw_targets)
 
 

gage_timeseries_gif/src/gw_utils.R

@@ -45,17 +45,15 @@ extract_gw_sites <- function(area_name, site_info){
       filter(! state %in% names(crs_map)) |>
       points_sp() |> # (this function is in map_utils.R)
       sf::st_as_sf() |>
-      sf::st_transform(crs = 5070) |>
-      left_join(site_info, by = "site_no")
+      sf::st_transform(crs = 5070) 
 
   } else if(area_name != "CONUS") {
 
     sites_out <- site_info %>% 
       filter(state == area_name) |>
       points_sp() |>
       sf::st_as_sf() |>
-      sf::st_transform(crs = crs_map[[area_name]]) |>
-      left_join(site_info, by = "site_no")
+      sf::st_transform(crs = crs_map[[area_name]]) 
   }
 
   return(sites_out)

gage_timeseries_gif/src/map_utils.R

@@ -3,12 +3,13 @@
 #' 
 #' @param locations coordinates of the spatial points for mapping
 #' 
-points_sp <- function(locations){
-  points <- cbind(locations$dec_long_va, locations$dec_lat_va) 
-  points_sp_obj <- sp::SpatialPoints(points, sp::CRS("+proj=longlat +datum=WGS84")) 
-  points_transform <- sp::spTransform(points_sp_obj, proj.string) %>% 
-    sp::SpatialPointsDataFrame(data = locations[c('site_no')])
-  return(points_transform)
+points_sp <- function(locations, crs_out){
+  sf::st_as_sf(
+    locations,
+    coords = c("dec_long_va", "dec_lat_va"),
+    crs = 4326  # WGS84
+  ) #|>
+    #sf::st_transform(crs = crs_out)
 }
 
 
@@ -22,12 +23,12 @@ extract_states <- function(area_name){
   crs_map <- c(AK = 3338, HI = 6633, PR = 4139)
 
   if(area_name == "CONUS") {
-    state_map <- tigris::states(cb = TRUE) %>%
-      filter(!(STUSPS %in% c('AS', 'MP','GU', 'VI', 'HI', 'PR', 'AK'))) %>%
+    state_map <- tigris::states(cb = TRUE) |>
+      filter(!(STUSPS %in% c('AS', 'MP','GU', 'VI', 'HI', 'PR', 'AK'))) |>
       sf::st_transform(crs = 5070) 
   } else {
-    state_map <- tigris::states(cb = TRUE) %>%
-      filter(STUSPS == area_name) %>%
+    state_map <- tigris::states(cb = TRUE) |>
+      filter(STUSPS == area_name) |>
       sf::st_transform(crs = crs_map[area_name]) 
   }
 
@@ -39,18 +40,22 @@ extract_states <- function(area_name){
 #' @param gage_data data frame of gages read from .rds file
 #' 
 fetch_gage_info <- function(gage_data){
-  gages_info <- gage_data %>%
-    pull(site) %>% 
-    dataRetrieval::readNWISsite() %>% 
-    group_by(site_no) %>% 
-    # parse huc_cd to 2 digits, and rename to huc to stay consistent
-    # 3/19/2020 - discovered that some sites have more than one HUC code, not sure how that
-    # is possible, but filtering to the last one (most recent) for now
-    #   E.g. site number `11434500` had huc_cd=="18020129" for year == 2014 and then 
-    #     huc_cd `16050101` for year == 2016
-    summarize(huc = stringr::str_sub(tail(unique(huc_cd),1), 1L, 2L), 
-              # huc = paste(stringr::str_sub(unique(huc_cd)[[1]], 1L, 2L), collapse = "|"),
-              dec_lat_va = mean(dec_lat_va), dec_long_va = mean(dec_long_va)) %>% 
+  site_ids <- gage_data |>
+    pull(site) |>
+    unique()
+
+  # need to chunk sites to pull info from dataRetrieval
+  site_chunks <- split(site_ids, ceiling(seq_along(site_ids) / 100))
+
+  gages_info <- site_chunks |>
+    map_dfr(possibly(readNWISsite, otherwise = NULL)) |>
+    group_by(site_no) |>
+    summarize(
+      huc = stringr::str_sub(tail(unique(huc_cd), 1), 1L, 2L),
+      dec_lat_va = mean(dec_lat_va, na.rm = TRUE),
+      dec_long_va = mean(dec_long_va, na.rm = TRUE),
+      .groups = "drop"
+    ) |> 
     filter(dec_long_va < -65.4,
            dec_lat_va > 0)  # remove US virgin Islands and other things we won't plot
 
@@ -71,15 +76,16 @@ extract_sites <- function(area_name, gage_info){
 
   if(area_name == "CONUS") {
 
-    sites_out <- gage_info %>% 
-      filter(!huc %in% huc_map) %>% 
+    sites_out <- gage_info |> 
+      filter(!huc %in% huc_map) |> 
+      filter(!is.na(dec_lat_va), !is.na(dec_long_va)) |>  # exclude NA coordinates
       points_sp() |>
       sf::st_as_sf() |>
       sf::st_transform(crs = 5070)
 
   } else if(area_name != "CONUS") {
 
-    sites_out <- gage_info %>% 
+    sites_out <- gage_info |> 
       filter(huc == huc_map[[area_name]]) |>
       points_sp() |>
       sf::st_as_sf() |>
@@ -89,4 +95,34 @@ extract_sites <- function(area_name, gage_info){
   return(sites_out)
 }
 
-
+#' 
+#' Bind all sites together, reprojecting them to ortho/global projection
+#' 
+#' @param in_CONUS sf object with sites in CONUS
+#' @param in_HI sf object with sites in Hawaii
+#' @param in_PR sf object with sites in PR
+#' @param in_AK sf object with sites in Alaska
+#' @param crs the ortho projection as string
+#' 
+harmonize_sites <- function(in_CONUS, in_HI, in_PR, in_AK, crs){
+
+  temp_conus <- in_CONUS |>
+    sf::st_transform(crs = crs) |>
+    mutate(location = "CONUS")
+
+  temp_AK <- in_AK |>
+    sf::st_transform(crs = crs) |>
+    mutate(location = "AK")
+
+  temp_HI <- in_HI |>
+    sf::st_transform(crs = crs) |>
+    mutate(location = "HI")
+
+  temp_PR <- in_PR |>
+    sf::st_transform(crs = crs) |>
+    mutate(location = "PR")
+
+  out_sf <- bind_rows(temp_conus, temp_AK, temp_HI, temp_PR)
+
+  return(out_sf)
+}

gage_timeseries_gif/src/plot_utils.R

@@ -36,7 +36,7 @@ plot_gage_age <- function(gage_melt, yr, font_fam){
       axis.text.y.left = element_blank(),
       axis.text.x = element_text(vjust = 0, margin = margin(r = 0, t = -2))
     )
-    
+
 
 }
 
@@ -105,7 +105,7 @@ plot_gage_map <- function(gage_melt, yr, site_map, state_map){
   yr_gages <- gage_melt |>
     filter(year == as.numeric(yr)) |>
     distinct(site)
-                                  
+
   # sites active in a given year
   gages_active <- site_map |>
     dplyr::filter(site_no %in% yr_gages$site)
@@ -145,7 +145,7 @@ compose_chart <- function(bar_chart,
                           standalone_logic,
                           gage_melt){
 
-
+  
   # summary stats for standalone viz
   average_age <- gage_melt |>
     filter(year == yr) |>
@@ -155,7 +155,7 @@ compose_chart <- function(bar_chart,
 
   # base composition for both standalone and gif
   base_plot <- ggdraw(xlim = c(0, 1), 
-         ylim = c(0,1)) +
+                      ylim = c(0,1)) +
     # create background canvas
     draw_grob(grid::rectGrob(
       x = 0, y = 0, 
@@ -173,7 +173,7 @@ compose_chart <- function(bar_chart,
       x = -0.02,
       y = 0.12,
       width = 1.03
-      )+
+    )+
     draw_text("Alaska",
               x = 0.17, y = 0.27,
               color = grey_dark,
@@ -255,7 +255,7 @@ compose_chart <- function(bar_chart,
                 y = 0.43,
                 hjust = 1)
   }
-
+  
 
   ggsave(png_out, 
          width = 5, height = 5, dpi = 300, units = 'in')
@@ -337,3 +337,77 @@ optimize_gif <- function(out_file, frame_delay_cs) {
   return(out_file)
 }
 
+#' Create global map
+#' 
+#' @param gage_melt gages by year 
+#' @param sites_in_sf global sf of all streamgages
+#' @param globe_in_sf time spent on each frame
+#' @param focal_year year for plotting
+#' @param png_out chr string of save location for final png
+#' @param crs crs for global map layout
+#' 
+compose_global_map <- function(gage_melt, 
+                               sites_in_sf, 
+                               globe_in_sf,
+                               focal_year,
+                               crs,
+                               ortho_crs,
+                               png_out){
+
+  # filter gage data to given year
+  yr_gages <- gage_melt |>
+    filter(year == as.numeric(focal_year)) |>
+    distinct(site)
+
+  # sites active in a given year
+  gages_active <- sites_in_sf |>
+    dplyr::filter(site_no %in% yr_gages$site)
+
+  # generate lat/lon gridlines
+  graticules <- sf::st_graticule(lat = seq(-90, 90, by = 15),
+                                 lon = seq(-180, 180, by = 15))
+
+  # transform graticules
+  graticules_ortho <- sf::st_transform(graticules, crs = crs) 
+
+  # create clipping circle in projected space
+  circle <- sf::st_sfc(sf::st_point(c(0, 0)), crs = crs) |>
+    sf::st_buffer(dist = 6371000)  # match Earth's radius in meters
+
+  # clip graticules
+  graticules_clipped <- sf::st_intersection(graticules_ortho, circle) 
+
+  # transform world to same projection
+  world_ortho <- sf::st_transform(globe_in_sf, crs = crs)
+
+  # background for behind globe
+  globe_background <- grid::grid.circle(x=0.5, y=0.5, r=0.5, default.units="npc", name=NULL,
+                                        gp=grid::gpar(fill = "#CDE0EA"), draw=TRUE, vp=NULL) #A8BBCD
+
+  # plot with orthographic projection
+  globe_map <- ggplot(world_ortho) +
+    geom_sf(data = graticules_clipped, 
+            linewidth = 5) +
+    #geom_sf(data = graticules_clipped, color = grey_dark, linewidth = 0.2, fill = NA) +  # lat/lon lines
+    geom_sf(fill = grey_dark, color = "#CDE0EA", alpha = 0.5, linewidth = 0.5) +
+    geom_sf(data = gages_active, color = blue_color, size = 1, alpha = 0.7) +
+    #coord_sf(crs = "+proj=ortho +lat_0=43 +lon_0=290") +  # adjust lat/lon center
+    theme_void()+
+    theme(plot.background = element_rect(fill = NA, color = NA),
+          panel.background = element_rect(fill = NA, color = NA))
+
+  composition <- ggdraw(ylim = c(0, 16), # 0-16 scale makes it easy to place viz items on canvas
+                        xlim = c(0, 16))  +
+    draw_plot(globe_background,
+              x = 0.8, y = -0 + 0.68, width = 16*0.912, height = 16*0.912) +
+    draw_plot(globe_map, 
+              x = 0.1, y = -0, width = 16, height = 16) +
+    draw_text("Alaska",
+              x = 0.17, y = 0.27,
+              color = grey_dark,
+              family = font_fam, size = 24)
+
+
+  ggsave(png_out, width = 16, height = 16, units = "in", dpi = 300, bg = "white")
+
+}