gdalio_rayshader_tests.R

# example...
library(gdalio)
library(rayshader)
library(topography)
library(gdalwebsrv)

rotate <- function(x) t(apply(x, 2, rev))


gdalio_rayshader_matrix <- function(dsn, ...) {
  v <- gdalio_data(dsn, ...)
  g <- gdalio_get_default_grid()

  m <- matrix(v[[1]], g$dimension[1])[,g$dimension[2]:1, drop = F]

  rotate(rotate(m))|>
    apply(2,rev)
}

gdalio_rayshader_array <- function(dsn, alpha=1,...) {
  v <- gdalio_data(dsn, ...)
  g <- gdalio_get_default_grid()

  matrix_thing <- function(.v){
    m <- matrix(.v, g$dimension[1])#[,g$dimension[2]:1, drop = F]
    rotate(m)
  }

  v2 <- lapply(v, matrix_thing)
  a <- matrix(NA, g$dimension[2],g$dimension[1])
  aa <- array(c(unlist(v2, use.names = FALSE), a), c(g$dimension[2], g$dimension[1], 4))[,g$dimension[1]:1, , drop = FALSE]
  aa <- aa%>%
    scales::rescale(.,to=c(0,1))
  aa[,,4] <- alpha
  return(aa)
}

# Various gdalio grid setting options...
gdalio_set_default_grid(list(extent=c(-180,  180,  -90,   90),
                             dimension=c(180*4, 90*4),
                             projection=c("+proj=longlat +datum=WGS84")))

gdalio_set_default_grid(list(extent=c(-1e8, 1e8, -5e7, 5e7),
                             dimension=c(180*4,90*4),
                             projection=c('+proj=laea +lon_0=-90')))

# gdalio_set_default_grid(list(extent = c(0.0000, 849024.0785, 0.0000, 4815054.8210),
#                              dimension = c(3024, 3024),
#                              projection = sf::st_crs(54030)))

# grid0 <- list(extent = c(-1e6, 1e6, -5e5, 5e5 ),
#               dimension = c(512, 256),
#               projection = "+proj=laea +lon_0=147 +lat_0=-42")
# gdalio_set_default_grid(grid0)


gdalio_set_default_grid(list(extent=c(-2e7,  2e7, -8e6,  8e6 ),
                             dimension=c(180*4, 90*4),
                             projection=c('+proj=mbtfpq')))

# For comparison with known working function
arr <- gdalio_graphics(server_file("wms_virtualearth"))
grid1 <- gdalio_get_default_grid()
plot(matrix(grid1$extent, 2), type = "n", asp = 1)
graphics::rasterImage(arr, grid1$extent[1], grid1$extent[3], grid1$extent[2], grid1$extent[4])


# DO some rayshader stuff
topo_mat <- gdalio_rayshader_matrix(topography_source('aws'))

drape_array <- gdalio_rayshader_array(server_file("wms_virtualearth"),
                                      alpha=0.7, bands=1:3)
plot_map(drape_array)
topo_mat %>%
  sphere_shade()%>%
  add_overlay(., drape_array,rescale_original=T) %>%
  plot_map()