project raster via GDALWarp() app lib

[mdsumner]

(rewritten PR to get past merge issue, previous were #1120 and #1216)

Propose to use GDALWarp utility library for project() raster rather than with ChunkAndWarpImage.

The benefit is automatic choice of a reasonable overview for multi-zoom sources (like large GeoTIFFs or image tile servers), (and hence very fast read, with data resolution suitable for target raster even over large regions.

The problem remains that we only access bands from the first source. see below

(This doesn't address the capability for collections, but I think it's worth sharing as is now because that will require a solution for the multiple sources, and maybe that can just be done at the R level).

Currently is a lot slower in terra for online image servers and very large COGs because there's no overview detection and it defaults to having to scan the highest resolution tiles, for small local regions the current capability is fast because only a small number of tiles are visited.

The PR does the following

• replaces ChunkAndWarp with GDALWarp()
• no longer uses set_warp_options(), I've left the working inline in warper()
• sets WRITE_FLUSH, and INIT_DEST, and NUM_THREADS as before
• adds a warp_by_util/warper_by_util module
• modifies project() to add a by_util argument to allow invoking this pathway from R

GDALwarp() was librarified in 2.1.0 so this should be supported by all terra installs.

Problem

I don't think GDALWarp() can work this way for multi-source SpatRaster (at least, not until version 3.7.0), because with multiple sources you are writing to specific bands in the destination dataset. In 3.7.0 you can explicitly set up 'srcbands' and 'dstbands' as options for the GDALWarpAppOptions. I don't think this can be achieved until then without reworking how the multple bands in a multi-source SpatRaster are collated.

So we could set project() to warn or error on multi-source inputs, which I suggest is a good way to go. I think it would be rare to generate multi-source rast objects from very large online sources.

Examples

## GEBCO 2023 thanks to Philippe Massicotte
dsn <- "/vsicurl/https://gebco2022.s3.valeria.science/gebco_2022_complete_cog.tif"
library(terra)
#> terra 1.7.41

## basic
project(rast(system.file("ex/meuse.tif",  package = "terra")), "OGC:CRS84", by_util = TRUE)
#> class       : SpatRaster 
#> dimensions  : 95, 105, 1  (nrow, ncol, nlyr)
#> resolution  : 0.0004381799, 0.0004381799  (x, y)
#> extent      : 5.720655, 5.766664, 50.95453, 50.99616  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
#> source(s)   : memory
#> name        : meuse 
#> min value   :   138 
#> max value   :  1736
project(rast(system.file("ex/meuse.tif",  package = "terra")), "+proj=laea +lon_0=5 +lat_0=50", by_util = TRUE)
#> class       : SpatRaster 
#> dimensions  : 115, 81, 1  (nrow, ncol, nlyr)
#> resolution  : 40.00032, 40.00032  (x, y)
#> extent      : 50615.37, 53855.39, 106469.4, 111069.5  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=laea +lat_0=50 +lon_0=5 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#> source(s)   : memory
#> name        : meuse 
#> min value   :   138 
#> max value   :  1736

## big source (fast because the warper (and not just the COG driver) auto-picks the right overview for the target)
src <- rast(dsn)
project(src, rast(), by_util = TRUE)
#> class       : SpatRaster 
#> dimensions  : 180, 360, 1  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
#> source(s)   : memory
#> name        : gebco_2022_complete_cog 
#> min value   :                   -8106 
#> max value   :                    5528

## align works
src <- rast(dsn)
project(src, rast(ext(-180, 0, -90, 0), ncols = 90, nrows = 45), align = TRUE, by_util = TRUE)
#> class       : SpatRaster 
#> dimensions  : 90, 180, 1  (nrow, ncol, nlyr)
#> resolution  : 2, 2  (x, y)
#> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
#> source(s)   : memory
#> name        : lyr.1 
#> min value   : -7159 
#> max value   :  5500

## mask works as per project 
plot(project(src, rast(ext(-20037508.34, 20037508.34, -1e7, 1e7 ), ncols = 512, nrows = 256, crs = "+proj=sinu"), mask = TRUE, by_util = TRUE))

##  something more interesting
(laea <- project(src, rast(ext(c(-1, 1, -1, 1) * 1e6), res = 5000, crs = "+proj=laea +lon_0=147 +lat_0=-42"), by_util = TRUE))
#> class       : SpatRaster 
#> dimensions  : 400, 400, 1  (nrow, ncol, nlyr)
#> resolution  : 5000, 5000  (x, y)
#> extent      : -1e+06, 1e+06, -1e+06, 1e+06  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=laea +lat_0=-42 +lon_0=147 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#> source(s)   : memory
#> name        : gebco_2022_complete_cog 
#> min value   :                   -5948 
#> max value   :                    1968

## another source (Mercator tiles)
imgsrc <- sprintf("<GDAL_WMS><Service name=\"TMS\"><ServerUrl>http://services.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer/tile/${z}/${y}/${x}</ServerUrl></Service><DataWindow><UpperLeftX>-20037508.34</UpperLeftX><UpperLeftY>20037508.34</UpperLeftY><LowerRightX>20037508.34</LowerRightX><LowerRightY>-20037508.34</LowerRightY><TileLevel>17</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>10</MaxConnections><Cache /><UserAgent>%s</UserAgent></GDAL_WMS>", 
                                    getOption("HTTPUserAgent"))

project(rast(imgsrc), rast(res = 0.5), method = "bilinear", by_util = TRUE)
#> class       : SpatRaster 
#> dimensions  : 360, 720, 3  (nrow, ncol, nlyr)
#> resolution  : 0.5, 0.5  (x, y)
#> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
#> source(s)   : memory
#> names       : GDAL_WMS>_1, GDAL_WMS>_2, GDAL_WMS>_3 
#> min values  :          61,         120,         107 
#> max values  :         255,         255,         255


meuse <- rast(system.file("ex/meuse.tif",  package = "terra"))
res(meuse) <- res(meuse)/4
plotRGB(project(rast(imgsrc), meuse, method = "cubic", by_util = TRUE))

## sadly, multi DOES NOT WORK (because we can't line up arbitrary bands in the GDALwarp destination object)
project(rast(rep(system.file("ex/meuse.tif",  package = "terra"), 2)), "OGC:CRS84", by_util = TRUE)
#> class       : SpatRaster 
#> dimensions  : 95, 105, 2  (nrow, ncol, nlyr)
#> resolution  : 0.0004381799, 0.0004381799  (x, y)
#> extent      : 5.720655, 5.766664, 50.95453, 50.99616  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
#> source(s)   : memory
#> names       : meuse, meuse 
#> min values  :   138,     0 
#> max values  :  1736,     0

^{Created on 2023-07-06 with reprex v2.0.2}

[mdsumner]

awesome, thanks !