GDALWarp applib for project()

[mdsumner]

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). ChunkAndWarpImage scales to out of memory usage but is not currently being used for that atm, I think this is equivalent to current capability but can quickly read overviews from massive remote sources (without workarounds).

(off-topic but relevant, GDALWarp can also be directed to an out of memory data store but does its own determination of the chunking required afaik)

The examples below are reasonably quick for large regions with this PR. It's currently a lot slower in terra 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

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

Example

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

(src <- rast(dsn))
#> class       : SpatRaster 
#> dimensions  : 43200, 86400, 1  (nrow, ncol, nlyr)
#> resolution  : 0.004166667, 0.004166667  (x, y)
#> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : gebco_2022_complete_cog.tif 
#> name        : gebco_2022_complete_cog
project(src, rast())
#> 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


##  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")))
#> 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")
#> 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"))

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

Problems

I cannot get it to work with multiple sources, I will update if I figure it out

The second source doesn't have valid values after warping e.g.

 project(rast(rep(system.file("ex/meuse.tif",  package = "terra"), 2)), "OGC:CRS84")
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 

[rhijmans]

Given that this changes a lot and breaks, could you perhaps implement this as a separate function? That would allow for easy comparison with the current implementation; and perhaps for using one over the other for particular cases.

[mdsumner]

ok for sure, I wanted to flag interest with it working - and getting there was a minor personal accomplishment 😅👌

I've been using it for a long time to read into mem and create files and I can't work without being able to do that with these sources now. One other thing terra doesn't allow here is multiple sources that have different grid, and simply letting GDAL pick the output grid.

I'll think about how that might work, because it's really a single step - define a grid (or lean on gdal suggested) and send in source/s, so a new R function makes sense, and would encapsulate my other #1175 .

[mdsumner]

I've made a version with

• a new internal function warper_by_util() (analogous for warp-case to warper())
• new module method warp_by_util (to be called by project method as alternative pathway)
• I left out resample, align, and mask for now while I explore

https://github.com/mdsumner/terra/tree/mdsumner-warp-by-utils-app1

Examples here: https://gist.github.com/mdsumner/872a167b5b69898e96611768cbd0e1d9

still needs thought for the align/resample/mask options (if they belong), and for filename and driver choice which should be trivial to do - the use-cases I see are just

project(rast(dsn), rast(), method = <>)
project(rast(dsn), rast(), filename =<some.file>, driver = <file>)

but other options could be added by loading up options = '-cutline', 'polygon.gpkg' so will probably take me a while to figure out the context for embedding in terra, thank you

[rhijmans]

Thanks (for what is forthcoming)

One other thing terra doesn't allow here is multiple sources that have different grid, and simply letting GDAL pick the output grid.

It rarely makes much sense to me to let GDAL pick the output raster geometry; but I think it is a good idea to add a project<SpatRasterCollection> (and resample<SpatRasterCollection> method. That should be relatively straightforward.

[mdsumner]

oh right, yes, I had meant to explore that, that covers enough

the thing that does make sense about gdal picking is that you can provide one (or two or three or none) of dim, res, ext, crs, and because some sources don't have a regular grid already (the rectify case, or GCPs, or #1175 and **) - but I'm keen to get this going as you suggested and have learnt a lot about terra internals now!

**gdalwarp will output a reasonable resolution and extent in longlat for any source with geolocation arrays with no options set, for example

[rhijmans]

I said that it might be a good idea to add a project<SpatRasterCollection> (and resample<SpatRasterCollection> method.
This actually exists as terra::impose. But you have to provide a template SpatRaster for the output. I suppose it would be possible to allow for more flexibility.

[mdsumner]

oh cool that's good, I'll explore - GDAL's suggested warp output could be a lot better, and your align model is very cool - impose is a better verb than project I think 👌