This is an R Markdown Notebook. I would to share the R code chunks of the visualization of elevation in Leaflet.

It is advisable to start cleaning the memory:

rm(list=ls())

Then we will do loading the libraries:

library(elevatr)
library(sf)
library(leaflet)
library(leaflet.extras)
library(terra)
library(leaflet)
library(rasterVis)

To access raster elevation data (e.g. a DEM) the elevatr package uses the Amazon Web Services Terrain Tiles. We will explore this functionality in this notebook.

Import department boundaries. As this notebook creates a new folder named “ElevationR”, it is simple to get a list of the available data there.

#dir.create("./ElevationR")
#list.files("./ElevationR")

Get Raster Elevation Data

I used a function get_elev_raster() in the elevatr pacakge which was written to standarize access to elevation data from web APIs. It can access the USGS Elevation Point Query Service.

Example usage of each is included below. For these examples, we can create a dataset to use.

Hokkaido

# Define bounding box for Hokkaido
hokkaido_bbox <- st_as_sf(st_as_sfc(
  st_bbox(c(xmin = 139.0, ymin = 41.0, xmax = 145.5, ymax = 45.5), crs = 4326)
))

# Download DEM (returns RasterLayer, so we convert to SpatRaster)
dem_hokkaido_raster <- get_elev_raster(
  locations = hokkaido_bbox,
  z = 8,
  clip = "bbox",
  prj = "EPSG:4326"
)
## Mosaicing & Projecting
## Clipping DEM to bbox
## Note: Elevation units are in meters.
# Convert to SpatRaster
dem_hokkaido <- terra::rast(dem_hokkaido_raster)

# Mask sea elevation (set values < 1 m to NA)
dem_hokkaido[dem_hokkaido < 1] <- NA

# Get elevation range safely
range_hokkaido <- terra::global(dem_hokkaido, range, na.rm = TRUE)[1, ]

# Convert to RasterLayer for Leaflet
dem_hokkaido_r <- terra::as.raster(dem_hokkaido)

# Define color palette
pal_hokkaido <- colorNumeric(
  palette = terrain.colors(25),
  domain = c(range_hokkaido[1], range_hokkaido[2]),
  na.color = "transparent"
)
# Write the raster to a temporary GeoTIFF
#tempfile_hokkaido <- tempfile(fileext = ".tif")
terra::writeRaster(dem_hokkaido, "./ElevationR/hokkaido.tif", overwrite = TRUE)

# Load it as a RasterLayer for Leaflet
dem_hokkaido_r <- rast("./ElevationR/hokkaido.tif")
# Visualize in Leaflet
circles <- data.frame(
  lng = c(141.35, 143.2044, 144.3814),
  lat = c(43.0667, 42.9172, 42.9849),
  popup = c("Sapporo", "Obihiro", "Kushiro")
)

leaflet() |>
  addTiles() |>
  addRasterImage(dem_hokkaido_r, colors = pal_hokkaido, opacity = 0.8) |>
  addCircleMarkers(
    data = circles,
    lng = ~lng,
    lat = ~lat,
    color = "red",
    popup = ~popup
  ) |>
  addLegend(
    pal = pal_hokkaido,
    values = range_hokkaido,
    title = "Elevation (m)"
  ) |> 
  #addSimpleGraticule()
  addMiniMap(width = 120, height = 120)

Kyushu

# Define bounding box for Hokkaido
kyushu_bbox <- st_as_sf(st_as_sfc(
  st_bbox(c(xmin = 129.5, ymin = 31.0, xmax = 132.5, ymax = 34.5), crs = 4326)
))


# Download DEM (returns RasterLayer, so we convert to SpatRaster)
dem_kyushu_raster <- get_elev_raster(
  locations = kyushu_bbox,
  z = 8,
  clip = "bbox",
  prj = "EPSG:4326"
)
## Mosaicing & Projecting
## Clipping DEM to bbox
## Note: Elevation units are in meters.
# Convert to SpatRaster
dem_kyushu <- terra::rast(dem_kyushu_raster)

# Mask sea elevation (set values < 1 m to NA)
dem_kyushu[dem_kyushu < 1] <- NA

# Get elevation range safely
range_kyushu <- terra::global(dem_kyushu, range, na.rm = TRUE)[1, ]

# Convert to RasterLayer for Leaflet
dem_kyushu_r <- terra::as.raster(dem_kyushu)

# Define color palette
pal_kyushu <- colorNumeric(
  palette = terrain.colors(25),
  domain = c(range_kyushu[1], range_kyushu[2]),
  na.color = "transparent"
)
# Write the raster to a temporary GeoTIFF
#tempfile_kyushu <- tempfile(fileext = ".tif")
terra::writeRaster(dem_kyushu, "./ElevationR/kyushu.tif", overwrite = TRUE)

# Load it as a RasterLayer for Leaflet
dem_kyushu_r <- rast("./ElevationR/kyushu.tif")

# Visualize in Leaflet
leaflet() %>%
  addTiles() %>%
  addRasterImage(dem_kyushu_r, colors = pal_kyushu, opacity = 0.8) %>%
  addLegend(
    pal = pal_kyushu,
    values = unlist(range_kyushu),
    title = "Elevation (m)"
  )