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.
## 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 identitysg <-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 plottingtor <-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 sizetg <-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 onabline(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 pixelsunique(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 northcrs <-"+proj=laea +lon_0=5 +lat_0=0"## reproject the extent, we use a densified boundary to find the new extentprex <- 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))