Skip to content

Instantly share code, notes, and snippets.

@brfitzpatrick
Created September 13, 2018 11:54
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save brfitzpatrick/5934680cd8861e1588a358d61bb1bd9f to your computer and use it in GitHub Desktop.
Save brfitzpatrick/5934680cd8861e1588a358d61bb1bd9f to your computer and use it in GitHub Desktop.

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')

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:

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))

plot of chunk unnamed-chunk-5

So I guess what we're seeing here is that raster::plot is adjusting the aspect ratio of the plot in accordance with the coordinate reference system while there is no such information associated with the matrix for image or rayshader to use (though I don't think either of these have the functionality to estimate such an adjustment from the CRS).

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

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

plot of chunk unnamed-chunk-6

##    user  system elapsed 
##   3.192   0.108   3.300
@mdsumner
Copy link

mdsumner commented Sep 13, 2018

try

... %>%
 graphics::image(x = ., col = viridis(n = 1e3, asp = 1)

and

... %>%
 rayshader::plot_map( asp = 1)

But, the "correct" aspect ratio does depend on the cell sizes of the source data if they are different - the .asc format in use here technically doesn't support uneven cell sizes, but it can be done and some software supports it.

Can you show?

print(dem.200m.rast)

It's unfortunate that this format is used, because other better formats actually can store the CRS - so I'd find an alternative if you can.

Otherwise, my general advice is to stay in raster context as much as possible. I'm experimenting with raster inputs for rayshader 3D plots here, for example: https://github.com/hypertidy/cartilage

@mdsumner
Copy link

Can you possibly share the data, or a link to it? I might be able to explain a bit more with it, and I'd like to see it!

@brfitzpatrick
Copy link
Author

Cheers! I tried to link to the data above it's freely available here: https://shop.swisstopo.admin.ch/en/products/height_models/dhm25200

@mdsumner
Copy link

@brfitzpatrick
Copy link
Author

That's the one. I picked it to make this example reproducible. I have access to a 5m DEM that I'm not allowed share which uses the same CRS and has the projection info included in the file, visualising that (probably aggregated lightly) is the end goal for me at the moment.

@brfitzpatrick
Copy link
Author

brfitzpatrick commented Sep 13, 2018

sorry this is my first time using gistr I can authenticate to create new gists but not to edit an existing gist apparently. I've tried your suggestions here: https://gist.github.com/brfitzpatrick/05bf7e518477a69a815b790864f49a1b

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