Creating a new raster image from averages

Introduction

In this article, we’ll explore how to create a new raster image that represents the average of a certain number of rasters in a GIS (Geographic Information System). This process is commonly used in remote sensing and geospatial analysis, where large datasets need to be processed efficiently. We’ll walk through the steps involved in creating such an image using RasterStack, a package for working with raster data in R.

Background

Raster images are composed of pixel values that represent various characteristics of the data, such as elevation, vegetation index, or temperature. When working with rasters, it’s often necessary to perform operations that involve averaging or aggregating values across a spatial area. In this article, we’ll focus on creating a new raster image that represents the average of two adjacent rasters.

Setting up the workspace

Before diving into the code, let’s set up our R workspace. We assume that all the rasters are stored in a single folder named rasdir. To start, we need to:

Set the rasdir environment variable
List the files (rasters) in the rasdir folder using list.files
Create a vector of raster paths (raspaths) with full names

# Load necessary libraries
library(raster)

# Set the rasdir directory
rasdir <- "myrasters/"

# Get the list of files (rasters) in the rasdir directory
raspaths <- list.files(path = rasdir, full.names = TRUE)

Stacking rasters

Once we have our raster paths, we can stack them using stack(). This creates a single raster object (rascube) that contains all the individual rasters.

# Stack the rasters
rascube <- stack(raspaths)

Creating the averaging function

The next step is to create a function that calculates the moving average of adjacent rasters. We’ll use a matrix (mat) to determine how the rasters are aggregated. The for loop will iterate over each row in the matrix, computing the moving average and writing it to a new raster file.

# Create the averaging function
rasfun <- function(x = rascube, win = 2, outdir = getwd()) {
  
  # Sanity check: Ensure the number of rasters is divisible by the window size
  if (!(length(raspaths) / win) %% 1 == 0) {
    stop("Number of rasters must be divisible by window size")
  }
  
  # Create an index matrix that determines how rasters are aggregated
  mat <- matrix(data = 1:length(raspaths), ncol = win, byrow = TRUE)
  
  # Loop over the index matrix to calculate the moving average
  for (i in 1:nrow(mat)) {
    
    # Compute the ith moving mean
    res_i <- sum(x[[mat[i,1]:mat[i,win]]], na.rm = TRUE) / win
    
    # Write the output to a file with the respective input raster name
    writeRaster(x = res_i, filename = paste(mat[i,1:win], collapse = "", sep = "_"),
                format = "GTiff", overwrite = TRUE)
  }
}

Running the averaging function

Finally, we can run the rasfun function on our stacked rasters with a window size of 2 to create new raster images that represent the average of adjacent rasters.

# Run the averaging function
rasfun(x = rascube, win = 2, outdir = "moving_mean_rasters/")

Conclusion

In this article, we explored how to create a new raster image that represents the average of two adjacent rasters. We used RasterStack to stack our individual rasters and created a function (rasfun) that calculates the moving average. By following these steps, you can efficiently process large datasets in your GIS environment.

Notes

The number of rasters must be divisible by the window size to avoid errors.
You can modify the rasfun function to suit your specific needs and processing requirements.
Always ensure that the output raster files are properly labeled with their respective input raster names.

Last modified on 2025-04-17