Skip to content

Instantly share code, notes, and snippets.

@brfitzpatrick
Created September 13, 2018 13:47

Preparing Elevation Data for Visualisation with rayshader

library(magrittr)
library(viridis)
library(raster)
library(rayshader)

Using the DEM available here which the source states is using the LV03 (a.k.a. EPSG 21781) coordinate reference system:

dem.200m.rast <- raster('~/Downloads/data/DHM200.asc')

crs(dem.200m.rast) <- crs('+init=epsg:21781')

print(dem.200m.rast)
## class       : RasterLayer 
## dimensions  : 1201, 1926, 2313126  (nrow, ncol, ncell)
## resolution  : 200, 200  (x, y)
## extent      : 479900, 865100, 61900, 302100  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:21781 +proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +towgs84=674.4,15.1,405.3,0,0,0,0 +units=m +no_defs 
## data source : /home/ben/Downloads/data/DHM200.asc 
## names       : DHM200

Thus is appears that cells in this raster are squares.

Visualising the elevation raster with raster::plot

plot(dem.200m.rast, col = viridis(n = 1e3))

plot of chunk unnamed-chunk-4

My extraction of the data from the raster to a matrix results in a distorted image despite attempting to set the aspect ratio to 1:

dem.200m.mat <- raster::as.matrix(dem.200m.rast)

apply(X = dem.200m.mat, MARGIN = 2, FUN = rev) %>%
  t() %>%
    graphics::image(x = ., col = viridis(n = 1e3), asp = 1)

plot of chunk unnamed-chunk-5

dem.200m.mat.t <- t(dem.200m.mat)

system.time(
  rayshader::sphere_shade(dem.200m.mat.t, texture = 'bw') %>%
  rayshader::plot_map(asp = 1)
)

plot of chunk unnamed-chunk-6

##    user  system elapsed 
##   3.203   0.260   3.463
@mdsumner
Copy link

mdsumner commented Sep 14, 2018

Oh I see, the matrix is implicitly mapped to 0, 1 in image, I just never would do that bare matrix vis. rayshader also needs to keep the original axes but I don't yet understand why it's not 0, ncol; 0, nrow already

@brfitzpatrick
Copy link
Author

Thanks again for looking at this. In my first go at this I plotted the DEM with raster::plot then tried rayshader. The image plot was just a diagnostic to see where in the process the distortion was occurring. Could it be that raster::plot is actually using an aspect ratio different from one that is estimated from the CRS and that I need to supply this to image or rayshader?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment