dakvid/kapiti_neon.R

## kapiti_neon.R
# A 3d animated neon contour lines plot of Kapiti Island
# For #30DayMapChallenge 2021 - Day 21 - Elevation
#  -- David Friggens, 1 December 2021

# Adapting the amazing tutorial from Tyler:
# https://www.tylermw.com/pathtracing-neon-landscapes-in-r/

library(magrittr)
library(glue)
library(stringr)
library(dplyr)
library(tidyr)
library(purrr)
library(sf)
library(elevatr)
library(isoband)
library(rayshader)
library(rayrender)
library(av)


# Elevation -------------------------------------------------------------

# Luckily the island is its own Statistical Area 2, & I have those from
# https://datafinder.stats.govt.nz
ki_outline <-
  read_sf("statsnz/statistical-area-2-2021-clipped-generalised.gpkg") %>%
  filter(SA22021_V1_00_NAME_ASCII == "Kapiti Island")
ki_outline <-
  ki_outline %>%
  st_transform(crs = 4326L) %>%
  select(geom) %>%
  st_cast("POLYGON") %>%
  slice(4) # Just take the main island and ignore the pesky tiny ones

# Get elevation data
ki_raster <-
  get_elev_raster(
    locations = ki_outline,
    z = 12,
    clip = "locations")

# Convert to matrix for plotting
ki_matrix <-
  raster_to_matrix(ki_raster)


# Contours ----------------------------------------------------------------

# Note: I had to manually examine the matrix to see that the elevation range
# was -74 to 515, and then choose desired steps
# ki_matrix %>% min(na.rm = TRUE)
# ki_matrix %>% max(na.rm = TRUE)

ki_contours <-
  isolines(x = 1:ncol(ki_matrix),
           y = 1:nrow(ki_matrix),
           z = ki_matrix,
           levels = seq(-50, 500, by = 25))

# Helpful to check the contour lines at this point:

# contours <-
#   ki_contours %>%
#   iso_to_sfg()
# sf_contours <-
#   st_sf(level = names(contours),
#         geometry = st_sfc(contours))
# library(ggplot2)
# ggplot(sf_contours) + geom_sf(aes(color = level))


# Neon testing ------------------------------------------------------------

# This is just for testing, but necessary for fine tuning the details

NUM_CONTOURS <- length(ki_contours)

# Didn't have time to fully understand these, but did
# increase HEIGHT_REDUCE to shrink the gaps between
HEIGHT_LOWER <- 80
HEIGHT_REDUCE <- 10
HEIGHT_OFFSET <- 3

# If you want it to spin around you probably want to have
# (0,0) somewhere aroud your middle, rather than in a corner
X_OFFSET <- 250
Y_OFFSET <- 250

RADIUS <- 0.6

create_segment_test <-
  function(start_x, start_y, start_z,
           end_x, end_y, end_z,
           heat_color) {
    segment(start = c(start_x, start_y, start_z),
            end = c(end_x, end_y, end_z),
            radius = RADIUS,
            material = diffuse(color = heat_color))
  }

process_each_contour_test <-
  function(the_contour, the_height, the_color) {

    # (h - L) / R - O
    transformed_height <-
      the_height %>%
      as.numeric() %>%
      subtract(HEIGHT_LOWER) %>%
      divide_by(HEIGHT_REDUCE) %>%
      subtract(HEIGHT_OFFSET)

    # Don't know about you, but I find x,y + h => x,y,z confusing...
    contour_segments <-
      the_contour %>%
      as_tibble() %>%
      group_by(id) %>%
      mutate(start_x = x - X_OFFSET,
             start_y = transformed_height,
             start_z = y - Y_OFFSET,
             end_x = lead(x) - X_OFFSET,
             end_y = transformed_height,
             end_z = lead(y) - Y_OFFSET,
             heat_color = the_color) %>%
      ungroup() %>%
      drop_na() %>%
      select(start_x, start_y, start_z,
             end_x, end_y, end_z,
             heat_color) %>%
      pmap_df(create_segment_test)

    return(contour_segments)
  }

ki_scene_test <-
  pmap_df(list(the_contour = ki_contours,
               the_height = names(ki_contours),
               the_color = heat.colors(NUM_CONTOURS)),
          process_each_contour_test)

# Generate for manual inspection

# Where are we looking from and to, using the offset/transformed coordinates?
# Remember, these are x,y,z but y is height
FROM_TEST <- c(0, 110, 1000)
AT_TEST <- c(15, -1, 15)
FOV_TEST <- 25

APERTURE_TEST <- 0
WIDTH_TEST <- 1200
HEIGHT_TEST <- 800

generate_ground(material = diffuse(color="grey20")) %>%
  add_object(ki_scene_test) %>%
  render_scene(fov = FOV_TEST,
               lookfrom = FROM_TEST, lookat = AT_TEST,
               samples = 200, aperture = APERTURE_TEST,
               width = WIDTH_TEST, height = HEIGHT_TEST)


# Neon static -------------------------------------------------------------

# Now that you've figured out how the scene is layed out we can make it look nice.
# The extra bling takes longer to render, of course.

# Reusing the variables defined in the testing section
INTENSITY <- 3

create_segment <-
  function(start_x, start_y, start_z,
           end_x, end_y, end_z,
           heat_color) {
    segment(start = c(start_x, start_y, start_z),
            end = c(end_x, end_y, end_z),
            radius = RADIUS,
            material = light(intensity = INTENSITY,
                             color = heat_color))
  }

create_sphere <-
  function(start_x, start_y, start_z,
           heat_color,
           ...) {
    sphere(x = start_x,
           y = start_y,
           z = start_z,
           radius = RADIUS,
           material = light(intensity = INTENSITY,
                            color = heat_color))
  }

process_each_contour <-
  function(the_contour, the_height, the_color) {

    # (h - L) / R - O
    transformed_height <-
      the_height %>%
      as.numeric() %>%
      subtract(HEIGHT_LOWER) %>%
      divide_by(HEIGHT_REDUCE) %>%
      subtract(HEIGHT_OFFSET)

    # Don't know about you, but I find x,y + h => x,y,z confusing...
    contour_coordinates <-
      the_contour %>%
      as_tibble() %>%
      group_by(id) %>%
      mutate(start_x = x - X_OFFSET,
             start_y = transformed_height,
             start_z = y - Y_OFFSET,
             end_x = lead(x) - X_OFFSET,
             end_y = transformed_height,
             end_z = lead(y) - Y_OFFSET,
             heat_color = the_color) %>%
      ungroup() %>%
      select(start_x, start_y, start_z,
             end_x, end_y, end_z,
             heat_color)

    # Same segments as before
    contour_segments <-
      contour_coordinates %>%
      drop_na() %>%
      pmap_df(create_segment)

    # Add spheres at corners
    contour_spheres <-
      contour_coordinates %>%
      pmap_df(create_sphere)

    return(bind_rows(contour_segments,
                     contour_spheres))
  }

ki_scene <-
  pmap_df(list(the_contour = ki_contours,
               the_height = names(ki_contours),
               the_color = heat.colors(NUM_CONTOURS)),
          process_each_contour)

# Where are we looking from and to, using the offset/transformed coordinates?
# Remember, these are x,y,z but y is height
FROM_STATIC <- c(0, 110, 1000)
AT_STATIC <- c(15, -1, 15)
FOV_STATIC <- 25

APERTURE_STATIC <- 0
BLOOM_STATIC <- 5
WIDTH_STATIC <- 1200
HEIGHT_STATIC <- 800

FILE_STATIC <- "Day_21/Day_21_kapiti_neon.png"

generate_ground(material = metal(color = "grey20",
                                 fuzz = 0.05)) %>%
  add_object(ki_scene) %>%
  render_scene(fov = FOV_STATIC, bloom = BLOOM_STATIC,
               tonemap = "reinhold",
               lookfrom = FROM_STATIC, lookat = AT_STATIC,
               samples = 200, aperture = APERTURE_STATIC,
               width = WIDTH_STATIC, height = HEIGHT_STATIC,
               filename = FILE_STATIC)


# Animation ---------------------------------------------------------------

AT_ANIM <- c(15, -1, 15)
FROM_XZ_ANIM <- 800
FROM_X_OFFSET <- 0
FROM_Z_OFFSET <- 10
FROM_Y_ANIM <- 120

animation_steps <-
  tibble(frame_num = 1:360) %>%
  mutate(frame_x = FROM_XZ_ANIM * sinpi(frame_num / 180) - FROM_X_OFFSET,
         frame_z = FROM_XZ_ANIM * cospi(frame_num / 180) - FROM_Z_OFFSET)

FOV_ANIM <- 25

APERTURE_ANIM <- 0.1
BLOOM_ANIM <- 5
WIDTH_ANIM <- 1200
HEIGHT_ANIM <- 800

## to use with glue::glue
FILE_ANIM_STEP <- "Day_21/temp/kapiti_frame_{str_pad(frame_num, 3, pad = '0')}.png"
FILE_ANIM <- "Day_21/Day_21_kapiti_neon.mp4"

generate_animation_step <-
  function(frame_num, frame_x, frame_z) {
    generate_ground(material = metal(color = "grey20",
                                     fuzz = 0.05)) %>%
      add_object(ki_scene) %>%
      render_scene(lookfrom = c(frame_x, FROM_Y_ANIM, frame_z),
                   lookat = AT_ANIM,
                   samples = 200, fov = FOV_ANIM, aperture = APERTURE_ANIM,
                   tonemap = "reinhold",
                   width = WIDTH_ANIM, height = HEIGHT_ANIM,
                   filename = glue(FILE_ANIM_STEP))
  }

# Make the stills
animation_steps %>%
  pwalk(generate_animation_step)

# Combine into a video
av_encode_video(list.files(path = "Day_21/temp", pattern = "png", full.names = TRUE),
                framerate = 30,
                output = FILE_ANIM)
	# A 3d animated neon contour lines plot of Kapiti Island
	# For #30DayMapChallenge 2021 - Day 21 - Elevation
	# -- David Friggens, 1 December 2021

	# Adapting the amazing tutorial from Tyler:
	# https://www.tylermw.com/pathtracing-neon-landscapes-in-r/

	library(magrittr)
	library(glue)
	library(stringr)
	library(dplyr)
	library(tidyr)
	library(purrr)
	library(sf)
	library(elevatr)
	library(isoband)
	library(rayshader)
	library(rayrender)
	library(av)


	# Elevation -------------------------------------------------------------

	# Luckily the island is its own Statistical Area 2, & I have those from
	# https://datafinder.stats.govt.nz
	ki_outline <-
	read_sf("statsnz/statistical-area-2-2021-clipped-generalised.gpkg") %>%
	filter(SA22021_V1_00_NAME_ASCII == "Kapiti Island")
	ki_outline <-
	ki_outline %>%
	st_transform(crs = 4326L) %>%
	select(geom) %>%
	st_cast("POLYGON") %>%
	slice(4) # Just take the main island and ignore the pesky tiny ones

	# Get elevation data
	ki_raster <-
	get_elev_raster(
	locations = ki_outline,
	z = 12,
	clip = "locations")

	# Convert to matrix for plotting
	ki_matrix <-
	raster_to_matrix(ki_raster)



	# Contours ----------------------------------------------------------------

	# Note: I had to manually examine the matrix to see that the elevation range
	# was -74 to 515, and then choose desired steps
	# ki_matrix %>% min(na.rm = TRUE)
	# ki_matrix %>% max(na.rm = TRUE)

	ki_contours <-
	isolines(x = 1:ncol(ki_matrix),
	y = 1:nrow(ki_matrix),
	z = ki_matrix,
	levels = seq(-50, 500, by = 25))

	# Helpful to check the contour lines at this point:

	# contours <-
	# ki_contours %>%
	# iso_to_sfg()
	# sf_contours <-
	# st_sf(level = names(contours),
	# geometry = st_sfc(contours))
	# library(ggplot2)
	# ggplot(sf_contours) + geom_sf(aes(color = level))



	# Neon testing ------------------------------------------------------------

	# This is just for testing, but necessary for fine tuning the details

	NUM_CONTOURS <- length(ki_contours)

	# Didn't have time to fully understand these, but did
	# increase HEIGHT_REDUCE to shrink the gaps between
	HEIGHT_LOWER <- 80
	HEIGHT_REDUCE <- 10
	HEIGHT_OFFSET <- 3

	# If you want it to spin around you probably want to have
	# (0,0) somewhere aroud your middle, rather than in a corner
	X_OFFSET <- 250
	Y_OFFSET <- 250

	RADIUS <- 0.6

	create_segment_test <-
	function(start_x, start_y, start_z,
	end_x, end_y, end_z,
	heat_color) {
	segment(start = c(start_x, start_y, start_z),
	end = c(end_x, end_y, end_z),
	radius = RADIUS,
	material = diffuse(color = heat_color))
	}

	process_each_contour_test <-
	function(the_contour, the_height, the_color) {

	# (h - L) / R - O
	transformed_height <-
	the_height %>%
	as.numeric() %>%
	subtract(HEIGHT_LOWER) %>%
	divide_by(HEIGHT_REDUCE) %>%
	subtract(HEIGHT_OFFSET)

	# Don't know about you, but I find x,y + h => x,y,z confusing...
	contour_segments <-
	the_contour %>%
	as_tibble() %>%
	group_by(id) %>%
	mutate(start_x = x - X_OFFSET,
	start_y = transformed_height,
	start_z = y - Y_OFFSET,
	end_x = lead(x) - X_OFFSET,
	end_y = transformed_height,
	end_z = lead(y) - Y_OFFSET,
	heat_color = the_color) %>%
	ungroup() %>%
	drop_na() %>%
	select(start_x, start_y, start_z,
	end_x, end_y, end_z,
	heat_color) %>%
	pmap_df(create_segment_test)

	return(contour_segments)
	}

	ki_scene_test <-
	pmap_df(list(the_contour = ki_contours,
	the_height = names(ki_contours),
	the_color = heat.colors(NUM_CONTOURS)),
	process_each_contour_test)

	# Generate for manual inspection

	# Where are we looking from and to, using the offset/transformed coordinates?
	# Remember, these are x,y,z but y is height
	FROM_TEST <- c(0, 110, 1000)
	AT_TEST <- c(15, -1, 15)
	FOV_TEST <- 25

	APERTURE_TEST <- 0
	WIDTH_TEST <- 1200
	HEIGHT_TEST <- 800

	generate_ground(material = diffuse(color="grey20")) %>%
	add_object(ki_scene_test) %>%
	render_scene(fov = FOV_TEST,
	lookfrom = FROM_TEST, lookat = AT_TEST,
	samples = 200, aperture = APERTURE_TEST,
	width = WIDTH_TEST, height = HEIGHT_TEST)




	# Neon static -------------------------------------------------------------

	# Now that you've figured out how the scene is layed out we can make it look nice.
	# The extra bling takes longer to render, of course.

	# Reusing the variables defined in the testing section
	INTENSITY <- 3

	create_segment <-
	function(start_x, start_y, start_z,
	end_x, end_y, end_z,
	heat_color) {
	segment(start = c(start_x, start_y, start_z),
	end = c(end_x, end_y, end_z),
	radius = RADIUS,
	material = light(intensity = INTENSITY,
	color = heat_color))
	}

	create_sphere <-
	function(start_x, start_y, start_z,
	heat_color,
	...) {
	sphere(x = start_x,
	y = start_y,
	z = start_z,
	radius = RADIUS,
	material = light(intensity = INTENSITY,
	color = heat_color))
	}

	process_each_contour <-
	function(the_contour, the_height, the_color) {

	# (h - L) / R - O
	transformed_height <-
	the_height %>%
	as.numeric() %>%
	subtract(HEIGHT_LOWER) %>%
	divide_by(HEIGHT_REDUCE) %>%
	subtract(HEIGHT_OFFSET)

	# Don't know about you, but I find x,y + h => x,y,z confusing...
	contour_coordinates <-
	the_contour %>%
	as_tibble() %>%
	group_by(id) %>%
	mutate(start_x = x - X_OFFSET,
	start_y = transformed_height,
	start_z = y - Y_OFFSET,
	end_x = lead(x) - X_OFFSET,
	end_y = transformed_height,
	end_z = lead(y) - Y_OFFSET,
	heat_color = the_color) %>%
	ungroup() %>%
	select(start_x, start_y, start_z,
	end_x, end_y, end_z,
	heat_color)

	# Same segments as before
	contour_segments <-
	contour_coordinates %>%
	drop_na() %>%
	pmap_df(create_segment)

	# Add spheres at corners
	contour_spheres <-
	contour_coordinates %>%
	pmap_df(create_sphere)

	return(bind_rows(contour_segments,
	contour_spheres))
	}

	ki_scene <-
	pmap_df(list(the_contour = ki_contours,
	the_height = names(ki_contours),
	the_color = heat.colors(NUM_CONTOURS)),
	process_each_contour)

	# Where are we looking from and to, using the offset/transformed coordinates?
	# Remember, these are x,y,z but y is height
	FROM_STATIC <- c(0, 110, 1000)
	AT_STATIC <- c(15, -1, 15)
	FOV_STATIC <- 25

	APERTURE_STATIC <- 0
	BLOOM_STATIC <- 5
	WIDTH_STATIC <- 1200
	HEIGHT_STATIC <- 800

	FILE_STATIC <- "Day_21/Day_21_kapiti_neon.png"

	generate_ground(material = metal(color = "grey20",
	fuzz = 0.05)) %>%
	add_object(ki_scene) %>%
	render_scene(fov = FOV_STATIC, bloom = BLOOM_STATIC,
	tonemap = "reinhold",
	lookfrom = FROM_STATIC, lookat = AT_STATIC,
	samples = 200, aperture = APERTURE_STATIC,
	width = WIDTH_STATIC, height = HEIGHT_STATIC,
	filename = FILE_STATIC)



	# Animation ---------------------------------------------------------------

	AT_ANIM <- c(15, -1, 15)
	FROM_XZ_ANIM <- 800
	FROM_X_OFFSET <- 0
	FROM_Z_OFFSET <- 10
	FROM_Y_ANIM <- 120

	animation_steps <-
	tibble(frame_num = 1:360) %>%
	mutate(frame_x = FROM_XZ_ANIM * sinpi(frame_num / 180) - FROM_X_OFFSET,
	frame_z = FROM_XZ_ANIM * cospi(frame_num / 180) - FROM_Z_OFFSET)

	FOV_ANIM <- 25

	APERTURE_ANIM <- 0.1
	BLOOM_ANIM <- 5
	WIDTH_ANIM <- 1200
	HEIGHT_ANIM <- 800

	## to use with glue::glue
	FILE_ANIM_STEP <- "Day_21/temp/kapiti_frame_{str_pad(frame_num, 3, pad = '0')}.png"
	FILE_ANIM <- "Day_21/Day_21_kapiti_neon.mp4"

	generate_animation_step <-
	function(frame_num, frame_x, frame_z) {
	generate_ground(material = metal(color = "grey20",
	fuzz = 0.05)) %>%
	add_object(ki_scene) %>%
	render_scene(lookfrom = c(frame_x, FROM_Y_ANIM, frame_z),
	lookat = AT_ANIM,
	samples = 200, fov = FOV_ANIM, aperture = APERTURE_ANIM,
	tonemap = "reinhold",
	width = WIDTH_ANIM, height = HEIGHT_ANIM,
	filename = glue(FILE_ANIM_STEP))
	}

	# Make the stills
	animation_steps %>%
	pwalk(generate_animation_step)

	# Combine into a video
	av_encode_video(list.files(path = "Day_21/temp", pattern = "png", full.names = TRUE),
	framerate = 30,
	output = FILE_ANIM)