stuff to fix from readme

[mdsumner]

Get imagery and DEM for use in 3D visualization

This code downloads a specific region as elevation and imagery at specific zoom levels to build a 3D scene.

library(ceramic); library(quadmesh)
library(raster); library(rgl);
library(reproj); library(htmlwidgets)
clear3d()

## longlat extent
ex0 <- c(147.15, 147.45, -42.9, -42.6)
ex <- extent(ex0)

## local LAEA projection, based on the centre of the extent
prj <- sprintf("+proj=laea +lon_0=%f +lat_0=%f +datum=WGS84", mean(ex0[1:2]), mean(ex0[3:4]))
## Mapbox elevation
dem <- cc_elevation(ex, zoom = 8)
## Mapbox satellite imagery
im <- cc_location(ex, zoom = 13)

## quadmesh with texture for rgl, in local projection
qm <- reproj::reproj(quadmesh(dem, texture = im), prj)

## plot with rgl, set the aspect ratio and backround
shade3d(qm, lit = FALSE);
aspect3d(1, 1, .1)
bg3d(grey(0.8))

The zoom levels were chosen by first reading an automatic level for the extents, and then tweaking the zoom - higher resolution for the imagery, lower resolution for the elevation. This provides a compelling visualization in 3D as the imagery is textured onto the elevation data, using rgl's mesh3d type and shade3d() function.

Local caching of tiles

A key feature of ceramic is caching, all data is downloaded in a systematic way that is suitable for later re-use. Many tools for imagery services treat the imagery as transient, but here we take control over the raw data itself. All file names match exactly the address URL of the original source data.

There is a helper function to find existing tiles.

aa <- cc_location(loc = cbind(0, 0), buffer = 330000, type = "mapbox.satellite")
ceramic_tiles(zoom = 7, type = "mapbox.satellite")

and every row has the extent values useable directly by raster:

ceramic_tiles(zoom = 7, type = "mapbox.satellite") %>% 
  dplyr::slice(1:5) %>% 
   purrr::transpose()  %>% 
  purrr::map(~raster::extent(unlist(.x[c("xmin", "xmax", "ymin", "ymax")])))

Another example

my_bbox <-
  st_bbox(c(xmin = 144,
            xmax = 147.99,
            ymin = -44.12,
            ymax = -40),
          crs = st_crs("+proj=longlat +ellps=WGS84"))
im <- cc_location(cbind(145.5, -42.2), buffer = 5e5)
plotRGB(im)
plot(st_transform(ozmaps::abs_lga$geometry, projection(im)), add = TRUE, lwd = 2, border = "white")

An internal function sets up a plot of tiles at particular zoom levels.

ceramic::plot_tiles(ceramic_tiles(zoom = c(7, 9)))

And we can add the tiles to an existing plot.

plotRGB(im)
ceramic::plot_tiles(ceramic_tiles(zoom = 7), add = TRUE)

Future improvements

See the Issue tab and please make suggestions and give feedback!

There's an ongoing process to uncouple from the "spatial packages", currently we have an unclassy mix of sp, terra, raster, and even sf
in the readme. We'll fix it.