Resample() cuts off pixels when extents are different

[zoeschindler]

Hello, I am not sure if this is a user-side error, but I wanted to report it in case it is not.
I have two rasters, a fine raster with resolution 0.1m (link) and a rough raster with resolution 0.5m (link). I want to resample the fine raster to the rough raster using the maximum value per cell. Also, I want the cells of the new raster to overlap with the cells of the old rough raster. I tried this using resample(), but it cuts off cells at the edges which contain only few cells in the fine raster:

# plot raster
plot(rough, col = "red")
plot(fine, add = T)
  
# maximum resample
resampled <- resample(fine, rough, "max")

# plot raster
plot(resampled , col = "red")

Plotting the rough and fine raster on top of each other:

Plotting the new raster which is missing 4 pixels:

I noticed I don't have this issue when changing the extent of the fine raster previously:

# change extent
new_ext <- ext(floor(xmin(fine)), ceiling(xmax(fine)), floor(ymin(fine)), ceiling(ymax(fine)))
fine <- extend(fine, new_ext)

(Though here, I end up with too many cells, most likely because of #844.)

I don't know if this is intended behaviour or not. I think I was mislead by the documentation of resample() because it says "If the origin and extent of the input and output are the same, you should consider using these other functions instead: aggregate, disagg, extend or crop.", which sounds like resample would be able to handle the different extents.)

[rhijmans]

Here is a reproducible example of this issue that does not depend on files. It shows that cells that cross the boundaries of cells they are resampled to are not treated consistently.

Example data

library(terra)
r <- rast(ncol=5, nrow=5, xmin=0, xmax=1, ymin=0, ymax=1, crs="local")
f <- rast(ncol=50, nrow=50, xmin=0, xmax=1, ymin=0, ymax=1, crs="local")
p <- as.polygons(r)
set.seed(515)
xy <- lapply(1:ncell(r), \(i) spatSample(crop(f, p[i]), i, values=FALSE, xy=TRUE))
xy <- do.call(rbind, xy)
f <- rasterize(vect(xy), f)

There are 25 cells of f in each cell of r. In the first cell of r, only 1 cell of f has a value. In the second cell of r, two cells of f have a value. etc.

Aggregate and resample give the same result if the rasters are aligned

a <- aggregate(f, 10, sum, na.rm=TRUE)
#plot(a); text(a)

x <- resample(f, r, "sum")
# plot(x); text(x)

plot(as.lines(r),  axes=F)
plot(f, col="red", add=T, leg=F, axes=F)
text(x, col="blue")

The red blocks are the cells in f that have values (always 1). The lines are the cells of r and the text shows the values of r.

Now shift r half a cell such that the rasters are not aligned; and things do not look good because of the way that cells are treated that cross a boundary.

rr <- shift(r, .01, .01)
xx <- resample(f, rr, "sum")
plot(as.lines(rr), axes=FALSE)
plot(f, col="red", add=T, leg=F, axes=F)
text(xx , col="blue")

But with the "right" tweaking, it can be really terrible

rrr <- extend(r, c(1,1))
rrr <- shift(rrr, .025, .025)
xxx <- resample(f, rrr, "sum")
plot(as.lines(rrr), axes=F)
plot(f, col="red", add=T, axes=F, leg=F)
text(xxx, col="blue")

At this point, I am just trying to understand the problem. I do not know where this goes wrong (in terra or the GDALwarp method that is called in this case)

[rhijmans]

Continuing on the above, I get the same result with "stars", suggesting that the issue is with GDALwarp

library(stars)
src <- st_as_stars(f)
dst <- st_as_stars(rrr)
y <- st_warp(src, dst, use_gdal=T, method="sum", no_data_value=-99)
y <- rast(y)
plot(y); text(y)
lines(rrr)

[rhijmans]

The same happens when directly calling the gdalwarp tool (via sf). I created the input file with the code above and

writeRaster(f, "f.tif")

And you can download it herehttps://github.com/rspatial/files/blob/master/f.tif

library(sf)
#Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE

gdal_utils("warp", "f.tif", "warped.tif", options=c("-te", 0.025, 0.025, 1.025, 1.025, "-tr", 0.2, 0.2, "-r", "sum", "-overwrite"))
r <- rast("warped.tif")
r
#class       : SpatRaster 
#dimensions  : 5, 5, 1  (nrow, ncol, nlyr)
#resolution  : 0.2, 0.2  (x, y)
#extent      : 0.025, 1.025, 0.025, 1.025  (xmin, xmax, ymin, ymax)
#coord. ref. : Cartesian (Meter) 
#source      : warped.tif 
#name        : lyr.1 

plot(r, axes=F, legend=F)
lines(rast(r))
text(r)

[rhijmans]

The same problem with using "max", but with an older version of gdalwarp

$ gdalinfo --version
# GDAL 2.2.3, released 2017/11/20

And

$ gdalwarp f.tif -te 0.025 0.025 1.025 1.025 -tr .2 .2 -r max -overwrite out.tif
Creating output file that is 5P x 5L.
Processing input file f.tif.
Using internal nodata values (e.g. nan) for image f.tif.
Copying nodata values from source f.tif to destination out.tif.
0...10...20...30...40...50...60...70...80...90...100 - done.

The results are good. In this case, each cell has the value 1. The option "-r sum" is not available in this version of GDAL,

y = rast("out.tif")
plot(y, axes=F, col="light blue")
lines(rast(y))

But not with GDAL 3.4.3

gdal_utils("warp", "f.tif", "maxed.tif", options=c("-te", 0.025, 0.025, 1.025, 1.025, "-tr", 0.2, 0.2, "-r", "max", "-overwrite"))
x <- rast("maxed.tif")
plot(x, axes=F, col="light blue")
lines(rast(x))

And with the R script I get a good result too on the system with GDAL 2.2.3

gdal()
#[1] "2.2.3"
rrr <- extend(r, c(1,1))
rrr <- shift(rrr, .025, .025)
xxx <- resample(f, rrr, "max")
plot(rrr, legend=F)
plot(f, add=T, col="red", leg=F)
lines(rast(rrr))

[zoeschindler]

Thank you very much for looking into this so thoroughly!

[rhijmans]

This should now work with the development version of GDAL. The current version is 3.5.2; it was released in Sept 2022. It may take a while (several months at least) before the new version of GDAL is released and incorporated into the windows version of terra.

[rhijmans]

"fixes #847" should have been "fixes #874" --- it is not related to this issue.