Last active
August 16, 2016 05:44
-
-
Save ozjimbob/fe9ea62634e73ef6ee3053af94cec864 to your computer and use it in GitHub Desktop.
Example code to great a hillshade raster from a DEM in R using spherical environment mapping for the colour scheme.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
library(raster) | |
library(jpeg) | |
# Spherical environment mapping hillshade function | |
getv=function(i,a,s){ | |
ct = dim(i)[1:2]/2 | |
sx = values(s)/90 * ct[1] | |
sy = values(s)/90 * ct[2] | |
a = values(a) * 0.01745 | |
px = floor(ct[1] + sx * -sin(a)) | |
py = floor(ct[2] + sy * cos(a)) | |
template = brick(s,s,s) | |
values(template)=NA | |
cellr = px + py * ct[1]*2 | |
cellg = px + py * ct[1]*2 + (ct[1]*2*ct[2]*2) | |
cellb = px + py * ct[1]*2 + 2*(ct[1]*2*ct[2]*2) | |
template[[1]] = i[cellr] | |
template[[2]] = i[cellg] | |
template[[3]] = i[cellb] | |
template = template * 256 | |
template | |
} | |
# Load a digital elevation raster | |
# An example can be found here: https://www.dropbox.com/s/utup77tpyrioqio/dem.tif?dl=0 | |
dem = raster("dem.tif") | |
# Generate slope and aspect rasters | |
slope=terrain(dem,opt='slope',unit='degrees') | |
aspect=terrain(dem,opt='aspect',unit='degrees') | |
# Load an environment map image | |
# An example can be found at http://i.imgur.com/9pvbHjN.jpg | |
map=readJPEG("9pvbHjN.jpg") | |
# Do the mapping | |
out = getv(map, aspect, slope) | |
# Plot pretty mountains | |
plotRGB(out) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment