TL;DR: I have to downsample, project and align a raster to perfectly fit another one. Is it better to a) use
aggregate followed by
projectRaster or b) first use
projectRaster for reprojection only followed by
resample for aggregation and alignment?
Similar questions have been asked before, but as far as I can tell, none of them covers all three aspects of downsampling, projecting and aligning, nor do they consider all potentially reasonable options.
My source file is SEDAC's Human Influence Index v2, more specifically the global data "in geographic coordinate system at 30 arc-second grid cell size" (SEDAC download page).
> library(raster) > hii <- raster("./hii_v2geo/") > crs(hii) CRS arguments: +proj=longlat +ellps=clrk66 +no_defs > res(hii)  0.00833333 0.00833333 > extent(hii) class : Extent xmin : -180 xmax : 179.9999 ymin : -65.23565 ymax : 89.89762
My target raster is in ETRS 89 LCC, at a coarser resolution and cropped to the German State of Bavaria (I know there's a Europe-only version of the HII, I've got my reasons to use the global file). I've put a reference grid filled with random values up for download here.
ref <- raster("reference.tif") > crs(ref) CRS arguments: +proj=lcc +lat_1=35 +lat_2=65 +lat_0=52 +lon_0=10 +x_0=4000000 +y_0=2800000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs > res(ref)  5000 5000 > extent(ref) class : Extent xmin : 3928000 xmax : 4273000 ymin : 2290000 ymax : 2645000 >
The task is now to perfectly align the hii raster with the reference (so that CRS, resolution and extent are identical) for use as an explanatory variable in a species distribution model.
The first step is obvious, I crop the hii raster to the extent of my study area, leaving some additional data at the edges:
hii <- crop(hii, extent(8.4, 14.4, 46.7, 51.1))
After that, there are three things to do:
- aggregate to the coarser resolution
- project to ETRS 89 LCC
- make sure raster cells are aligned
As far as I'm aware, there are three tools available for this in the raster package:
aggregate()can aggregate raster cells by a given aggregation factor (integer). It can not align cells and it can not change the projection.
aggregate()(source code) to first aggregate to "approximately equal resolution as the target", then uses bilinear interpolation (or nearest neighbor) to align the cells. It can not change the projection.
projectRaster()will do all three tasks if you let it, but can also do a simple projection without touching cell size or alignment
The question is now: which tools to use in which order, and why? There are several options available:
Option 1 is the quick and dirty way of just using projectRaster for the whole thing:
hii.new.1 <- projectRaster(from=hii, to=ref, method="bilinear")
However, this is discouraged in the command's documentation "if the resolution of the output is much larger than that of the input" (which i guess is the case here).
Option 2 is what is being suggested by
?projectRaster: first use
aggregate() to get close to the target resolution and then use
projectRaster() for projection and alignment of cells:
hii.new.2 <- aggregate(x = hii, fact=8, fun=mean) # value for "fact" found via trial and error: slightly finer than target resolution hii.new.2 <- projectRaster(from=hii.new.2, to=ref, method="bilinear")
Option 3 is quite the opposite – project first without any aggregation or resampling to retain as much information as possible for as long as possible, followed by
hii.new.3 <- projectRaster(from=hii, crs=crs(ref)) hii.new.3 <- resample(x=hii.new.3, y=ref, method="bilinear")
Let's have a look at the results from the different methods:
- Option 1 looks quite noisy
- Option 2 by comparison lacks detail
- Option 3 looks like a good compromise to me.
So Option 1 apparently is not really an option, but I wonder why Option 2 is suggested by the projectRaster() documentation and not Option 3. In digital photography, where warping/stretching a picture in order to perform "perspective correction" poses some similar problems, people tend to upsample (disaggregate) first, them transform and then downsample (aggregate) back to the original resolution or coarser in order to preserve as much detail as possible.
It was pointed out below that averaging and interpolating may not be the best choice for this kind of data, so I ran a second "B" set of tests using
fun='modal' instead. The results are interesting:
- Option B1 is very similar to Option 1
- Option B2 looks really bad
- Option B3 is very similar to Option B1 now.
Meanwhile, I've been trying to find some more objective measures. First of all, plotting the frequencies of values as histograms reveals interesting patterns if the bins are small enough:
The "bad" Option 1 comes closest to to the original distribution – confirmed by a quick Kruskal-Wallace test. Then, I generated 1k randomly placed points, extracted values from all rasters, and ran a simple correlation analysis (repeated a couple of times with different points and basically the same results):
+----------------+----------------+-------------+------+ | transformation | rho (Spearman) | R (Pearson) | R² | +================+================+=============+======+ | original | 1.00 | 1.00 | 1.00 | | option 1 | 0.76 | 0.79 | 0.62 | | option 2 | 0.77 | 0.79 | 0.62 | | option 3 | 0.81 | 0.83 | 0.69 | | option B1 | 0.73 | 0.77 | 0.59 | | option B2 | 0.62 | 0.64 | 0.41 | | option B3 | 0.68 | 0.71 | 0.50 | +----------------+----------------+-------------+------+
This one quite clearly suggests that my old favorite Option 3 may actually represent the original data best.