Data Visualization - Geospatial Data
Luke Ghallahorne
2024-12-09
The Data
For this assignment, I figured it was a good opportunity to highlight a mapping loop that I wrote this Fall for my project. This code pulls my larval tracking data and overlays it on a map. I copied and re-annotated the code below. A lot of it is pulling from files I have produced using a larval tracking model developed by Jake Lawlor and myself.
# library
library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.4.3## Warning: package 'ggplot2' was built under R version 4.4.3## Warning: package 'purrr' was built under R version 4.4.3## Warning: package 'lubridate' was built under R version 4.4.3## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(sf)## Warning: package 'sf' was built under R version 4.4.3## Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE# the here package makes calling files easier when using paste() or paste0()
library(here)## here() starts at C:/Users/lukeg/OneDrive/Desktop/WWU MESP/ESCI599_DataSeminar/ESCI_599_Data_Viz/9_Geospatial_Datalibrary(PNWColors)
library(ggrepel)## Warning: package 'ggrepel' was built under R version 4.4.3# load necessary objects
## this is the starting dataframe with my different release locations
load(here("data","releasesites_WADFW.RDATA"))
## this is the coordinate locations for the PNNL output data
### basically, these are where the currents and temperatures, etc., are stored that the larval dispersal model uses
### loaded here for use in the basemap
load(here("data","WA_nodes_coords.Rdata"))
load(here("data","WA_neles_coords.Rdata"))Step 1: Create Basemap
# WA coast Shorezone detailed coastline
WA_coastline_2 <- read_sf("data/ShoreZone_Inventory_-_Shoreline_type.shp")
# transform reference to NAD83 zone 10
WA_coastline_2 <- st_transform(WA_coastline_2, crs = st_crs(26910))
# find trimming coordinates
bm_xmin <- min(WA_nodes_coords$x)
bm_xmax <- max(WA_nodes_coords$x)
bm_ymin <- min(WA_nodes_coords$y)
bm_ymax <- max(WA_nodes_coords$y)
# make bounding box to limit coast shapefile
coast_box <- st_polygon(
list(
cbind(c(bm_xmin, bm_xmax, bm_xmax, bm_xmin, bm_xmin),
c(bm_ymin, bm_ymin, bm_ymax, bm_ymax, bm_ymin))
)
)
coast_poly <- st_sfc(coast_box, crs = st_crs(26910))
# Trim coastline polygon to WA_nodes_coords limits
WA_coastline_3 <- st_intersection(WA_coastline_2, coast_poly)## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries# build basemap with WA coastline and PNNL SSM nodes
basemap_nodes <- ggplot() +
geom_point(data = WA_nodes_coords, aes(x = x, y = y),
color="turquoise3", size = 0.3, alpha = 0.5) +
geom_sf(data = WA_coastline_2, alpha = 0.6) +
labs(title = "SSM High-Resolution Nodes (2014)") +
theme(plot.title = element_text(hjust = 0.5))
basemap_nodes
Step 2: Map Loop
This section is set up as a loop so I can easily change one thing and pull data for my different release locations. For these purposes, I’m going to stick with two release sites: Padilla Bay and Kilisut Harbor, because they look really nice.
# set year, as it will eventually be from 2014, 2017, or 2095
year = 14
# Start loop per release site
## KH = Kilisut Harbor
## PB = Padilla Bay
for(site in c("KH", "PB")) {
# set release hours (part of the dispersal model)
# high and low tides during neap of early July
release_hours <- c(477, 485)
# set site name and ID
site_id <- site
site_name <- as.character(filter(releasesites_WADFW, id == site_id)$site)
# scaling factor for modifying basemap axis limits
map_scale <- 1000
# setup empty vectors
map_df <- vector()
map_set <- vector()
map_set_mids <- vector()
map_set_tracks <- vector()
map_ns <- vector()
# progress bar for inner loop
pb <- txtProgressBar(min = 1, max = length(release_hours) * 3, style = 3)
pb_i <- 0
# Inner Loop 1: Combine Larval track data across release sites and behaviors
for(release in release_hours) {
for (behavior in c("none", "photo", "onto")) {
# Load disperal model output
name_long <- paste0("larvatrack_",year,"_",release,
"_",behavior,"_",site,"_df")
load(here("data", paste0(name_long,".Rdata")))
load(here("data", paste0(site,"_",release,"_",behavior,"_settled.Rdata")))
load(here("data", paste0(site,"_",release,"_",behavior,"_settled_mids.Rdata")))
# add release hour and behavior labels
n <- nrow(get(paste0(name_long)))
assign(paste0(name_long), get(paste0(name_long)) %>%
add_column(release = rep(release,
times = n)) %>%
add_column(behavior = rep(behavior,
times = n)))
# some tracks did not have any settlement - this if statement accounts for that
if(!is.null(get(paste0(site,"_",release,"_",behavior,"_set")))){
# add release hour and behavior labels
n <- nrow(get(paste0(site,"_",release,"_",behavior,"_set")))
assign(paste0(site,"_",release,"_",behavior,"_set"),
get(paste0(site,"_",release,"_",behavior,"_set")) %>%
add_column(release = rep(release,
times = n)) %>%
add_column(behavior = rep(behavior,
times = n)))
# add release hour and behavior labels
n <- nrow(get(paste0(site,"_",release,"_",behavior,"_set_mids")))
assign(paste0(site,"_",release,"_",behavior,"_set_mids"),
get(paste0(site,"_",release,"_",behavior,"_set_mids")) %>%
add_column(release = rep(release,
times = n)) %>%
add_column(behavior = rep(behavior,
times = n)))
# find only tracks that settled
loop_set_ind <- unique(get(paste0(site,"_",release,"_",behavior,"_set"))$track)
loop_set_tracks <- get(paste0(name_long)) %>%
filter(track %in% loop_set_ind)
# find number of larvae that settled
loop_n_set <- length(loop_set_ind)
# combine into looped dataframes
map_set <- rbind(map_set, get(paste0(site,"_",release,"_",behavior,"_set")))
map_set_mids <- rbind(map_set,
get(paste0(site,"_",release,"_",behavior,"_set_mids")))
map_set_tracks <- rbind(map_set_tracks, loop_set_tracks)
# cleanup environment (these are big files)
rm(list = paste0(site,"_",release,"_",behavior,"_set"))
rm(list = paste0(site,"_",release,"_",behavior,"_set_mids"))
rm(loop_set_ind, loop_set_tracks)
} else {
loop_n_set <- 0
}
# get max size and other n labels for map
loop_n_size <- round(max(get(paste0(name_long))$size), digits = 2)
map_ns <- rbind(map_ns, c(loop_n_set,
paste0("n settled = ", loop_n_set),
loop_n_size,
paste0("max size = ", loop_n_size),
release,
behavior))
map_df <- rbind(map_df, get(paste0(name_long)))
rm(list = paste0(name_long))
rm(loop_n_size)
# set progress bar
pb_i <- pb_i + 1
setTxtProgressBar(pb, pb_i)
}
}
# set map bounds to the extremes of all larval tracks being mapped
map_xmin <- min(map_df$x) - 2*map_scale
if(map_xmin < min(WA_nodes_coords$x)) {
map_xmin <- min(WA_nodes_coords$x)
}
map_xmax <- max(map_df$x) + 2*map_scale
if(map_xmax > max(WA_nodes_coords$x)) {
map_xmax <- max(WA_nodes_coords$x)
}
map_ymin <- min(map_df$y) - map_scale
if(map_ymin < min(WA_nodes_coords$y)) {
map_ymin <- min(WA_nodes_coords$y)
}
map_ymax <- max(map_df$y) + map_scale
if(map_ymax > max(WA_nodes_coords$y)) {
map_ymax <- max(WA_nodes_coords$y)
}
# date labels
SSM_start <- as_date("2014-06-15 00:00:00")
tz(SSM_start) <- "America/Los_Angeles"
map_ns <- tibble(as.data.frame(map_ns)) %>%
rename(n_settled = V1,
n_settled_lab = V2,
max_size = V3,
max_size_lab = V4,
release = V5,
behavior = V6) %>%
mutate(n_settled = as.numeric(n_settled),
max_size = as.numeric(max_size),
release = as.numeric(release))
# rearrange where labels appear based on release site (printing test maps required to find where is a good spot on each map)
if(site %in% c("DH", "UR", "BI")) {
map_ns <- map_ns %>%
add_column(x = map_xmax - 2*map_scale,
y = map_ymax - 0.5*map_scale)
} else if(site %in% c("FB", "KH", "PB", "SimB")) {
map_ns <- map_ns %>%
add_column(x = map_xmin + 2*map_scale,
y = map_ymax - 0.5*map_scale)
} else if(site %in% c("SamB")) {
map_ns <- map_ns %>%
add_column(x = map_xmin + 2*map_scale,
y = map_ymin + 0.5*map_scale)
} else if(site %in% c("UR")) {
map_ns <- map_ns %>%
add_column(x = (map_xmin + map_xmax)/2,
y = map_ymax - 0.5*map_scale)
}
# labels and color palette
map_ns$release_time <- SSM_start + as.period(map_ns$release, unit = "hour")
release_times <- SSM_start + as.period(release_hours, unit = "hour")
n_settled <- sum(map_ns$n_settled)
pal1 <- pnw_palette("Bay", n = n_settled)
end_hours <- c(max(filter(map_df, release == release_hours[1])$hour),
max(filter(map_df, release == release_hours[2])$hour))
end_times <- SSM_start + as.period(end_hours, unit = "hour")
lab_release <- c("477" = paste0(format(release_times[1],
"%Y %b %d %H:%M"), " to \n",
format(end_times[1],
"%Y %b %d %H:%M")),
"485" = paste0(format(release_times[2],
"%Y %b %d %H:%M"), " to \n",
format(end_times[2],
"%Y %b %d %H:%M")))
lab_behavior <-c("none" = "Neutral",
"photo" = "Phototactic",
"onto" = "Ontogenetic")
# create map, first pulling basemap
psmall <- basemap_nodes +
# all larval tracks
geom_path(data = map_df, aes(x = x, y = y, group = track),
alpha = 0.075) +
# just larval tracks that settled (colored)
geom_path(data = map_set_tracks, aes(x = x, y = y, color = track),
alpha = 0.5) +
# mid-points of each settlement region
geom_point(data = map_set_mids, aes(x = x, y = y, fill = track),
shape = 21, size = 0.5, alpha = 0.5, stroke = 0.5) +
# release position
geom_point(data = filter(map_df, hour == release &
sec == min(map_df$sec) &
track == 1),
aes(x = x, y = y),
color = "darkred", shape = 8) +
# number settled label
geom_text_repel(data = map_ns, aes(x = x, y = y, label = n_settled_lab),
size = 3) +
# max size label
geom_text_repel(data = map_ns, aes(x = x, y = y, label = max_size_lab),
vjust = 2, size = 3) +
# adjust coordinate limits
coord_sf(xlim = c(map_xmin, map_xmax),
ylim = c(map_ymin, map_ymax)) +
# titles with n, max size, and total settled
labs(title = paste0("Olympia Oyster Larval Dispersal Model \n",
site_name, " Larva Tracks | 2014 SSM"),
subtitle = paste0("n = 1000 per release | ",
"max size = ", max(map_df$size)," | ",
"total settled = ", n_settled)) +
# remove legends
guides(color = "none", fill = "none") +
# apply color palettes
scale_color_manual(values = pal1) +
scale_fill_manual(values = pal1) +
# center plot titles
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5)) +
# remove x and y axis labels
xlab(NULL) +
ylab(NULL) +
# facet by release hour and by behavior
facet_grid(rows = vars(release), cols = vars(behavior),
labeller = labeller(release = lab_release,
behavior = lab_behavior))
print(psmall)
# cleanup environment
rm(psmall, map_scale, n_settled,
map_xmin, map_xmax, map_ymax, map_ymin,
map_df, map_set, map_set_mids, map_set_tracks,
map_ns, lab_release, lab_behavior,
end_hours, end_times, release_times, SSM_start)
}## | | | 0% | |============== | 20% | |============================ | 40% | |========================================== | 60% | |======================================================== | 80% | |======================================================================| 100%
## | | | 0% | |============== | 20% | |============================ | 40% | |========================================== | 60% | |======================================================== | 80% | |======================================================================| 100%
There are still some finishing touches I want to make to this before they are fully complete (like the density of the points when there are a lot of settlers), but I’m rather proud of this overall.