The goal of wmts is to obtain imagery from online image servers. You can
provide any URL with a GetCapabilities.xml
and quickly get an image
for any region at any resolution.
Any URL? Only Pseudo Mercator services tested so far, see examples below for success stories. Please file an issue or ping me on twitter if you have wants/needs/suggestions/bug-fixes/high-fives.
Use loc
to specify a location (cbind(lon, lat)
) and buffer
to give
a radius around that point in metres. Alternatively, any spatial object
can be used as loc
and its extent (in the Mercator tile server
projection) will be used.
There is an argument max_tiles
which is 9 by default (3x3 256x256
tiles) and ensures a reasonably small sized result. Increase max_tiles
to obtain greater resolution (squares make sense, i.e. 16, 25, 36, 49, …
but the actual shape depends on loc
and buffer
).
The zoom
argument can be used but please take care, the function will
report the zoom in use and it’s advisable not to increase this much (use
max_tiles to get a known maximum amount of data).
See GDAL WMTS documentation for technical details on this GDAL-based approach.
## remotes::install_github("mdsumner/wmts")
centre <- c(-80.888, 32.332) ## lonlat
radius <- 4000 ## metres
u <- "WMTS:https://basemap.nationalmap.gov/arcgis/rest/services/USGSTopo/MapServer/WMTS/1.0.0/WMTSCapabilities.xml,layer=USGSTopo,tilematrixset=default028mm"
library(wmts)
## note that input can be centre (longlat pt), buffer (km) or a spatial object (ignore buffer)
## output is a raster in World Mercator
x <- wmts(u, centre, buffer = radius) ## zoom determined interactively, don't use a hardcoded zoom
#> zoom: 13
raster::plotRGB(x, interpolate = TRUE)
#f <- system.file("gpkg/nc.gpkg", package = "sf", mustWork = TRUE)
#sf <- sf::read_sf(f)
#x <- wmts(u, sf)
tab <- tibble::as_tibble(raster::as.data.frame(x, xy = TRUE))
names(tab) <- c("x", "y", "red", "green", "blue")
##tab <- dplyr::filter(tab, !is.na(red)) ##... when we have missing values, we should drop them or rgb() will error
tab$hex <- rgb(tab$red, tab$green, tab$blue, maxColorValue = 255)
library(ggplot2)
ggplot(tab, aes(x, y, fill = hex)) +
geom_raster(interpolate = TRUE) +
coord_equal() +
scale_fill_identity()
Get a Tassie example.
ext <- raster::extent(146, 148, -43, -40.5)
u <- "WMTS:https://services.thelist.tas.gov.au/arcgis/rest/services/Basemaps/Topographic/MapServer/WMTS/1.0.0/WMTSCapabilities.xml,layer=Basemaps_Topographic,tilematrixset=default028mm"
img <- wmts(u, loc = ext)
#> zoom: 8
raster::plotRGB(img, interpolate = TRUE)
ortho <- wmts("https://services.thelist.tas.gov.au/arcgis/rest/services/Basemaps/Orthophoto/MapServer/WMTS/1.0.0/WMTSCapabilities.xml", raster::extent(147.125, 147.528, -43.054, -42.717), max_tiles = 16)
#> zoom: 11
raster::plotRGB(ortho, interpolate = TRUE)
## note that we can use an spatial object for the extent, saving much pain
url_shade <- "https://services.thelist.tas.gov.au/arcgis/rest/services/Basemaps/HillshadeGrey/MapServer/WMTS/1.0.0/WMTSCapabilities.xml"
greyshade <- wmts(url_shade, ortho)
#> zoom: 10
raster::plotRGB(greyshade, interpolate = TRUE)
prefix_shade <- "WMTS:https://services.thelist.tas.gov.au/arcgis/rest/services/Basemaps/HillshadeGrey/MapServer/WMTS/1.0.0/WMTSCapabilities.xml"
shade <- wmts(prefix_shade, ortho)
#> zoom: 10
# same: raster::plotRGB(shade, interpolate = TRUE)
prefix_shade <- "WMTS:https://services.thelist.tas.gov.au/arcgis/rest/services/Basemaps/Hillshade/MapServer/WMTS/1.0.0/WMTSCapabilities.xml"
hillshade <- wmts(prefix_shade, ortho)
#> zoom: 10
raster::plotRGB(hillshade)
zlocal <- raster::raster(raster::extent(c(16391003, 16394978, -5301061, -5298768)), crs = raster::projection(ortho))
zz <- wmts(prefix_shade, zlocal)
#> zoom: 14
raster::plotRGB(zz)
TasmapRaster <- "https://services.thelist.tas.gov.au/arcgis/rest/services/Basemaps/TasmapRaster/MapServer/WMTS/1.0.0/WMTSCapabilities.xml"
xx <- wmts(TasmapRaster, zlocal, max_tiles = 16)
#> zoom: 15
raster::plotRGB(xx)
street <- "https://services.thelist.tas.gov.au/arcgis/rest/services/Raster/TTSA/MapServer/WMTS/1.0.0/WMTSCapabilities.xml"
st <- wmts(street, raster::shift(zlocal, 5000, 5000), max_tiles = 36)
#> zoom: 16
raster::plotRGB(st, interpolate = TRUE)
A list of Tasmanian servers is here: https://services.thelist.tas.gov.au/arcgis/rest/services/Basemaps
You want the REST GetCapabilities(.xml), and we have to prefix with “WMTS:” and include the layer name (WIP).
u <- "WMTS:https://gibs.earthdata.nasa.gov/wmts/epsg3857/best/1.0.0/WMTSCapabilities.xml,layer=MODIS_Terra_L3_Land_Surface_Temp_Monthly_Day"
im <- wmts(u, cbind(0, 0), 20037508, zoom = 0)
#> zoom: 0
library(raster)
#> Loading required package: sp
plotRGB(im, interpolate = TRUE)
im2 <- wmts(u, extent(100, 180, -55, 15), max_tiles= 36)
#> zoom: 4
plotRGB(im2, interpolate = TRUE)
Please note that the wmts project is released with a Contributor Code of Conduct. By contributing to this project, you agree to abide by its terms.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.