Light, Elevation and Mountains

Posted by Mike Birdgeneau

Category: Data

This is the start of a project I’ve been wanting to put together for quite a while: Using public digital elevation data in the mountains to figure out when I should head out to take photos of different mountain faces. (I also plan to extend this further to look for good slopes for skiing, with consideration to safety and overhead hazard!)

Luckily, R has functionality to handle digital elevation maps, as well as functionality to calculate aspect, slope, and hill-shading (which is exactly what we are looking for).

Digital Elevation Maps

In order to plot the elevation data, I needed to find a good source of data for the areas I’m interested in. Fortunately, Natural Resources Canada (NRCAN) has a free service called GeoBase which contains digital elevation maps (DEM) that provide exactly this source of information.

library(raster)
library(RColorBrewer)
library(ggplot2)
library(manipulate) # Included with RStudio
library(lubridate)

if(!file.exists("082o03_0201_demw.dem")){
  download.file("ftp://ftp2.cits.rncan.gc.ca/pub/geobase/official/cded/50k_dem/082/082o03.zip","082o03.zip")
  unzip("082o03.zip")
}

You can find this data here: GeoBase Digital Elevation Data

Plotting Elevation

Once we have this data, we can load the data into R by using the raster package. This is pretty simple, and we can plot it using the image function within the raster package, but it’s also pretty nice to use Hadley’s ggplot2 package to plot the elevation data.

dem<-mosaic(raster("082o03_0201_demw.dem"),raster("082o03_0201_deme.dem"),fun=mean)
DEM<-data.frame(rasterToPoints(dem))
colnames(DEM) <- c("X","Y","Elevation")

image(dem)

Using image, we get this plot showing us that this all worked properly:

Digital Elevation Model

This is much prettier in ggplot with some terrain colours, but you can see it is definitely missing the lighting to give the mountains some depth. We’ll look at this next.

ggplot(DEM,aes(x=X,y=Y))+
    geom_raster(aes(fill=layer),alpha=0.75)+
    scale_fill_gradientn(name="Altitude",colours = terrain.colors(100))+
    theme_bw()+coord_equal()+xlab("Longitude")+ylab("Latitude")

Digital Elevation Model in ggplot2

Determining Sun Position

In order to make the hill-shading work in an intuitive way, we want to be able to plot the shading on the hills based on a date and time we select. Fortunately, someone has spent the time creating a function for us. A quick google search finds a useful function on StackOverflow. I’ll spare you the details here, but if you’d like to see the function, you can find it here: sunPosition.R

Plotting Hillshade

Next we want to layer on the way the light hits the mountains to see where the sun would strike them given the position of the sun. This can be done using the hillShade function within the raster R package.

test<-function(mon,dy,hr){
  sunpos<-sunPosition(year = 2014,month = mon,day = dy,hour = hr+7,min = 40,sec = 0,lat = mean(DEM$Y),long=mean(DEM$X))
  hs<-hillShade(slope,aspect,angle = sunpos$elevation,direction = sunpos$azimuth,normalize = TRUE)/255.0
  #plot(hs,col=grey(1:100/100),legend=F)
  #plot(dem,col=terrain.colors(100),alpha=0.0,add=T,legend=F)
  HS<-data.frame(rasterToPoints(hs))
  DEM<-data.frame(rasterToPoints(dem))
  b.dem <- seq(min(DEM$layer),max(DEM$layer),length.out=100)

  print(ggplot(HS,aes(x=x,y=y))+
    geom_raster(data=DEM,aes(fill=layer),alpha=0.75)+
    geom_raster(aes(alpha=1-layer),fill="gray20")+
    scale_alpha(guide=FALSE,range = c(0,1.00))+
    scale_fill_gradientn(name="Altitude",colours = terrain.colors(100))+
    theme_bw()+coord_equal()+xlab("Longitude")+ylab("Latitude")+ggtitle(paste0("Canmore, AB Canada - DEM w Lighting for ",mon,"-",dy," @ ",hr,":00")))

}

manipulate(test(mon,dy,hr),
           mon=slider(1,12,initial = month(Sys.Date())),
           dy=slider(1,31,day(Sys.Date())),
           hr=slider(1,23,initial = hour(Sys.time())))

Final Result

All of this put together gives us the following plot:

Alt Text

Check this out in a Shiny app: Digital Elevation Map Lighting

Comments

Python Development with Docker

I've been trying to brush up my skills in Python after spending a lot of time in R. One of the challenges I've run into is maintaining a common development environment on my Mac. Enter Docker.

Force Carbonation Charts

Force carbonation charts are readily available in either metric or imperial units, mostly the latter. The problem was, I use Celcius, and my regulator is in psi. This was a quick solution to the problem.

Leap Year Birthday?

If you're born on Feb 29th, when's your Birthday? Good question.

Other ways to reach me