heaphy.R

library(leaflet)
library(rayshader)
library(raster)
library(geoviz)

# define bounding box with longitude/latitude coordinates
bbox <- list(
  p1 = list(long = 172, lat = -41.0),
  p2 = list(long = 173.1, lat = -40.4)
)

leaflet() %>%
  addTiles() %>% 
  addRectangles(
    lng1 = bbox$p1$long, lat1 = bbox$p1$lat,
    lng2 = bbox$p2$long, lat2 = bbox$p2$lat,
    fillColor = "transparent"
  ) %>%
  fitBounds(
    lng1 = bbox$p1$long, lat1 = bbox$p1$lat,
    lng2 = bbox$p2$long, lat2 = bbox$p2$lat,
  )

# -41.2317104,171.99075,10.29
# -40.5374054,172.9658431,12.29

# tmp <- rgdal::readGDAL("nztm.tif")
elev_file <- "nztm.tif"
elev_img <- raster::raster(elev_file)
elev_img
# class      : RasterLayer 
# dimensions : 19712, 13261, 261400832  (nrow, ncol, ncell)
# resolution : 80, 80  (x, y)
# extent     : 1062532, 2123412, 4705791, 6282751  (xmin, xmax, ymin, ymax)
# crs        : +proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
# source     : /Users/baptiste/Downloads/kx-nz-80m-digital-elevation-model-GTiff/nztm.tif 
# names      : nztm 
# values     : -32768, 32767  (min, max)

lat <- c(bbox$p1$lat, bbox$p2$lat)
long <- c(bbox$p1$long, bbox$p2$long)
dem <- crop_raster_track(elev_img, lat, long, width_buffer = 2)
# plot(dem)

raster_sml <- dem
raster_sml$nrow=100
raster_sml$ncol=100
elev_sml <- resample(elev_img, raster_sml, method='bilinear')

elev_matrix <- matrix(
  elev_sml, 
  nrow = ncol(elev_sml), ncol = nrow(elev_sml)
)

elev_matrix[is.na(elev_matrix)] <- 0

# calculate rayshader layers
raymat <- ray_shade(elev_matrix, zscale = 30, lambert = TRUE)
watermap <- detect_water(elev_matrix)

overlay_image <-
  slippy_overlay(raster_sml, image_source = "mapbox", image_type = 'satellite', api_key = api_key, return_png=TRUE, max_tiles=10)


scene <- elev_matrix %>%
    sphere_shade(texture = "imhof4") %>%
    add_water(watermap, color = "imhof4") %>%
    add_shadow(raymat, max_darken = 0.5) %>%
    # add_shadow(ambmat, max_darken = 0.5) %>%
  add_overlay(overlay_image)


rayshader::plot_3d(
  scene,
  elev_matrix,
  zscale = raster_zscale(raster_sml)*0.2,
  solid = TRUE,
  shadow = TRUE,
)