resample() function in R's "terra" library does not work properly in specific situati-CodePudding

In the code below, I am trying to resample a high resolution raster into a lower resolution raster using resample(x, y, method = "sum"). However, the resampled raster shows NA on some edges.

library(terra)
set.seed(42)

low_res <- rast(xmin = -1.05, xmax = 1.05, ymin = -0.05, ymax = 2.05, res = 0.5)

high_res <- rast(xmin = -1, xmax = 1, ymin = 0, ymax = 2, res = 0.01)
high_res[] <- runif(ncell(high_res))
plot(high_res, colNA = "darkblue")

resampled <- resample(high_res, low_res, method = "sum")
plot(resampled, colNA = "darkblue")
plot(as.polygons(low_res), add=TRUE, border='black', lwd=1)

The high resolution raster:

The resampled raster (the dark-blue cells are NAs):

But, if the extent of the low resolution raster rounded (i.e., deleting the _.05), everything looks good:

library(terra)
set.seed(42)

##################################
# only changed extent here
low_res <- rast(xmin = -1, xmax = 1, ymin = -0, ymax = 2, res = 0.5) 
##################################

high_res <- rast(xmin = -1, xmax = 1, ymin = 0, ymax = 2, res = 0.01)
high_res[] <- runif(ncell(high_res))
plot(high_res, colNA = "darkblue")

resampled <- resample(high_res, low_res, method = "sum")
plot(resampled, colNA = "darkblue")
plot(as.polygons(low_res), add=TRUE, border='black', lwd=1)

The resampled raster:

CodePudding user response：

Okay, I think, if a cell in the low_res raster is not within the extent of the high_res raster, it will get NA value (I wish it was not like this). So, as a solution, I tried adding cells around the high_res raster using extend() (which, I think, the conservative size for extending high_res would be the size of the resolution of low_res raster). The code below works well:

library(terra)

set.seed(42)
low_res <- rast(xmin = -1.05, xmax = 1.05, ymin = -0.05, ymax = 2.05, res = 0.5)

high_res <- rast(xmin = -1, xmax = 1, ymin = 0, ymax = 2, res = 0.01)
high_res[] <- runif(ncell(high_res))
plot(high_res, colNA = "darkblue")

#######################################
# add paddings around high_res raster
# with at least the resolution size of low_res raster
add_n_cells <- ceiling(res(low_res) / res(high_res))
high_res <- extend(high_res, add_n_cells)
plot(high_res, colNA = "darkblue")
#######################################

resampled <- resample(high_res, low_res, method = "sum")
plot(resampled, colNA = "darkblue")
plot(as.polygons(low_res), add=TRUE, border='black', lwd=1)

High resolution raster extended:

Final resampled raster:

CodePudding user response：

What you experience is due to a bug in GDAL that has been fixed.

I do not see the problem (in your first example) with GDAL 3.6

library(terra)
gdal()
#[1] "3.6.0"

Here is a post about GDAL versions with R packages.

Right now, you can get GDAL version 3.6.2 on windows by installing from source in R-devel and RTools43 and on OSX by installing from source with the current versions of brew.