Get basic SVG map with states & a few HUCs

.gitignore

@@ -19,6 +19,7 @@
 
 # RStudio files
 .Rproj.user/
+*.Rproj
 
 # produced vignettes
 vignettes/*.html
@@ -37,3 +38,11 @@ vignettes/*.pdf
 
 # R Environment Variables
 .Renviron
+
+# Pipeline files to ignore
+*/out/*
+*/tmp/*
+_targets/*
+
+# Exclude the empty files
+!*.empty

1_fetch.R

@@ -0,0 +1,34 @@
+# Get spatial data into sf objects
+
+source("1_fetch/src/maps_to_sf.R")
+
+p1_targets <- list(
+
+  # Albers Equal Area
+  tar_target(p1_proj_str, "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"),
+
+  tar_target(
+    p1_conus_sf,
+    generate_conus_sf(p1_proj_str)
+  ),
+
+  tar_target(
+    p1_conus_states_sf,
+    generate_conus_states_sf(p1_proj_str)
+  ),
+
+  # Get basins
+  # TODO: add more than the one IWS basin and propogate these
+  # labels through to the SVG additions.
+  tar_target(
+    p1_huc8s, c("07120001", "07120002", "07120003")
+  ),
+  tar_target(
+    p1_huc8s_sf,
+    get_huc8(id = p1_huc8s) %>% 
+      st_buffer(0) %>% 
+      st_union() %>%
+      st_make_valid() %>% 
+      st_transform(p1_proj_str)
+  )
+)

1_fetch/src/maps_to_sf.R

@@ -0,0 +1,121 @@
+# Utils for using spatial data from `maps`
+# to create state and county sf objects
+
+generate_conus_sf <- function(proj_str) {
+
+  usa_sf <- maps::map("usa", fill = TRUE, plot=FALSE) %>%
+    st_as_sf() %>% 
+    st_transform(proj_str) %>% 
+    st_buffer(0) 
+
+  return(usa_sf)
+}
+
+generate_conus_states_sf <- function(proj_str) {
+
+  usa_sf <- maps::map("usa", fill = TRUE, plot=FALSE) %>%
+    st_as_sf() %>% 
+    st_transform(proj_str) %>% 
+    st_buffer(0) 
+
+  # Need to remove islands from state outlines and then add back in 
+  # later so that they can be drawn as individual polygons. Otherwise,
+  # drawn with the state since the original state maps data only has 1
+  # ID per state.
+
+  usa_islands_sf <- usa_sf %>% filter(ID != "main")
+  usa_addl_islands_sf <- generate_addl_islands(proj_str)
+  usa_mainland_sf <- usa_sf %>% 
+    filter(ID == "main") %>% 
+    st_erase(usa_addl_islands_sf) 
+
+  # Have to manually add in CO because in `maps`, it is an incomplete
+  # polygon and gets dropped somewhere along the way.
+  co_sf <- maps::map("state", "colorado", fill = TRUE, plot=FALSE) %>%
+    st_as_sf() %>%
+    # st_buffer(0) %>% # Hmm thought it would fix the weird line but doesn't 
+    st_transform(proj_str)
+
+  maps::map("state", fill = TRUE, plot=FALSE) %>%
+    st_as_sf() %>%
+    st_transform(proj_str) %>%
+    st_buffer(0) %>% 
+    # Get rid of islands from state outline data
+    st_intersection(usa_mainland_sf) %>%
+    select(-ID.1) %>% # st_intersection artifact that is unneeded
+    # Add islands back in as separate polygons from states
+    bind_rows(usa_islands_sf) %>%
+    bind_rows(usa_addl_islands_sf) %>% 
+    st_buffer(0) %>%
+    st_cast("MULTIPOLYGON") %>% # Needed to bring back to correct type to use st_coordinates
+    rmapshaper::ms_simplify(0.5) %>%
+    bind_rows(co_sf) # bind CO after bc otherwise it gets dropped in st_buffer(0)
+
+}
+
+generate_addl_islands <- function(proj_str) {
+  # These are not called out specifically as islands in the maps::map("usa") data
+  # but cause lines to be drawn across the map if not treated separately. This creates those shapes.
+
+  # Counties to be considered as separate polygons
+
+  separate_polygons <- list(
+    `upper penninsula` = list(
+      state = "michigan",
+      counties = c(
+        "alger",
+        "baraga",
+        "chippewa",
+        "delta",
+        "dickinson",
+        "gogebic",
+        "houghton",
+        "iron",
+        "keweenaw",
+        "luce",
+        "mackinac",
+        "marquette",
+        "menominee",
+        "ontonagon",
+        "schoolcraft"
+      )),
+    `eastern shore` = list(
+      state = "virginia",
+      counties = c(
+        "accomack",
+        "northampton"
+      )),
+    # TODO: borders still slightly wonky bc it doesn't line up with counties perfectly. 
+    `nags head` = list(
+      state = "north carolina",
+      counties = c(
+        "currituck"
+      )),
+    # This + simplifying to 0.5 took care of the weird line across NY
+    `staten island` = list(
+      state = "new york",
+      counties = c(
+        "richmond"
+      )))
+
+  purrr::map(names(separate_polygons), function(nm) {
+    maps::map("county", fill = TRUE, plot=FALSE) %>%
+      sf::st_as_sf() %>%
+      st_transform(proj_str) %>% 
+      st_buffer(0) %>%
+      filter(ID %in% sprintf("%s,%s", separate_polygons[[nm]][["state"]],
+                             separate_polygons[[nm]][["counties"]])) %>% 
+      mutate(ID = nm) 
+  }) %>% 
+    bind_rows() %>% 
+    group_by(ID) %>% 
+    summarize(geom = st_union(geom))
+}
+
+# List counties to use to query `maps()`
+list_state_counties <- function(state_abbr) {
+  tolower(gsub(" County", "", countyCd$COUNTY_NAME[which(countyCd$STUSAB == state_abbr)]))
+}
+
+# Function to remove a state
+st_erase <- function(x, y) st_difference(x, st_union(st_combine(y)))

2_process.R

@@ -0,0 +1,46 @@
+# Steps for converting spatial features (sf) objects into SVG land
+
+source("2_process/src/sf_to_coords.R")
+source("2_process/src/coords_to_svg_path.R")
+
+p2_targets <- list(
+
+  tar_target(svg_width, 1000),
+  tar_target(p2_view_bbox, st_bbox(p1_conus_states_sf)),
+
+  # Prepare CONUS states for SVG
+
+  tar_target(
+    p2_conus_states_names,
+    p1_conus_states_sf %>% 
+      st_drop_geometry() %>% 
+      pull(ID)
+  ),
+
+  tar_target(
+    p2_conus_states_coords,
+    p1_conus_states_sf %>% 
+      filter(ID %in% p2_conus_states_names) %>% 
+      sf_to_coords(svg_width, view_bbox = p2_view_bbox),
+    pattern = map(p2_conus_states_names)
+  ),
+
+  tar_target(
+    p2_conus_states_paths,
+    coords_to_svg_path(p2_conus_states_coords, close_path = TRUE),
+    pattern = map(p2_conus_states_coords)
+  ),
+
+  # Prepare HUCs for SVG
+
+  tar_target(
+    p2_huc8s_coords,
+    p1_huc8s_sf %>% 
+      sf_to_coords(svg_width, view_bbox = p2_view_bbox)
+  ),
+
+  tar_target(
+    p2_huc8s_paths,
+    coords_to_svg_path(p2_huc8s_coords, close_path = TRUE)
+  )
+)

2_process/src/coords_to_svg_path.R

@@ -0,0 +1,25 @@
+coords_to_svg_path <- function(xy_df, close_path = FALSE) {
+
+  x <- xy_df$x
+  y <- xy_df$y
+
+  # Build path
+  first_pt_x <- head(x, 1)
+  first_pt_y <- head(y, 1)
+
+  all_other_pts_x <- tail(x, -1)
+  all_other_pts_y <- tail(y, -1)
+  path_ending <- ""
+  if(close_path) {
+    # Connect path to start to make polygon
+    all_other_pts_x <- c(all_other_pts_x, first_pt_x)
+    all_other_pts_y <- c(all_other_pts_y, first_pt_y)
+    path_ending <- " Z"
+  }
+
+  d <- sprintf("M%s %s %s%s", first_pt_x, first_pt_y,
+               paste0("L", all_other_pts_x, " ", 
+                      all_other_pts_y, collapse = " "),
+               path_ending)
+  return(d)
+}

2_process/src/sf_to_coords.R

@@ -0,0 +1,33 @@
+# Converting sf polygons into SVG coordinates
+# This function will work if `sf_obj` is an individual polygon
+sf_to_coords <- function(sf_obj, svg_width, view_bbox = NULL) {
+
+  coords <- st_coordinates(sf_obj)
+  x_dec <- coords[,1]
+  y_dec <- coords[,2]
+
+  # Using the whole view, figure out coordinates
+  # If view_bbox isn't provided, assume sf_obj is the whole view
+  if(is.null(view_bbox)) view_bbox <- st_bbox(sf_obj)
+
+  x_extent <- c(view_bbox$xmin, view_bbox$xmax)
+  y_extent <- c(view_bbox$ymin, view_bbox$ymax)
+
+  # Calculate aspect ratio
+  aspect_ratio <- diff(x_extent)/diff(y_extent)
+
+  # Figure out what the svg_height is based on svg_width, maintaining the aspect ratio
+  svg_height <- svg_width / aspect_ratio
+
+  # Convert longitude and latitude to SVG horizontal and vertical positions
+  # Remember that SVG vertical position has 0 on top
+  x_extent_pixels <- x_extent - view_bbox$xmin
+  y_extent_pixels <- y_extent - view_bbox$ymin
+  x_pixels <- x_dec - view_bbox$xmin # Make it so that the minimum longitude = 0 pixels
+  y_pixels <- y_dec - view_bbox$ymin # Make it so that the maximum latitude = 0
+
+  data.frame(
+    x = round(approx(x_extent_pixels, c(0, svg_width), x_pixels)$y, 6),
+    y = round(approx(y_extent_pixels, c(svg_height, 0), y_pixels)$y, 6)
+  )
+}

3_build.R

@@ -0,0 +1,57 @@
+# Build an SVG using XML components
+
+source("3_build/src/svg_xml_helpers.R")
+
+p3_targets <- list(
+
+  tar_target(
+    root_svg,
+    init_svg("3_build/tmp/root.svg", 
+             viewbox_dims = c(0, 0, svg_width=svg_width, svg_height=700)),
+    format = "file"
+  ),
+
+  # Add states
+  # TODO: groups don't seem to actually be working
+  tar_target(
+    g_conus_state_svg,
+    add_grp(out_svg = "3_build/tmp/g_conus_state.svg", 
+            in_svg = root_svg, 
+            grp_nm = 'conus-states', trans_x = 0, trans_y = 0),
+    format = "file"
+  ),
+
+  tar_target(
+    state_paths_svg, 
+    add_child_paths(out_svg = "3_build/tmp/state_paths.svg",
+                    in_svg = g_conus_state_svg,
+                    paths = p2_conus_states_paths,
+                    path_nms = sprintf('state-%s', p2_conus_states_names)),
+    format = "file"
+  ),
+
+  # Add in HUCs
+  tar_target(
+    g_huc8s_svg,
+    add_grp(out_svg = "3_build/tmp/g_huc8s.svg", 
+            in_svg = state_paths_svg, 
+            grp_nm = 'huc8s', trans_x = 0, trans_y = 0),
+    format = "file"
+  ),
+
+  tar_target(
+    huc8s_paths_svg, 
+    add_child_paths(out_svg = "3_build/tmp/huc8s_paths.svg",
+                    in_svg = g_huc8s_svg,
+                    paths = p2_huc8s_paths,
+                    path_nms = rep('huc8s', length(p2_huc8s_paths))),
+    format = "file"
+  ),
+
+  tar_target(
+    map_svg,
+    build_final_svg("3_build/out/map.svg", huc8s_paths_svg),
+    format = "file"
+  )
+
+)

3_build/src/svg_xml_helpers.R

@@ -0,0 +1,42 @@
+# Each of the steps has to read and write a file or you will get an
+# error about an invalid external pointer (this is because of how xml2
+# edits the global var, see more at https://github.com/tidyverse/rvest/issues/181)
+
+init_svg <- function(out_svg, viewbox_dims) {
+  # Create the main "parent" svg node. This is the top-level part of the svg
+  svg_root <- xml_new_root('svg', 
+               viewBox = paste(viewbox_dims, collapse=" "), 
+               preserveAspectRatio="xMidYMid meet",
+               xmlns="http://www.w3.org/2000/svg", 
+               `xmlns:xlink`="http://www.w3.org/1999/xlink") 
+  write_xml(svg_root, out_svg)
+  return(out_svg)
+}
+
+add_grp <- function(out_svg, in_svg, grp_nm, trans_x, trans_y) {
+
+  current_svg <- read_xml(in_svg)
+
+  current_svg %>% 
+    xml_add_child('g', id = grp_nm, 
+                  transform = sprintf("translate(%s %s) scale(1, 1)", trans_x, trans_y))
+
+  write_xml(current_svg, out_svg)
+  return(out_svg)
+}
+
+add_child_paths <- function(out_svg, in_svg, paths, path_nms) {
+  svg_state <- read_xml(in_svg)
+  for(i in 1:length(paths)) {
+    xml_add_child(svg_state, 'path', d = paths[i], 
+                  class = path_nms[i], 
+                  style = "stroke:#9fabb7;stroke-width:0.5;fill:none")
+  }
+  write_xml(svg_state, out_svg)
+  return(out_svg)
+}
+
+build_final_svg <- function(out_svg, in_svg) {
+  read_xml(in_svg) %>% write_xml(file = out_svg)
+  return(out_svg)
+}