Skip to content

Instantly share code, notes, and snippets.

@mdsumner
Last active January 20, 2025 13:46
Show Gist options
  • Save mdsumner/cceef5bef48694905c22348c7436f869 to your computer and use it in GitHub Desktop.
Save mdsumner/cceef5bef48694905c22348c7436f869 to your computer and use it in GitHub Desktop.

WIP, trying to get the right resolution from a web tile server

library(tmap)
library(terra)


## reliable world polygons
psrc <- "/vsizip//vsicurl/https://github.com/wmgeolab/geoBoundaries/raw/main/releaseData/CGAZ/geoBoundariesCGAZ_ADM0.zip"
psrcsql <- "SELECT shapeGroup FROM geoBoundariesCGAZ_ADM0 WHERE shapeGroup IN ('AUS')"
d <- vect(psrc, query = psrcsql)
prj <- "EPSG:3577"

d <- project(d, prj)
ex <- c(xmin = 940936, xmax = 1560393, ymin = -4936503, ymax = -4356586)
e <- ext(ex)
d <- crop(d, e)


#tm_basemap() + tm_rgb() +  tm_shape(d, main = T) + tm_lines()
tm_shape(d, main = T) + tm_lines()

## what if we want an image server version
## openstreetmap tiles: 
dsn <- "<GDAL_WMS><Service name=\"TMS\"><ServerUrl>https://tile.openstreetmap.org/${z}/${x}/${y}.png</ServerUrl></Service><DataWindow><UpperLeftX>-20037508.34</UpperLeftX><UpperLeftY>20037508.34</UpperLeftY><LowerRightX>20037508.34</LowerRightX><LowerRightY>-20037508.34</LowerRightY><TileLevel>19</TileLevel><TileCountX>1</TileCountX><TileCountY>1</TileCountY><YOrigin>top</YOrigin></DataWindow><Projection>EPSG:3857</Projection><BlockSizeX>256</BlockSizeX><BlockSizeY>256</BlockSizeY><BandsCount>3</BandsCount><UserAgent>R user</UserAgent><Cache /></GDAL_WMS>"
#dsn <- "<GDAL_WMS><!-- Data is subject to term of use detailed at http://code.google.com/intl/nl/apis/maps/terms.html andhttp://www.google.com/intl/en_ALL/help/terms_maps.html --><Service name=\"TMS\"><!-- ServerUrl>http://mt.google.com/vt/lyrs=m&amp;x=${x}&amp;y=${y}&amp;z=${z}</ServerUrl> --><!-- Map --><!-- <ServerUrl>http://mt.google.com/vt/lyrs=s&amp;x=${x}&amp;y=${y}&amp;z=${z}</ServerUrl> --> <!-- Satellite --><ServerUrl>http://mt.google.com/vt/lyrs=y&amp;x=${x}&amp;y=${y}&amp;z=${z}</ServerUrl> <!-- Hybrid --><!-- <ServerUrl>http://mt.google.com/vt/lyrs=t&amp;x=${x}&amp;y=${y}&amp;z=${z}</ServerUrl> --> <!-- Terrain --><!-- <ServerUrl>http://mt.google.com/vt/lyrs=p&amp;x=${x}&amp;y=${y}&amp;z=${z}</ServerUrl> --> <!-- Terrain, Streets and Water  --></Service><DataWindow><UpperLeftX>-20037508.34</UpperLeftX><UpperLeftY>20037508.34</UpperLeftY><LowerRightX>20037508.34</LowerRightX><LowerRightY>-20037508.34</LowerRightY><TileLevel>20</TileLevel><TileCountX>1</TileCountX><TileCountY>1</TileCountY><YOrigin>top</YOrigin></DataWindow><Projection>EPSG:900913</Projection><BlockSizeX>256</BlockSizeX><BlockSizeY>256</BlockSizeY><BandsCount>3</BandsCount><MaxConnections>5</MaxConnections><Cache /></GDAL_WMS>"
tree <- grid::grid.ls(viewports = TRUE, flatten = TRUE, grobs = FALSE)
tree$name
vpname <- grep("vp_map", tree$name, value = TRUE)[1]  ## how to get the right vp ... but this is the one we want
grid::seekViewport(vpname)
vp <- grid::current.viewport()
## not sure how to reset our seek viewport ...


pix <- lapply(grid::deviceLoc(grid::unit(c(0, 1), "npc"), grid::unit(c(0, 1), "npc"), device = T), as.numeric)
## this is the "ideal dimension" for a background because this is the sea of pixels we have
dm <- unlist(lapply(pix, \(.x) abs(ceiling(diff(.x)))))
ex <- c(vp$xscale, vp$yscale)

library(terra)
target <- rast(ext(ex), ncols = dm[1], nrows = dm[2], crs = prj)
## project multiband by util (without it we can't hit an overview without extra work, project() on its on targets zoom=0)
project_multi <- function(x, target, ...) {
  project(x, rep(target, length.out = terra::nlyr(x)), by_util = TRUE, ...)
}
im <- project_multi(rast(dsn), target, method = "cubic")
tm_shape(im) + tm_rgb() + tm_shape(d) + tm_lines()

grid::grid.lines(pix$x[1], c(0, dev.size("px")[2]), default.units = "native")
grid::grid.lines(pix$x[2], c(0, dev.size("px")[2]), default.units = "native")
grid::grid.lines(c(0, dev.size("px")[1]), pix$y[1], default.units = "native")
grid::grid.lines(c(0, dev.size("px")[1]), pix$y[2], default.units = "native")

image

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment