Conservative regridding with GDAL (?)

news
code
Author

Michael Sumner

Published

December 11, 2024

Can GDAL do conservative re-gridding? For cases of regular grid to regular grid, yes I think it can.

Please note that I’m using tools I’m comfortable with, because I wrote them. I will reframe in other tools and other languages in time. For some reason I’d been blocked on understanding this issue.

Simple grid with four values

Take a grid, 2x2 with values 1,2,3,4 that sums to 10 and warp it to a new size.

## target
dm <- c(2, 2)

g1 <- matrix(c(1, 2, 3, 4), dm[2], dm[1])
e1 <- c(-1, 1, -1, 1)

library(vapour)
library(terra)
terra 1.7.83
## we need this so we can use MEM: via dsn::mem()
Sys.setenv(GDAL_MEM_ENABLE_OPEN = "YES")


## dsn::mem generates a Memory raster, 
## gdal_raster_data is the warper, here we get identity
sg <- gdal_raster_data(dsn::mem(g1, extent = e1, projection = "EPSG:4326"))

unique(sg[[1]])
[1] 1 2 3 4
sum(sg[[1]])
[1] 10
## spatialize for easy plotting
tor <- function(x) {
 
  dm <- attr(x, "dimension")[2:1]
  r <- terra::rast(terra::ext(attr(x, "extent")), ncols = dm[1], nrows = dm[2], crs = attr(x, "projection"), 
                   vals = x[[1]])
  r
}

r <- "sum"

##   (always use mem() "live" to avoid the garbage collector, and only for Float64 I'm afraid)


## now use the same extent but reduce pixel size
tg <- gdal_raster_data(dsn::mem(g1, extent = e1, projection = "EPSG:4326"), target_ext = e1,
                         target_dim = dm * 8, resample = r)

str(tg)
List of 1
 $ : num [1:256] 0.0156 0.0156 0.0156 0.0156 0.0156 ...
 - attr(*, "dimension")= num [1:2] 16 16
 - attr(*, "extent")= num [1:4] -1 1 -1 1
 - attr(*, "projection")= chr "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AU"| __truncated__

So we have a lot more values as we resized from 2x2 to 16x16.

It still looks the same, but the overall quantity has been distributed.

plot(tor(tg), main = r)
## draw boundaries on
abline(v = vaster::x_corner(dm * 8, e1), h = vaster::y_corner(dm * 8, e1))

print(sum(tg[[1]]))
[1] 10
## our sum of 10 was distributed across 256 pixels
unique(tg[[1]])
[1] 0.015625 0.031250 0.046875 0.062500

Now a different example, an actual map projection change.

But, gawd … there’s a bug here. WIP

## now try reprojecting our unit 1x1 longlat grid to LAEA centred a few degrees east and north
crs <- "+proj=laea +lon_0=5 +lat_0=0"
## reproject the extent, we use a densified boundary to find the new extent
prex <- reproj::reproj_extent(c(-1, 1, -1, 1), crs, source = "EPSG:4326")
newg <- gdal_raster_data(dsn::mem(g1, extent = e1, projection = "EPSG:4326"), 
                         target_ext = prex, target_res = c(10000, 10000), target_crs = crs, resample = "sum", options = "-tap")
plot(tor(newg))

sum(newg[[1]])
[1] 10

And what about trying to return, of course this can only be approximate (I want to know how this compares to other tools.)

## what if we go the other way, 

newdm <- attr(newg, "dimension")
newext <- attr(newg, "extent")
x <- gdal_raster_data(dsn::mem(matrix(newg[[1]],newdm[2], byrow = TRUE),  extent = newext, projection = crs), 
                 target_dim = dm, target_ext = e1, target_crs = "EPSG:4326", resample = "sum")
plot(tor(x))

sum(x[[1]])
[1] 9.496864