Skip to content
Permalink
master
Switch branches/tags
Go to file
 
 
Cannot retrieve contributors at this time
library(rayshader)
library(sp)
library(raster)
library(scales)
library(dplyr)
library(magick)
library(rgdal)
library(rgl)
library(av)
setwd("~/R/3dMapTutorial")
#SRTM DEMs can be downloaded from the SRTM tile Downloader:
# https://dwtkns.com/srtm30m/
#upload SRTM DEMs
elevation1 <- raster('N37W108.hgt')
elevation2 <- raster('N37W109.hgt')
#merge them together
SW_CO_Elevation <- merge(elevation1, elevation2)
#turn merged DEM into matrix for use in rayshader, then create a hieght
#shade map from them
height_shade(raster_to_matrix(SW_CO_Elevation)) %>% plot_map()
#stack the color bands into one raster and view the result
rgb_CO <- brick('Durango_raw.tif')
plotRGB(rgb_CO, scale = 255^2)
#the image is very dark
#correct Gamma Levels by taking the square root of the image
#view the result (scale is saved from last time)
rgb_CO_corrected <- sqrt(rgb_CO)
plotRGB(rgb_CO)
#check the coordinate ref. systems of the images
crs(rgb_CO_corrected)
crs(SW_CO_Elevation)
#both WGS84 longlat, that will work, but lets put it into UTM
rgb_CO_corrected_UTM <- projectRaster(rgb_CO_corrected, crs = '+proj=utm +zone=13 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0')
SW_CO_Elevation_UTM <- projectRaster(SW_CO_Elevation, crs = '+proj=utm +zone=13 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0')
#crop the DEM to the extent of the RGB
SW_CO_Elevation_UTM_Cropped <- crop(SW_CO_Elevation_UTM, rgb_CO_corrected_UTM)
SW_CO_Elevation_UTM_Cropped <- resample(SW_CO_Elevation_UTM_Cropped, rgb_CO_corrected_UTM)
#transpose the array becuase arrays work different than rasters in R
names(rgb_CO_corrected_UTM) <- c('R','G','B')
r_CO <- raster_to_matrix(rgb_CO_corrected_UTM$R)
g_CO <- raster_to_matrix(rgb_CO_corrected_UTM$G)
b_CO <- raster_to_matrix(rgb_CO_corrected_UTM$B)
SW_CO_Ele_Matrix <- raster_to_matrix(SW_CO_Elevation_UTM_Cropped)
rgb_CO_array <- array(0,dim=c(nrow(r_CO),ncol(r_CO),3))
rgb_CO_array[,,1] <- r_CO/255 #Red Layer
rgb_CO_array[,,2] <- g_CO/255 #Green Layer
rgb_CO_array[,,3] <- b_CO/255 #Blue Layer
rgb_CO_array <- aperm(rgb_CO_array, c(2,1,3))
plot_map(rgb_CO_array)
#fix contrast
rgb_CO_array <- rescale(rgb_CO_array, to=c(0,1))
plot_map(rgb_CO_array)
#plot Durango Point
findImgLat <- function(lat, ext, mat){
scaledLat <- ext[2] - ext[1]
LatProp <- (lat-ext[1])/scaledLat
return(LatProp*dim(mat)[1])
}
findImgLon <- function(lon, ext, mat){
scaledLon <- ext[4] - ext[3]
LonProp <- (lon-ext[3])/scaledLon
return(LonProp*dim(mat)[2])
}
DurLoc <- c(244631.9, 4129102)
DurImgLat <- findImgLat(lat = DurLoc[1],
ext = extent(SW_CO_Elevation_UTM_Cropped),
mat = SW_CO_Ele_Matrix)
DurImgLon <- findImgLon(lon = DurLoc[2],
ext = extent(SW_CO_Elevation_UTM_Cropped),
mat = SW_CO_Ele_Matrix)
label <- list(text = 'Durango')
#plot it 3D
plot_3d(rgb_CO_array, SW_CO_Ele_Matrix, windowsize = c(1600,900),
zscale = 10, shadowdepth = -50, zoom = 0.4, theta = -45,
fov = 10, background = '#F2E1D0', shadowcolor = '#523E2B',
solid = FALSE)
render_label(SW_CO_Ele_Matrix,x = DurImgLat, y = DurImgLon + 1050, z =1000, zscale =10,
text = label$text, textsize = 3, textcolor = 'white', linewidth = 3,
freetype = F)
render_snapshot(title_text = 'Southwest Colorado',
title_bar_color = '#1f5214', title_color = 'white',
title_bar_alpha = 1)
render_movie(filename ='Durango.mp4', title_text = 'Southwest Colorado | Landsat 7 Natural Color 30m | SRTM DEM',
title_color = 'white', title_bar_color = '#1f5214', title_bar_alpha = 1)