I would like to resample a high resolution raster to a coarser resolution, but in such a way that the maximum values of cells are retained for the coarser grid cells.
As there is no fun argument in the resample function in R's raster package, I have put together a simply custom function:
resampleCustom <- function(r1, r2) {
resRatio <- as.integer(res(r2) / res(r1))
ret <- aggregate(r1, fact = resRatio, fun = max)
if (!compareRaster(ret, r2, stopiffalse = FALSE)) {
ret <- resample(ret, r2, method = 'bilinear')
}
return(ret)
}
Basically, I use aggregate, where I can provide a custom function, to get close to the target raster, and then I use resample to apply some final adjustments.
I applied this to a raster that represents the projected distribution of a species of fish (where cell values represent suitability scores ranging from 0 to 1), and the odd thing is that the resulting raster has values that are greater than the max values in the original rasters.
The two rasters can be downloaded here and here.
library(raster)
# read in species raster and template
sp <- raster('Abalistes_filamentosus.tif')
template <- raster('rasterTemplate.tif')
> sp
class : RasterLayer
dimensions : 360, 720, 259200 (nrow, ncol, ncell)
resolution : 48243.14, 40790.17 (x, y)
extent : -17367530, 17367530, -7342230, 7342230 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs
data source : /Users/pascaltitle/Dropbox/Abalistes_filamentosus.tif
names : Abalistes_filamentosus
values : -5.684342e-14, 1 (min, max)
> template
class : RasterLayer
dimensions : 49, 116, 5684 (nrow, ncol, ncell)
resolution : 3e+05, 3e+05 (x, y)
extent : -17367530, 17432470, -7357770, 7342230 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs
data source : /Users/pascaltitle/Dropbox/rasterTemplate.tif
names : rasterTemplate
values : 1, 1 (min, max)
> resampleCustom(sp, template)
class : RasterLayer
dimensions : 49, 116, 5684 (nrow, ncol, ncell)
resolution : 3e+05, 3e+05 (x, y)
extent : -17367530, 17432470, -7357770, 7342230 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs
data source : in memory
names : Abalistes_filamentosus
values : -0.2061382, 1.206138 (min, max)
The max value is 1.2, but how can this be when the bilinear method should essentially be taking averages of cell values? I would expect all values of the resulting raster to be within the bounds of the original raster values.
The extreme values are for cells at the edge of the raster, where values are extrapolated, as there are no neighbors at one side. This shows where these values are:
x <- resampleCustom(sp, template)
a <- xyFromCell(x, which.max(x))
b <- xyFromCell(x, which.min(x))
plot(x)
points(a)
points(b)
Or
plot(Which(x < 0))
plot(Which(round(x, 15) > 0))
To remove these extreme values, you can use raster::clamp.
xc <- clamp(x, 0, 1)
By the way, what you do, first aggregate then resampling, is also what is done within raster::resample.
The fundamental problem is that your high-res raster data do not line up with the low resolution aggregation you are seeking. That suggests a mistake earlier on in your work flow. The best way to avoid this problem is probably to make the habitat suitability predictions with predictor raster data that are aligned with the high resolution raster. You perhaps did not consider that when you projected the predictor variables to +proj=cea?