hypertidy / cartilage

Adventures in rgl coordinates with raster and rayshader, WIP
10 stars 0 forks source link

overlay example #2

Open mdsumner opened 6 years ago

mdsumner commented 6 years ago

The TAS 25m DEM with list transport segments

<bump!> added some comments

## geocoding function to find a location by address
address <- function(x) {
  dismo::geocode(x)
}
## build a small raster extent from an address
address_extent <- function(x) {
  addr <- address(x)
  out <- raster::raster(raster::extent(unlist(addr[c("xmin", "xmax", "ymin", "ymax")])), 
                 crs = "+init=epsg:4326")
  out
}
## buffer out from a an address for a target raster 
address_radius <- function(x, target, radius) {
  ras <- raster::projectExtent(address_extent(x), raster::projection(target))
  out <- raster::raster(spex::buffer_extent(raster::extent(ras) + radius, raster::res(target)[1]), crs = raster::projection(ras))
  raster::res(out) <- raster::res(target)
  out
}

library(cartilage)

## read in data sets and get a location by address

## vector lines
roads <- sf::read_sf(".../List\\GISdata\\Roads\\transport_segments.shp")
roads <- sf::st_crop(roads, sf::st_bbox(raster::extent(a)))

## raster relief
f <- ".../tasDEM_2015.tif"
r <- raster::raster(f)

## a location, subset of the raster relief
a <- address_radius("<somewhere in Tasmania>", r, radius = 10000)

## crop to the address region and plot via rayshader in rgl
rr <- raster::crop(r, raster::extent(a))
ambient(rr)

## convert to 3D-ready format and copy Z on to the road segment vectors (+ a bit to avoid z-fighting)
library(anglr)  #hypertidy/anglr
library(silicate) #hypertidy/silicate
sc <- anglr::copy_down(silicate::SC(roads), rr + 5)

## unique colours for a vector attribute
cols <- sf:::sf.colors(length(unique(sc$object$USER_TYPE)))
## copy colours onto object level, to avoid default scheme
sc$object$color_ <- cols[as.integer(factor(sc$object$USER_TYPE))]

## add the vectors to the scene using anglr
plot3d(sc, add = TRUE)

## modify aspect ratio (especially if relief in metres, maps in longlat)
rgl::aspect3d(1, 1, .2)

image

mdsumner commented 6 years ago

image

jimjam-slam commented 6 years ago

I was playing with rayshader a couple of weeks ago and was thinking that it would be very cool to overlay other layers onto the 3D surface, so it's nice to know other people are thinking the same thing 😄