tylermorganwall / rayshader

R Package for 2D and 3D mapping and data visualization
https://www.rayshader.com/
2.06k stars 214 forks source link

raster_to_matrix() not working with terra SpatRaster #226

Closed wiesehahn closed 2 years ago

wiesehahn commented 2 years ago

When trying to run elmat = raster_to_matrix(dtm_terra) on a terra raster I get following error

Error in h(simpleError(msg, call)) : 
  error in evaluating the argument 'y' in selecting a method for function 'extract': unable to find an inherited method for function ‘extent’ for signature ‘"SpatRaster"’

with elmat = raster_to_matrix(raster::raster(dtm_terra)) it works

Is it the wanted behavior to not run with terra SpatRasters?

tylermorganwall commented 2 years ago

Right now that function specifically works with raster objects with the raster package, but I could look into adding terra support.

tylermorganwall commented 2 years ago

Check out the terra_support branch and see if that fixes your error.

https://github.com/tylermorganwall/rayshader/tree/terra_support

wiesehahn commented 2 years ago

Hey thanks for your effort. I am not that comfortable with R package development and git yet, but if I checked it correctly there is still a problem with your code. Running it on a sample dataset as follows

library(terra)
f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)

mat <- rayshader::raster_to_matrix(raster = r)

I get an error

Error in h(simpleError(msg, call)) : 
error in evaluating the argument 'y' in selecting a method for function 'extract': unable to find an inherited method for function ‘extent’ for signature ‘"SpatRaster"’

I guess adding something like the following would work

  return(matrix(terra::extract(raster, terra::ext(raster)), 
                nrow = ncol(raster), ncol = nrow(raster)))
tylermorganwall commented 2 years ago

Try the latest update to the branch, I believe it should work now.

wiesehahn commented 2 years ago

Yes, it works! Nice. As terra is the replacement for raster this should be a step in the right direction. Although I have not tested it, it should be faster with terra ([https://tmieno2.github.io/R-as-GIS-for-Economists/extract-speed.html]())

This also makes it possible to "stream" COG data (if available) without the step to download the raster. Although the data still has to be read in memory for matrix conversion. e.g.

cog.url <- "/vsicurl/https://srsp-open-data.s3.eu-west-2.amazonaws.com/lidar/phase-1/dtm/27700/gridded/NS77_1M_DTM_PHASE1.tif"

ras <- terra::rast(cog.url)
ras <- terra::crop(ras, terra::ext(ras)/10)

mat <- raster_to_matrix(raster = ras)

mat %>%
  sphere_shade(texture = "imhof1", zscale = 0.2) %>%
  plot_map()
mdsumner commented 2 years ago

nice, can go even further with the warper, and do

info <- vapour::vapour_raster_info(cog.url)
## not entirely sure what ext/10 does but it's something like this
ex10 <- matrix(info$extent, ncol = 2, byrow = TRUE)
ex10 <- rep(rowMeans(ex), each = 2) + rep(apply(ex, 1, diff)/10, each = 2) * c(-1, 1, -1, 1)
mat <- matrix(vapour::vapour_warp_raster_dbl(cog.url, dimension = as.integer(info$dimXY/10), extent = ex10, projection = info$projection), as.integer(info$dimXY[2]/10))

you can use any values there for dimension, extent, projection - but obviously a cap on the dimension is advisable, and it should geographically overlap (or not much use!), but the warper handles any CRS, with some configurations for weird situations etc.

use resample for controlling the sample (e.g. "bilinear", default is "nearest"

This also works via terra of course, but it won't work so well with non-COGs such as image tile servers. This is a slow COG for large regions (with decimation), as it doesn't have overviews.

To see what the speed can be like for a large region try a COG with overviews:

## using the same info as the cog.url to spec the matrix
cog <- "/vsicurl/https://public.services.aad.gov.au/datasets/science/GEBCO_2019_GEOTIFF/GEBCO_2019.tif"
mat <- matrix(vapour::vapour_warp_raster_dbl(cog, dimension = as.integer(info$dimXY/10), extent = info$extent, projection = info$projection, resample = "bilinear"), as.integer(info$dimXY[2]/10))

(obviously the data in the UK cog is better)

tylermorganwall commented 2 years ago

Terra support added in 2651da4d1faa20cbf88ba4544f9f164738370adb