terra_project_issue.md

Why does terra take forever (or never complete ...)?

## remotes::install_github("hypertidy/gdalio")
library(gdalio) 

## warp target 
extent <- c(-2e5, 2.4e5,  -2e5, 2.84e5)
dimension <- c(1024, 1024)
projection <- "+proj=lcc +lon_0=146 +lat_0=-42 +lat_1=-40 +lat_0=-44 +datum=WGS84"
target <- terra::rast(terra::ext(extent), nrows = dimension[2], ncols = dimension[1], crs = projection)

gdalio_set_default_grid(target)

img_src <- "WMTS:https://services.thelist.tas.gov.au/arcgis/rest/services/Basemaps/Topographic/MapServer/WMTS/1.0.0/WMTSCapabilities.xml,layer=Basemaps_Topographic,tilematrixset=default028mm"
## source data
terr <- terra::rast(img_src)

source(system.file("raster_format/raster_format.codeR", package = "gdalio", mustWork = TRUE))
## these should do exactly the same tasks through GDAL warper, but terra doesn't complete
v <- gdalio_terra(img_src, bands = 1:3, band_output_type = "Int32"); terra::plotRGB(v)
## why does this not complete ... don't run it, trust me it takes way too long and probably never finish
#terr <- terra::project(terr, target)

It does work well for an online GeoTIFF with overviews, so perhaps the COG driver does more for the auto-selection of the right zoon than the WMTS driver does. The GDAL WARP app lib does this well, for any driver - and I think terra might not use that.

update

add the zoom_level explicitly to the DSN and terra works (so the COG/GTiff driver is automatically selecting the right overview, and the way vapour/gdalio uses the warper already does this for other drivers)

img_src <- 
"WMTS:https://services.thelist.tas.gov.au/arcgis/rest/services/Basemaps/Topographic/MapServer/WMTS/1.0.0/WMTSCapabilities.xml,layer=Basemaps_Topographic,tilematrixset=default028mm"
## source data
(terr <- terra::rast(paste0(img_src, ",zoom_level=8")))

extent <- c(-2e5, 2.4e5,  -2e5, 2.84e5)
dimension <- dev.size("px")
projection <- "+proj=lcc +lon_0=146 +lat_0=-42 +lat_1=-40 +lat_0=-44 +datum=WGS84"
target <- terra::rast(terra::ext(extent), nrows = dimension[2], ncols = dimension[1], crs = projection)

x <- terra::project(terr, target, method = "bilinear")

terra::plotRGB(x)

REST servers with gdalio

src <-  "http://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer?f=json"

library(gdalio)
source(system.file("raster_format/raster_format.codeR", package = "gdalio", mustWork = TRUE))

x <- gdalio_terra(src, band_output_type = "integer", bands = 1:3)
terra::plotRGB(x)

## set your grid
v <- terra::vect(system.file("ex/lux.shp", package="terra"))
gdalio_set_default_grid(terra::rast(v, res = 0.001))

x <- gdalio_terra(src, band_output_type = "integer", bands = 1:3)
terra::plotRGB(x)

## set your grid
v1 <- terra::project(v, "EPSG:3857")
gdalio_set_default_grid(terra::rast(v1, nrows = 512, ncols = 768))
x1 <- gdalio_terra(src, band_output_type = "integer", bands = 1:3)
terra::plotRGB(x1)