py rasterio from an R user perspective

reprojecting-geotiff-zoom.md

working on the full-affine approach to warping a zoom, from the rasterio quickstart

https://rasterio.readthedocs.io/en/latest/topics/reproject.html#reprojecting-a-geotiff-dataset

this is very rough, just a bit of exploration and reference for me to expand elsewhere ...

  f <- '/vsicurl/https://github.com/rasterio/rasterio/raw/master/tests/data/RGB.byte.tif'
ri <- vapour::vapour_raster_info(f)

g <- list(extent = ri$extent, dimension = ri$dimXY, projection = ri$projection)

affine_from_gdal <- function(x) {
  ## see a,b,c,d,e,f https://github.com/rasterio/affine#usage
  idx <- c(2L, 3, 1, 5, 6, 4)  ## c, a, b, f, d, e
  ## https://github.com/rasterio/affine/blob/78c20a0cfbb5ece3dfadc1a5550465936a4fa731/affine/__init__.py#L178
  matrix(c(x[idx], 0, 0, 1), nrow = 3L, byrow = TRUE)
}

translation <- function(dx, dy) {
  matrix(c(1, 0, dx, 
           0, 1, dy, 
           0, 0, 1), nrow = 3L, byrow = TRUE)
}
scale <- function(x) {
  matrix(c(x, 0, 0, 
           0, x, 0, 
           0, 0, 1), nrow = 3L, byrow = TRUE)
}
translation(-ri$dimXY[1] / 2, -ri$dimXY[2] / 2)
#>      [,1] [,2]   [,3]
#> [1,]    1    0 -395.5
#> [2,]    0    1 -359.0
#> [3,]    0    0    1.0
scale(2)
#>      [,1] [,2] [,3]
#> [1,]    2    0    0
#> [2,]    0    2    0
#> [3,]    0    0    1
## phew c(0, 0, 1) is col, row, 1 in affine transformation form
affine_from_gdal(ri$geotransform) %*% c(0, 0, 1) == c(ri$extent[c(1, 4)], 1)
#>      [,1]
#> [1,] TRUE
#> [2,] TRUE
#> [3,] TRUE
## what about xmax, ymin
affine_from_gdal(ri$geotransform) %*% c(ri$dimXY[1], ri$dimXY[2], 1) == c(ri$extent[c(2, 3)], 1)
#>      [,1]
#> [1,] TRUE
#> [2,] TRUE
#> [3,] TRUE
## so, zoom out is
#src_transform*A.translation(-src.width/2.0, -src.height/2.0)*A.scale(2.0)
zoom2 <- affine_from_gdal(ri$geotransform) %*% translation(-ri$dimXY[1] / 2, -ri$dimXY[2] / 2) %*% scale(2) 

gdal_from_affine <- function(x) {
  t(x[1:2, ])[c(3, 1, 2, 6, 4, 5)]
}

gdal_from_affine(zoom2)
#> [1]  -16680.0000     600.0759       0.0000 2934630.0000       0.0000
#> [6]    -600.0836
g$extent <- affinity::gt_dim_to_extent(gdal_from_affine(zoom2), ri$dimXY)
library(gdalio)
gdalio_set_default_grid(g) 
source(system.file("raster_format/raster_format.codeR", package = "gdalio", mustWork = TRUE))
im <- gdalio_terra(f, bands = 1:3)
terra::plotRGB(im)

^{Created on 2021-12-07 by the reprex package (v2.0.1)}

the above approach is more general, but exactly what's happening is easily done simpler for a cascade of zooms like

f <- '/vsicurl/https://github.com/rasterio/rasterio/raw/master/tests/data/RGB.byte.tif'

library(terra)

r <- rast(f)
sq <- seq(1, 2.2, length.out = 4)
par(mfrow = c(3, 3))
for (i in 2^rev(sq)) {
  r1 <- crop(rast(r), ext(r) / i)
  res(r1) <- res(r) / i
  plotRGB(project(r, r1, method = "bilinear"))
}
plotRGB(r)
for (i in 2^sq) {
  r1 <- extend(rast(r), ext(r)  * i)
  res(r1) <- res(r) * i
  plotRGB(project(r, r1, method = "near"))
  plot(ext(r), add = TRUE)
}

not making any big point, just I happened to try this