ice.r

library(dplyr)
library(rayshader)
library(rayrender)
library(sf)
library(fasterize)
library(raster)
library(exactextractr)
library(rayrender)
library(terra)
library(antanym)

g <- an_read(cache = "session")
an_filter(g, query = "^Shackleton Fracture Zone")[, c("place_name", "longitude", "latitude")]
shakleton = c(-60, -60) 
an_filter(g, query = "^Amund")[, c("place_name", "longitude", "latitude")]
amundsen_sea = c(-112, -73) 
an_filter(g, query = "^Pine")[, c("place_name", "longitude", "latitude")]
pine_glacier = c(-101, -75) 
an_filter(g, query = "^Thwai")[, c("place_name", "longitude", "latitude")]
thwaites_glacier = c(-107, -75.5) 

flight = matrix(c(shakleton, amundsen_sea, pine_glacier, thwaites_glacier), ncol = 2, byrow = TRUE)
flight =  st_multipoint(flight) |> st_sfc(crs = 4326)  |>  st_set_crs(4326) |> st_transform(crs=3031)
poi <- st_coordinates(flight)

sea_bed = rast("seabed.tif")
ice_thickness = rast("thickness.tif")

plot(ice_thickness)

sea_bed <- project(sea_bed, ice_thickness)
sea_bed <- focal(sea_bed, w=911, fun=mean, na.policy="only", na.rm=T)


ice_thickness[ice_thickness < 1] <- NA
plot(ice_thickness)

plot(cover(ice_thickness, sea_bed))

ice_thickness <- cover(ice_thickness, sea_bed)

plot(ice_thickness)
plot(st_geometry(flight), add = TRUE)
flight = st_read('flight.gpkg')

raster_matrix <- raster::raster(ice_thickness) %>%
  rayshader::raster_to_matrix(verbose = FALSE)%>%
  rayshader::resize_matrix(0.1)

# set missing values to -1 (assume this is water)
# raster_matrix[raster_matrix < 1] <- NA

# this can be scaled to exaggerate elevation
scaled_zscale <- 2 * 1

base_map2 <- raster_matrix %>%
  height_shade(texture = scico::scico(30, palette = 'devon')) 

rayshader::plot_map(base_map2)

try(rgl::rgl.close())

#2.5 for higher res
base_map2 %>%
  rayshader::plot_3d(
    raster_matrix,
    windowsize = c(800, 800),
    zscale = scaled_zscale + 50,
    baseshape = "circle",
    solid = FALSE,
    soliddepth = 0,
    shadowdepth = -2,
    #water = TRUE,
    #wateralpha = 0.1,
    phi = 90, 
    zoom = 1,  
    theta = 0, 
    shadow = TRUE
  )

render_camera(zoom = .6)
rayrender::generate_camera_motion()

render_highquality(
  sample_method = "sobol_blue",
  samples = 16, 
  light = TRUE,
  lightdirection = rev(c(220, 220, 230, 230)),
  lightcolor = c("#c2d8eb", "white", "#c2ebea", "white"),
  lightintensity = c(750, 100, 1000, 100),
  lightaltitude = c(10, 80, 10, 80),
  width = 800, height = 800
)

own_flight = get_saved_keyframes()

camera_motion =  generate_camera_motion(positions = own_flight[,1:3], 
                                        offset_lookat = 1, 
                                        fovs = 80,
                                        type = "bezier",
                                        curvature_adjust = 'both',
                                        damp_motion = TRUE,
                                        damp_magnitude = 0.7,
                                        frames = 100) 
    
# own_flight_sf = st_cast(own_flight_sf, "LINESTRING")
# own_flight_sf = st_cast(own_flight_sf, "MULTILINESTRING")

plot(own_flight_sf)
plot(flight, add = TRUE)

render_camera(phi = 20, zoom = .6, theta = 0)
render_label(raster_matrix, lat = poi[2,2], long = poi[2,1], 
             extent = ext(ice_thickness),
             zscale = scaled_zscale + 160, text = "Amundsen Sea", textcolor = "black", linecolor="grey75",
             dashed = TRUE, clear_previous = TRUE)
render_label(raster_matrix, lat = poi[3,2], long = poi[3,1],
             extent = ext(ice_thickness),
             zscale = scaled_zscale + 120, text = "Pine Glacier", textcolor = "black", linecolor="grey75",
             dashed = TRUE)
render_label(raster_matrix, lat = poi[4,2], long = poi[4,1],
             extent = ext(ice_thickness),
             zscale = scaled_zscale + 80, text = "Thwaites Glacier", textcolor = "black", linecolor="grey75",
             dashed = TRUE)

# render_path(extent = ext(ice_thickness),
#             heightmap = raster_matrix,
#             lat = unlist(st_coordinates(flight))[,2],
#             long = unlist(st_coordinates(flight))[,1],
#             zscale=50,
#             color="red",
#             antialias=TRUE,
#             offset=100, linewidth=2)

camera_path = convert_path_to_animation_coords(extent = ext(ice_thickness),
                                               heightmap = raster_matrix,
                                               #altitude = 4000,
                                               lat = unlist(st_coordinates(flight))[,2],
                                               long = unlist(st_coordinates(flight))[,1],
                                               fovs = 80,
                                               zscale=50, 
                                               type = "bezier",
                                               curvature_adjust = 'both',
                                               damp_motion = TRUE,
                                               follow_camera = TRUE,
                                               damp_magnitude = 0.7,
                                               frames = 360)

camera_path = rbind(camera_motion, camera_path)

render_highquality(
  # We test-wrote to this file above, so we know it's good
  sprintf("%s/frame", 'plots'), 
  animation_camera_coords = camera_path,
  # See rayrender::render_scene for more info, but best
  # sample method ('sobol') works best with values over 256
  sample_method = "sobol",
  use_extruded_paths = TRUE,
  samples = 450, 
  preview = FALSE,
  light = TRUE,
  lightdirection = rev(c(220, 220, 230, 230)),
  lightcolor = c("#c2d8eb", "white", "#c2ebea", "white"),
  lightintensity = c(750, 100, 1000, 100),
  lightaltitude = c(10, 80, 10, 80),
  # 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 = 1.5,
  # # Rotate the light -- positive values move it counter-clockwise
  # rotate_env = 130,
  # 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),
  width = 2000, height = 2000
)

system("ffmpeg -framerate 30 -i plots/frame%d.png -pix_fmt yuv420p ice.mp4")

library(gifski)
gifski::gifski()
png_files <- list.files("plots/frame", pattern = ".*png$", full.names = TRUE)


images <- image_join(png_files)
animation <- image_animate(images, fps = 20)
image_write(animation, "test.gif")