I'm trying to improve the look of Rayshader by overlaying more recent (higher detail) satellite imagery (that I'm getting from the {leaflet} packages) but the overlay doesn't match with the 3D rendering.
Ideally I'm looking for a open-source solution that can get global satellite imagery. Bonus points if you find finer detail data for my area of interest - Hawaii.
One method using {geoviz} and {rayshader} uses the slippy_overlay() function to create a number of overlay images from either Mapbox (satellite, mapbox-streets-v8, mapbox-terrain-v2, mapbox-traffic-v1, terrain-rgb, mapbox-incidents-v1) or Stamen. Although I found mapbox-terrain-v2 the best it still lacks the detail I would like. Since it requires setting up an API for mapbox I just use stamen/watercolor below:
library(geoviz)
library(rayshader)
### Maui
lat = 20.785700
lon = -156.259204
square_km = 22
max_tiles = 10
dem <- mapzen_dem(lat, lon, square_km, max_tiles)
elev_matrix = matrix(
raster::extract(dem, raster::extent(dem), buffer=1000),
nrow = ncol(dem),
ncol = nrow(dem)
)
ambmat <- ambient_shade(elev_matrix, zscale = 30)
raymat <- ray_shade(elev_matrix, zscale = 30, lambert = TRUE)
watermap <- detect_water(elev_matrix)
overlay_img <-
slippy_overlay(dem,
image_source = "stamen",
image_type = "watercolor",
png_opacity = 0.3,
max_tiles = max_tiles)
elev_matrix %>%
sphere_shade(sunangle = 270, texture = "imhof4") %>%
add_water(detect_water(elev_matrix), color="imhof4") %>%
add_shadow(ray_shade(elev_matrix,zscale=3,maxsearch = 300),0.5) %>%
add_shadow(ambmat,0.5) %>%
add_overlay(overlay_img) %>%
plot_3d(elev_matrix,
solid = T,
water = T,
waterdepth = 0,
wateralpha = 0.5,
watercolor = "lightblue",
waterlinecolor = "white",
waterlinealpha = 0.5,
zscale= raster_zscale(dem) / 3,
fov=0,theta=135,zoom=0.75,phi=45, windowsize = c(1000,800))
I'm trying to adapt Will Bishop's workflow for getting overlays with the leaflet package but the result is very odd. Will's approach is a bit different as it fetches elevation data from USGS, which doesn't have baythmetric elevation which is must for me - so I used geoviz
library(leaflet)
# define bounding box with longitude/latitude coordinates
bbox <- list(
p1 = list(long = -156.8037, lat = 20.29737),
p2 = list(long = -155.7351, lat = 21.29577)
)
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,
)
What's the area of my hillshade from geoviz?
dim(dem)
780 780 1
Okay so the overlay image needs to be 780 x 780 so I modify the helper functions to download the overlay with the World_Imagery base map:
define_image_size <- function(bbox, major_dim = 780) {
# calculate aspect ration (width/height) from lat/long bounding box
aspect_ratio <- abs((bbox$p1$long - bbox$p2$long) / (bbox$p1$lat - bbox$p2$lat))
# define dimensions
img_width <- ifelse(aspect_ratio > 1, major_dim, major_dim*aspect_ratio) %>% round()
img_height <- ifelse(aspect_ratio < 1, major_dim, major_dim/aspect_ratio) %>% round()
size_str <- paste(img_width, img_height, sep = ",")
list(height = img_height, width = img_width, size = size_str)
}
get_arcgis_map_image <- function(bbox, map_type = "World_Imagery", file = NULL,
width = 780, height = 780, sr_bbox = 4326) {
require(httr)
require(glue)
require(jsonlite)
url <- parse_url("https://utility.arcgisonline.com/arcgis/rest/services/Utilities/PrintingTools/GPServer/Export%20Web%20Map%20Task/execute")
# define JSON query parameter
web_map_param <- list(
baseMap = list(
baseMapLayers = list(
list(url = jsonlite::unbox(glue("https://services.arcgisonline.com/ArcGIS/rest/services/{map_type}/MapServer",
map_type = map_type)))
)
),
exportOptions = list(
outputSize = c(width, height)
),
mapOptions = list(
extent = list(
spatialReference = list(wkid = jsonlite::unbox(sr_bbox)),
xmax = jsonlite::unbox(max(bbox$p1$long, bbox$p2$long)),
xmin = jsonlite::unbox(min(bbox$p1$long, bbox$p2$long)),
ymax = jsonlite::unbox(max(bbox$p1$lat, bbox$p2$lat)),
ymin = jsonlite::unbox(min(bbox$p1$lat, bbox$p2$lat))
)
)
)
res <- GET(
url,
query = list(
f = "json",
Format = "PNG32",
Layout_Template = "MAP_ONLY",
Web_Map_as_JSON = jsonlite::toJSON(web_map_param))
)
if (status_code(res) == 200) {
body <- content(res, type = "application/json")
message(jsonlite::toJSON(body, auto_unbox = TRUE, pretty = TRUE))
if (is.null(file))
file <- tempfile("overlay_img", fileext = ".png")
img_res <- GET(body$results[[1]]$value$url)
img_bin <- content(img_res, "raw")
writeBin(img_bin, file)
message(paste("image saved to file:", file))
} else {
message(res)
}
invisible(file)
}
Now download the file, then load it
image_size <- define_image_size(bbox, major_dim = 780)
# fetch overlay image
overlay_file <- "maui_overlay.png"
get_arcgis_map_image(bbox, map_type = "World_Imagery", file = overlay_file,
# width = image_size$width, height = image_size$height,
sr_bbox = 4326)
overlay_img <- png::readPNG("maui_overlay.png")
Okay let's make the plot
elev_matrix %>%
sphere_shade(sunangle = 270, texture = "imhof4") %>%
add_water(detect_water(elev_matrix), color="imhof4") %>%
add_shadow(ray_shade(elev_matrix,zscale=3,maxsearch = 300),0.5) %>%
add_shadow(ambmat,0.5) %>%
add_overlay(overlay_img, alphacolor = 1) %>%
plot_3d(elev_matrix,
solid = T,
water = T,
waterdepth = 0,
wateralpha = 0.5,
watercolor = "lightblue",
waterlinecolor = "white",
waterlinealpha = 0.5,
zscale= raster_zscale(dem) / 3,
fov=0,theta=135,zoom=0.75,phi=45, windowsize = c(1000,800))
As you can see the overlay image is rotated to the hillshade.
Now I'm also realizing that fetching satellite with a bounding box method isn't ideal when you're trying to show bathymatrix data. It would be ideal to subset this overlay somehow programmatically but I'll probably just end up using inkscape once I've figured out how to rotate the overlay.
I tried to use the {magick}'s image_rotate() function to no avail:
library(magick)
maui <- magick::image_read("maui_overlay.png")
image_rotate(maui, 30) # -> maui_30
# image_write(maui_30, path = "maui_overlay_30.png", format = "png")
But magick has changed the dimensions:
# A tibble: 1 x 7
format width height colorspace matte filesize density
<chr> <int> <int> <chr> <lgl> <int> <chr>
1 PNG 1068 1068 sRGB TRUE 0 38x38
And will give an error with rayshader:
overlay_img <- png::readPNG("maui_overlay_30.png")
elev_matrix %>%
sphere_shade(sunangle = 270, texture = "imhof4") %>%
add_water(detect_water(elev_matrix), color="imhof4") %>%
add_shadow(ray_shade(elev_matrix,zscale=3,maxsearch = 300),0.5) %>%
add_shadow(ambmat,0.5) %>%
add_overlay(overlay_img, alphacolor = 1) %>%
plot_3d(elev_matrix,
solid = T,
water = T,
waterdepth = 0,
wateralpha = 0.5,
watercolor = "lightblue",
waterlinecolor = "white",
waterlinealpha = 0.5,
zscale= raster_zscale(dem) / 3,
fov=0,theta=135,zoom=0.75,phi=45, windowsize = c(1000,800))
Error in add_overlay(., overlay_img, alpha = 0.8) : argument 3 matches multiple formal arguments
