It’s generally best to avoid resampling raster datasets.

This is because information can be lost during the mathematical reshaping process that happens when we convert grids from one projection system to another. Below, we will see that reprojecting a raster dataset can dramatically affect the data that it contains.

The setup

Load the relevant packages. Also, set a randomization seed so that every time you run this code from the top, it will reproduce the same results.

# Load packages, set seed
library(raster)
library(USAboundaries)
library(sf)
library(ggplot2)

# Set seed for replicability
set.seed(-356) # The start of the conquests of Alexander the Great

Make up some data, and create a raster object. Assign the raster a few relevant spatial properties so that we can contextualize what’s going on.

# A made-up raster dataset located in northeast Arizona
madeUpData <- matrix(runif(10000,0,100),
                     nrow=100,
                     ncol=100)

# Set up extent information of bounding box
xmn = 300000
xmx = 400000
ymn = 4000000
ymx = 4100000

# Start out with a common coordinate reference system: 
# The Universal Transverse Mercator (units are in meters)
# Zone 13 runs up through Arizona and New Mexico
AZ_UTM = "+proj=utm +zone=13 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"

# Compile the raster using the data, projection information, and extent information
baseRaster <- raster(madeUpData,
                     crs = AZ_UTM,
                     xmn = xmn, xmx = xmx,
                     ymn = ymn, ymx = ymx)
plot(baseRaster,
     main = "Original raster")

Let’s get ourselves oriented by plotting the extent of the data over the 4-corners region in the US Southwest.

# Access the 4 corners states using the USAboundaries library
myStates <- us_states(resolution = "low", 
                      states = c("Arizona",
                                 "New Mexico",
                                 "Utah",
                                 "Colorado"))

# Spatially transform the states lines to the baseRaster crs
myStates <- st_transform(myStates,crs(baseRaster))

# Look at where the made up data lives
ggplot() +
  geom_sf(data = myStates) +
  geom_rect(aes(xmin=extent(baseRaster)@xmin,
                xmax=extent(baseRaster)@xmax,
                ymin=extent(baseRaster)@ymin,
                ymax=extent(baseRaster)@ymax))

Reproject ad nauseum

Let’s generate a list of CRS systems to work with. These 9 epsg codes all specialize in Arizona, so are (approximately) comparable reference systems to the starting one. I searched “Arizona” in [epsg.io] to find them. Next, use a for loop to iterate the reprojection process using each of the CRS systems. Finally, project back to the starting point using the baseRaster object as a template.

# List of potential CRS systems to project to
mySystems <- c("+init=epsg:2222",
               "+init=epsg:2223",
               "+init=epsg:2224",
               "+init=epsg:2761",
               "+init=epsg:2762",
               "+init=epsg:2763",
               "+init=epsg:2867",
               "+init=epsg:2868",
               "+init=epsg:2869")

# Make a new raster based on the one we made above
newRaster <- baseRaster

# Use a `for` loop to iterate over projections
for(i in 1:length(mySystems)){ 
  # For each of the projections listed above...
  myCRS <- mySystems[i]
  # ... project the last raster to the new projection
  newRaster <- projectRaster(newRaster,
                             crs = crs(myCRS),
                             method = "bilinear")
}

# Convert newRaster back to baseRaster projection
newRaster <- projectRaster(newRaster,
                           crs = crs(baseRaster),
                           method = "bilinear")

Compare inputs and outputs

First use the all.equal() function to get a boolean value and any relevant warnings indicating whether the starting raster and reprojection result are equivalent objects.

# Check if all cells are equal
all.equal(newRaster, baseRaster) # Nope!
## Warning in compareRaster(target, current, ..., values = values, stopiffalse =
## stopiffalse, : different extent
## Warning in compareRaster(target, current, ..., values = values, stopiffalse =
## stopiffalse, : different number or columns
## Warning in compareRaster(target, current, ..., values = values, stopiffalse =
## stopiffalse, : different number or rows
## Warning in compareRaster(target, current, ..., values = values, stopiffalse =
## stopiffalse, : not all objects have the same values
## [1] FALSE

Clearly there are some differences: everything from extent to cell values has changed. Plot the starting raster next to the result of several reprojections. Are they the same image? What differences arose?

# Plot the two raster sets side-by-side:
par(mfrow=c(1,2))
plot(baseRaster,
     xlim = c(xmn,xmx),
     ylim = c(ymn,ymx),
     main = "Original raster")
plot(newRaster,
     xlim = c(xmn,xmx),
     ylim = c(ymn,ymx),
     main = "New raster")

In fact, if we had not included the xlim and ylim arguments in the above calls to plot(), we would see that the newRaster has a much larger extent than the original baseRaster. Compare histograms to check out differences in the distribution of values.

# Plot histograms
par(mfrow=c(1,2))
hist(baseRaster, 
     breaks = 2500, 
     main = "Base raster")
hist(newRaster, 
     breaks = 2500, 
     main = "New raster")

What happened?

Hint. Think about what bilinear interpolation is doing.

Can you imagine a real-life scenario where this would affect your work?