black canyon and great sand dunes

R/portraits/black_canyon/inset.R

@@ -0,0 +1,52 @@
+# inset map -- need to comment this code...
+header <- readRDS("R/portraits/black_canyon/header.rds")
+
+world <- rnaturalearth::ne_countries(scale = "small", returnclass = "sf")
+
+coords <- header$coords
+
+prj <- glue("+proj=ortho +lat_0={coords[2]} +lon_0={coords[1]} +x_0=0 +y_0=0 +a=6375000 +b=6375000 +units=m +no_defs")
+
+water <- st_sfc(st_point(c(0, 0)), crs = prj) %>%
+  st_buffer(., 6371000) %>%
+  st_transform(crs = 4326)
+
+circle_coords <- st_coordinates(water)[, c(1,2)]
+circle_coords <- circle_coords[order(circle_coords[, 1]),]
+circle_coords <- circle_coords[!duplicated(circle_coords),]
+
+rectangle <- list(rbind(circle_coords,
+                        c(X = 180, circle_coords[nrow(circle_coords), 'Y']),
+                        c(X = 180, Y = 90),
+                        c(X = -180, Y = 90),
+                        c(X = -180, circle_coords[1, 'Y']),
+                        circle_coords[1, c('X','Y')])) %>% 
+  st_polygon() %>% st_sfc(crs = 4326)
+
+
+rectangle %>%
+  ggplot()+
+  geom_sf(data = world) +
+  geom_sf(color = "red", fill = alpha("red", .5)) 
+
+sf::sf_use_s2(FALSE)
+w <- st_intersection(world, rectangle)
+w %>%
+  ggplot() +
+  geom_sf()
+
+spot <- tibble(x = coords[1], y = coords[2]) |> 
+  st_as_sf(coords = c("x", "y"), crs = 4326)
+
+colors <- header$colors
+water_color <- colors[7]
+text_color <- colors[1]
+
+loc_plot <- ggplot(data = w) +
+  geom_sf(data = water, color = NA, fill = alpha("white", .75)) +
+  geom_sf(fill = "red", size = .1, color = "white") +
+  coord_sf(crs = prj) +
+  theme_void() 
+
+loc_plot
+ggsave(loc_plot, filename = glue("images/{header$map}/{header$pal}_inset.png"), w = 4*1.5, h = 3*1.5)

R/portraits/black_canyon/markup.R

@@ -0,0 +1,39 @@
+# Load packages and `add_main_annotations` function
+for (f in list.files("R/utils")) {
+  source(paste0("R/utils/", f, collapse = ""))
+}
+
+# Load `header` list with needed data
+header <- readRDS("R/portraits/black_canyon/header.rds")
+
+# Take original graphic from `render_highquality` and
+# add annotations.
+
+add_main_annotations(map = header$map, pal = header$pal, 
+                     text_color = header$colors[1], 
+                     align = "center",
+                     base_coords = c(.725, .1),
+                     offset = 650,
+                     main_text = "Capitol Reef",
+                     main_size = 400,
+                     main_font = "Tapestry",
+                     secondary_text = "National Park",
+                     secondary_size = 200,
+                     secondary_font = "Tapestry",
+                     caption_size = 80,
+                     caption_font = "Tapestry",
+                     caption_coords = c(.5, .97),
+                     caption_align = "center",
+                     twitter_icon_coords = c(-500, 100),
+                     twitter_icon_size = 60,
+                     data_source = "USGS and AWS Terrain Tiles",
+                     original = header$outfile,
+                     crop = NULL,
+                     #crop_gravity = "north",
+                     crop_start = NULL,
+                     svg_file = NULL,
+                     svg_coords = NULL,
+                     svg_size = NULL,
+                     inset = glue("images/{header$map}/{header$pal}_inset.png"),
+                     inset_coords = c(.725, .3),
+                     inset_size = 1000)

R/portraits/black_canyon/render_graphic.R

@@ -0,0 +1,209 @@
+library(sf)
+library(tidyverse)
+library(elevatr)
+library(rayshader)
+library(glue)
+library(colorspace)
+library(NatParksPalettes)
+library(MetBrewer)
+library(scico)
+
+###################################
+# Set up polygon for clipping DEM #
+###################################
+
+# Set map name that will be used in file names, and 
+# to get get boundaries from master NPS list
+
+map <- "black_canyon"
+
+data <- st_read("data/nps_boundary/nps_boundary.shp") |>
+  filter(str_detect(PARKNAME, str_to_title(str_replace(map, "_", " "))))|> 
+  st_transform(crs = 5070)
+
+d_cent <- st_centroid(data) |> 
+  st_transform(crs = 4326)
+c <- d_cent |> 
+  mutate(long = unlist(map(d_cent$geometry,1)),
+         lat = unlist(map(d_cent$geometry,2)))  |> 
+  as_tibble()
+coords <- c(c[[1,"long"]], c[[1, "lat"]])
+
+# Example for using specific coordinates to create a geometry, here showing
+# coords for the Eye of the Sahara
+
+# data <- st_sfc(st_point(c(-11.401515762405115, 21.127630356438505)), crs = 4326) %>%
+#   st_transform(crs = 3081) %>%
+#   st_buffer(50000)
+
+# Plot to review
+data |>
+  ggplot() +
+  geom_sf()
+
+################
+# Download DEM #
+################
+
+# Get DEM using `elevatr` package. Larger Z value (max is 14)
+# results in greater resolution. Higher resolution takes more compute, though -- 
+# I can't always max `z` up to 14 on my machine. 
+
+z <- 14
+zelev <- get_elev_raster(data, z = z, clip = "location")
+mat <- raster_to_matrix(zelev)
+
+# When initially building your object to render, you'll want to work with
+# slimmed down data so you can iterate faster. I prefer to just start with
+# a `z` value of 10 above, but an alternative is to create a smaller matrix
+# with rayshader::resize_matrix().
+
+# small <- resize_matrix(mat, .25)
+
+# Set up color palette. The `pal` argument will be used in file names,
+# so it's important. `colors` will also be passed along. 
+
+
+pal <- "acadia"
+
+c1 <- natparks.pals("Acadia")
+
+colors <- c(c1)
+swatchplot(colors)
+
+
+
+# Calculate the aspect ratio of the plot so you can translate the dimensions
+
+w <- nrow(mat)
+h <- ncol(mat)
+
+# Scale so longer side is 1
+
+wr <- w / max(c(w,h))
+hr <- h / max(c(w,h))
+
+###################
+# Build 3D Object #
+###################
+
+# setting shadow to 500 feet below minimum value in DEM
+shadow_depth <- min(mat, na.rm = TRUE)
+
+# setting resolution to about 5x for height
+res <- mean(round(terra::res(zelev))) / 5
+
+# Keep this line so as you're iterating you don't forget to close the
+# previous window
+
+try(rgl::rgl.close())
+
+# Create the initial 3D object
+
+mat %>%
+  # This adds the coloring, we're passing in our `colors` object
+  height_shade(
+    texture = grDevices::colorRampPalette(colors)(256)
+  ) %>%
+  plot_3d(heightmap = mat,
+          # This is my preference, I don't love the `solid` in most cases
+          solid = FALSE,
+          # You might need to hone this in depending on the data resolution;
+          # lower values exaggerate the height
+          z = res,
+          # Set the location of the shadow, i.e. where the floor is.
+          # This is on the same scale as your data, so call `zelev` to see the
+          # min/max, and set it however far below min as you like.
+          shadowdepth = shadow_depth,
+          # Set the window size relatively small with the dimensions of our data.
+          # Don't make this too big because it will just take longer to build,
+          # and we're going to resize with `render_highquality()` below.
+          windowsize = c(800,800), 
+          # This is the azimuth, like the angle of the sun.
+          # 90 degrees is directly above, 0 degrees is a profile view.
+          phi = 90, 
+          zoom = 1, 
+          # `theta` is the rotations of the map. Keeping it at 0 will preserve
+          # the standard (i.e. north is up) orientation of a plot
+          theta = 0, 
+          background = "white") 
+
+# Use this to adjust the view after building the window object
+render_camera(phi = 80, zoom = 1, theta = 0)
+
+###############################
+# Create High Quality Graphic #
+###############################
+
+# You should only move on if you have the object set up
+# as you want it, including colors, resolution, viewing position, etc.
+
+# Ensure dir exists for these graphics
+if (!dir.exists(glue("images/{map}"))) {
+  dir.create(glue("images/{map}"))
+}
+
+# Set up outfile where graphic will be saved.
+# Note that I am not tracking the `images` directory, and this
+# is because these files are big enough to make tracking them on
+# GitHub difficult. 
+outfile <- str_to_lower(glue("images/{map}/{map}_{pal}_z{z}.png"))
+
+# Now that everything is assigned, save these objects so we
+# can use then in our markup script
+saveRDS(list(
+  map = map,
+  pal = pal,
+  z = z,
+  colors = colors,
+  outfile = outfile,
+  coords = coords
+), glue("R/portraits/{map}/header.rds"))
+
+{
+  # Test write a PNG to ensure the file path is good.
+  # You don't want `render_highquality()` to fail after it's 
+  # taken hours to render.
+  if (!file.exists(outfile)) {
+    png::writePNG(matrix(1), outfile)
+  }
+  # I like to track when I start the render
+  start_time <- Sys.time()
+  cat(glue("Start Time: {start_time}"), "\n")
+  render_highquality(
+    # We test-wrote to this file above, so we know it's good
+    outfile, 
+    # See rayrender::render_scene for more info, but best
+    # sample method ('sobol') works best with values over 256
+    samples = 450, 
+    # Turn light off because we're using environment_light
+    light = TRUE,
+    # lightdirection = rev(c(95, 95, 85, 85)),
+    lightdirection = rev(c(145, 145, 205, 205)),
+    lightcolor = c(colors[6], "white", colors[6], "white"),
+    lightintensity = c(750, 75, 750, 75),
+    lightaltitude = c(15, 80, 15, 80),
+    # All it takes is accidentally interacting with a render that takes
+    # hours in total to decide you NEVER want it interactive
+    interactive = FALSE,
+    preview = FALSE,
+    # HDR lighting used to light the scene
+    environment_light = "assets/env/phalzer_forest_01_4k.hdr",
+    # environment_light = "assets/env/small_rural_road_4k.hdr",
+    # Adjust this value to brighten or darken lighting
+    intensity_env = .5,
+    # Rotate the light -- positive values move it counter-clockwise
+    rotate_env = -20,
+    # This effectively sets the resolution of the final graphic,
+    # because you increase the number of pixels here.
+    # width = round(6000 * wr), height = round(6000 * hr),
+    ground_material = rayrender::diffuse(color = colors[1], noisecolor = colors[3]),
+    width = 8000, height = 8000
+  )
+  end_time <- Sys.time()
+  cat(glue("Total time: {end_time - start_time}"))
+}
+
+
+
+

R/portraits/great_sand_dunes/inset.R

@@ -0,0 +1,52 @@
+# inset map -- need to comment this code...
+header <- readRDS("R/portraits/great_sand_dunes/header.rds")
+
+world <- rnaturalearth::ne_countries(scale = "small", returnclass = "sf")
+
+coords <- header$coords
+
+prj <- glue("+proj=ortho +lat_0={coords[2]} +lon_0={coords[1]} +x_0=0 +y_0=0 +a=6375000 +b=6375000 +units=m +no_defs")
+
+water <- st_sfc(st_point(c(0, 0)), crs = prj) %>%
+  st_buffer(., 6371000) %>%
+  st_transform(crs = 4326)
+
+circle_coords <- st_coordinates(water)[, c(1,2)]
+circle_coords <- circle_coords[order(circle_coords[, 1]),]
+circle_coords <- circle_coords[!duplicated(circle_coords),]
+
+rectangle <- list(rbind(circle_coords,
+                        c(X = 180, circle_coords[nrow(circle_coords), 'Y']),
+                        c(X = 180, Y = 90),
+                        c(X = -180, Y = 90),
+                        c(X = -180, circle_coords[1, 'Y']),
+                        circle_coords[1, c('X','Y')])) %>% 
+  st_polygon() %>% st_sfc(crs = 4326)
+
+
+rectangle %>%
+  ggplot()+
+  geom_sf(data = world) +
+  geom_sf(color = "red", fill = alpha("red", .5)) 
+
+sf::sf_use_s2(FALSE)
+w <- st_intersection(world, rectangle)
+w %>%
+  ggplot() +
+  geom_sf()
+
+spot <- tibble(x = coords[1], y = coords[2]) |> 
+  st_as_sf(coords = c("x", "y"), crs = 4326)
+
+colors <- header$colors
+water_color <- colors[7]
+text_color <- colors[1]
+
+loc_plot <- ggplot(data = w) +
+  geom_sf(data = water, color = NA, fill = alpha("white", .75)) +
+  geom_sf(fill = "red", size = .1, color = "white") +
+  coord_sf(crs = prj) +
+  theme_void() 
+
+loc_plot
+ggsave(loc_plot, filename = glue("images/{header$map}/{header$pal}_inset.png"), w = 4*1.5, h = 3*1.5)

R/portraits/great_sand_dunes/markup.R

@@ -0,0 +1,39 @@
+# Load packages and `add_main_annotations` function
+for (f in list.files("R/utils")) {
+  source(paste0("R/utils/", f, collapse = ""))
+}
+
+# Load `header` list with needed data
+header <- readRDS("R/portraits/great_sand_dunes/header.rds")
+
+# Take original graphic from `render_highquality` and
+# add annotations.
+
+add_main_annotations(map = header$map, pal = header$pal, 
+                     text_color = header$colors[1], 
+                     align = "center",
+                     base_coords = c(.725, .1),
+                     offset = 650,
+                     main_text = "Capitol Reef",
+                     main_size = 400,
+                     main_font = "Tapestry",
+                     secondary_text = "National Park",
+                     secondary_size = 200,
+                     secondary_font = "Tapestry",
+                     caption_size = 80,
+                     caption_font = "Tapestry",
+                     caption_coords = c(.5, .97),
+                     caption_align = "center",
+                     twitter_icon_coords = c(-500, 100),
+                     twitter_icon_size = 60,
+                     data_source = "USGS and AWS Terrain Tiles",
+                     original = header$outfile,
+                     crop = NULL,
+                     #crop_gravity = "north",
+                     crop_start = NULL,
+                     svg_file = NULL,
+                     svg_coords = NULL,
+                     svg_size = NULL,
+                     inset = glue("images/{header$map}/{header$pal}_inset.png"),
+                     inset_coords = c(.725, .3),
+                     inset_size = 1000)