Skip to content

Instantly share code, notes, and snippets.

@h-a-graham
Created May 5, 2021 08:59
Show Gist options
  • Star 3 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save h-a-graham/7c4e968b7fc94908d8d3725bafad63f4 to your computer and use it in GitHub Desktop.
Save h-a-graham/7c4e968b7fc94908d8d3725bafad63f4 to your computer and use it in GitHub Desktop.
library(rayshader)
library(sf)
library(ggmap)
library(EAlidaR)
# ------------ User Inputs ---------------------------------
# Googple maps API KEY
register_google(key = "[PUT YOUR API KEY HERE!!]") # Create your own key as shown in the ggmap documentation: https://github.com/dkahle/ggmap
# WGS84 lat Long centroid location for your scene...
Loc_Coords <- c(54.45467816364232, -3.2118148812862377) #Scafell Pike
# ------ Get Aerial image and extent ---------------
X <- Loc_Coords[2]
Y <- Loc_Coords[1]
amap <- get_googlemap(center = c(lon = X , lat = Y), zoom = 15, scale = 2,
maptype ='satellite', color = 'color', size =c(640, 640))
#generate sf object of ggmap extent
bb <- attr(amap, "bb")
bbox <- as.numeric(unlist(bb2bbox(bb)))
boxcoord <- matrix(c(bbox[1],bbox[2], bbox[1],bbox[4], bbox[3],bbox[4],
bbox[3],bbox[2], bbox[1],bbox[2]), ncol=2, byrow = T)
poly <- st_polygon(list(boxcoord)) %>%
st_sfc()%>%
st_set_crs(value=4326) %>%
st_transform(crs=27700) %>%
st_as_sf()
# get basemap as raster
basemap <- ggmap(amap, extent = "device")
# ggmap(amap)
# save sat png
fname <- tempfile(fileext = '.png')
grDevices::png(filename = fname, width = 1280, height = 1280)
par(mar = c(0,0,0,0))
basemap
dev.off()
overlay_img <- png::readPNG(fname)
overlay_img_contrast = scales::rescale(overlay_img,to=c(0,1))
# ---- download Lidar ----
DTM <- get_area(poly, resolution = 1, model_type = 'DSM',
merge_tiles=TRUE, crop=TRUE)
DSM <- get_area(poly, resolution = 1, model_type = 'DTM',
merge_tiles=TRUE, crop=TRUE)
merge_DTM_DSM <- function(DEM1, DEM2){
raster::overlay(DEM1, DEM2, fun = function(x, y) {
x[is.na(x)] <- y[is.na(x)]
return(x)})
}
DEMmerge <- merge_DTM_DSM(DEM1 = DTM, DEM2 = DSM)
# ----- Rayshade --------
elmat = raster_to_matrix(DEMmerge)
plot_3d(overlay_img_contrast, elmat, zscale=1, theta=-45, phi=25, zoom=0.6, windowsize =2400,
solid=F)
Sys.sleep(2)
render_highquality(filename = 'ScafellPike_rayshadeGGMAP2.png',
lightdirection = c(60,120), lightaltitude=c(90,25),
lightintensity=c(300, 500), lightcolor = c("white", "#FF8D4B"))
render_movie(filename = 'Scafell_mov')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment