Raster resampling can be performed using the resample function. The required parameters of this function are as follows:

The resample function currently provides two resampling methods. The method that we use determines the way in which a cell in the new raster gets its value, based on the old raster's values:

The nearest-neighbor interpolation retains the original values, only transferring them to the new grid according to proximity, while resampling methods that use weighted averages from several cells (such as bilinear interpolation) inevitably involve modification of the original values thus producing a smoother image. Therefore, in cases where the original values should be preserved, the nearest-neighbor resampling method is recommended. In categorical rasters with integers representing different cover types, for example, it would make no sense to calculate a weighted average of categories. For instance, in a raster with two categories, say 1 for forest and 2 for pasture, the average of the four cell values is frequently not going to be equal to either 1 or 2, making the resulting image meaningless. On the other hand, when visual appearance is the primary consideration, bilinear interpolation and similar weighted-average methods are recommended.

For our resampling examples, we will use two NDVI images: one from Landsat and one from MODIS. As the Landsat and MODIS images are defined in the same CRS and they overlap (the MODIS image extent encompasses the Landsat image extent)—although they have very different resolutions (the Landsat image resolution is 30 m while the MODIS resolution is 500 m)—we can resample each one with the other. In the following three examples, we will perform these steps:

For a Landsat NDVI image, we will use the ndvi_00 raster that we previously calculated based on an image obtained on October 4, 2000. As for MODIS, the raster r contains 280 NDVI images for the period between February 18, 2000 and April 6, 2012 (refer to Chapter 4, Working with Rasters). To make the comparison more interesting, we will select the date of acquisition closest to October 4, 2000. Using the dates$date vector of MODIS acquisition dates (as we saw in Chapter 3, Working with Tables), we can find out the index of the date closest to October 4, 2000 by performing the following steps:

These three steps can be accomplished with a single line of code, as follows:

> l_date = which.min(abs(dates$date - as.Date("2000-10-04")))

The result, assigned to l_date, is the index of the MODIS image closest in its acquisition date to October 4, 2000. It is equal to 15:

> l_date
[1] 15

Thus, for the examples, we will use r[[l_date]], the fifteenth layer of the multiband raster r (the corresponding date of acquisition is September 29, 2000, as dates$date[l_date] will show).

In our first example, we will resample the Landsat image ndvi_00 to the MODIS grid of the r[[l_date]] image using the nearest-neighbor method as follows:

> l_resample = resample(ndvi_00, r[[l_date]], method = "ngb")

Let's compare the original Landsat image ndvi_00 and the resampled image l_resample by plotting both as follows:

> plot(ndvi_00, main = "Original Landsat image")
> plot(l_resample, ext = extent(ndvi_00),
+ main = "Resampled to MODIS")

Note that when plotting l_resample, we supply the ext parameter with the extent of the original Landsat image ndvi_00 to zoom in on the area of overlap. The l_resample image is, in fact, larger than ndvi_00 (it has the same extent as r[[l_date]]) but all cells that do not overlap with ndvi_00 are assigned with NA.

The resulting two images are shown in the following screenshot:

On the left, we can see the detailed 30 meter resolution ndvi_00 image and on the right we see the rough 500 meter resolution image where the values of the nearest pixels of ndvi_00 were assigned to each of the pixels in r[[l_date]].

The opposite operation is to resample the MODIS NDVI values from r[[l_date]] to the Landsat grid of the ndvi_00, again using the nearest-neighbor method:

> r_resample = resample(r[[l_date]], ndvi_00, method = "ngb")

In this case, the result r_resample is smaller in its extent than the original MODIS image r[[l_date]]; it has the same extent as ndvi_00. To show them side-by-side more conveniently, we will first extend r_resample to match the extent of r[[l_date]]. This can be done using the extend function that extends the raster according to and an Extent object (or an object from which an Extent object can be derived, such as a RasterLayer). The extend function adds NA rows and columns as necessary to increase the raster extent; it is therefore, in a way, the opposite of trim:

> r_resample = extend(r_resample, r[[l_date]])

We will now plot the results of our second resampling example: the original MODIS image r[[l_date]] and the resampled image r_resample. While plotting r[[l_date]], we will add a rectangle surrounding the extent of ndvi_00 to aid in the comparison:

> plot(r[[l_date]], main = "Original MODIS image")
> plot(extent(ndvi_00), add = TRUE)
> plot(r_resample, main = "Resampled to Landsat")

The two resulting images are shown in the following screenshot:

As we can see, the resampled image is identical in its appearance to the respective portion of the original image, although it has a 30 meter resolution unlike the original image that has a 500 meter resolution. The reason is that each set of 30 m Landsat pixels that coincides with a single 500 meter pixel in the MODIS image gets the value of that single pixel. Thus, it appears as if we have 500 meter pixels even if the underlying resolution is much higher.

What if we want to make the image look smoother? In such a case, we can use the bilinear interpolation method instead of the nearest-neighbor method. Let's perform the resampling operation from the last example twice, using the nearest-neighbor method (the same way we just did, assigning the result to r_resample_ngb this time) and the bilinear interpolation method (assigning the result to r_resample_bil). The results will then be combined in a two-band raster with layer names specifying resampling method:

> r_resample_ngb = resample(r[[l_date]], ndvi_00, method = "ngb")
> r_resample_bil = resample(r[[l_date]], ndvi_00,
+ method = "bilinear")
> resample_results = stack(r_resample_ngb, r_resample_bil)
> names(resample_results) = c("Nearest neighbor",
+ "Bilinear interpolation")

We will plot the results using the levelplot function of the rasterVis package. Contours are enabled with contour=TRUE to highlight the differences between methods:

> library(rasterVis)
> levelplot(resample_results,
+ par.settings = RdBuTheme,
+ contour = TRUE)

The following screenshot shows the graphical output:

As we saw in the previous example with nearest-neighbor resampling (the left image), sets of adjacent cells that overlap with a single 500 meter pixel get the same value, and thus the contours follow the coarse-grained pattern of MODIS pixels. With bilinear interpolation (the right image), on the other hand, each 30 m cell gets the average of four MODIS pixels closest to it, weighted according to their respective distances; thus, most (if not all) cell values are unique and the image is much smoother.

As mentioned earlier, raster reprojection—unlike vector layers reprojection—involves resampling. Unlike vector layers, where all points are independent and thus their coordinates can be individually transformed to a new CRS (refer to the previous chapter), raster pixel coordinates that are transformed to a different CRS will not form a legitimate grid on that CRS (that is, a grid with equal distances between all adjacent pixels and parallel to the CRS x and y axes); rasters therefore require resampling. Thus, raster reprojection consists of two steps: reprojection of the raster pixel coordinates to the new CRS (analogous to reprojecting a vector layer) and resampling of the pixel values to a grid defined in the new CRS.

There are two main ways of defining the new grid in raster reprojection: we can either provide only the new CRS definition and let the grid be generated automatically, or we can provide a specific grid of our own. In the latter case, the process is very similar to resampling; the only change is that the new grid has a different CRS from that of the original raster.

In R, the projectRaster function provides raster reprojection functionality and supports both previously mentioned ways of defining a new grid. In the first case, when we want the grid to be automatically generated, we need to supply:

In the second case, when we want to define the grid ourselves, we need to supply:

In both scenarios, we also need to specify the resampling method, either "ngb" or "bilinear" (the default).

As an example, we will reproject dem from its geographic CRS to the UTM Zone 36N projection. The corresponding PROJ.4 string, as we have already done in the previous chapters, will be obtained from another raster (r). Prior to reprojection, let's print out the current properties of dem in order to compare them later with the properties of the reprojection result:

> dem
class       : RasterLayer
dimensions  : 390, 384, 149760  (nrow, ncol, ncell)
resolution  : 0.0008333333, 0.0008333333  (x, y)
extent      : 34.83083, 35.15083, 32.63583, 32.96083  (xmin, xmax$
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,$
data source : in memory
names       : layer
values      : -14, 541  (min, max)

We see that, indeed, dem is defined in a geographic CRS, and has a resolution of 0.0008333333 degrees, longitude and latitude.

In case we do not have any special requirements for the projected dem to align with an existing layer, we can follow the first reprojection approach, supplying only the CRS and letting the grid be generated automatically. We can also supply the optional res argument, primarily to make sure the x and y axes resolution will be equal (90 m, for example). You can try and execute the following expression without specifying res, to verify that the default grid generated in this case will have a resolution of 78 m and 92.4 m on the x and y axes respectively:

> dem = projectRaster(from = dem,
+ crs = proj4string(r),
+ method = "ngb",
+ res = 90)

Printing the properties of the reprojected dem reveals several differences:

> dem
class       : RasterLayer
dimensions  : 417, 351, 146367  (nrow, ncol, ncell)
resolution  : 90, 90  (x, y)
extent      : 670666.3, 702256.3, 3611918, 3649448  (xmin, xmax, $
coord. ref. : +proj=utm +zone=36 +datum=WGS84 +units=m +no_defs +$
data source : in memory
names       : layer
values      : -14, 541  (min, max)

We can see that the CRS definition has changed to UTM Zone 36N. Also, the extent, resolution, and dimensions are now different and characterize the new grid that has been generated.

Plotting the new raster with plot(dem) produces the following graphical output:

UTM Zone 36N coordinates are now shown along the axes. In addition, you may have noticed that the projected dem image is not exactly parallel to the axes, as the UTM grid was not strictly parallel to the geographic one. Since a raster must be rectangular (refer to the preceding screenshot), empty areas around image margins, consisting of NA values, have been generated.