Resample Raster Images

This example demonstrates two R programming techniques for resampling (changing the cell size or spatial resolution) of a raster satellite image. The first technique uses data objects and methods from the sp package. The second technique uses the recently released (Fall, 2009) raster package.

We resample the same satellite image, and compare the results of all sampling methods by computing elementary statistics.

Download data and R Code for this example

When constructing multi-layer raster geospatial data bases from several data sources (e.g., satellite images, digital terrain models), analysts must often modify one or more of the layers so that each data set layer has the same spatial resolution. The technique used to change the spatial resolution of a raster image is known as image resampling

Two R packages enable image resampling: Within the sp package, the Spatial Pixel / Spatial Grid classes to translate an image raster with one spatial resolution into a new grid with a different resolution. The raster package s rasterStack and rasterBrick classes, which can contain multiple raster images, enable several images to be resampled with a single command.

Here are two examples.

Example 1: Using sp package

Satellite Image: MODIS 8-Day (Average) Land Surface Temperature Image, 500 meter resolution

Workflow:

1) Read input files into R Spatial data objects.

2) Display the satellite image, creating custom gray/color scales for the image .

3) If necessary, transform one or more of the input files into a new spatial frame of reference (i.e., map projection) to make it compatible with the other file.

4) Perform Point Overlay operation: Assign to each point location the spatially corresponding satellite image pixel value.

5) Write the point locations, along with assigned MODIS surface temperature values, to a new ESRI Shape File, as well as a separate Comma Separated Value file.

Discussion:

R Script:

 #########################################################################
# SubsampleImageSpGdal.r
#  
# This function reads a raster image file, and creates an output
# file with coarser spatial resolution. Technique used: Subsampling.
# Create a sampling grid at a different (usually coarser) resolution
# than the input image, 'overlay' the sampling grid on the input
# image, extract a point from the image at each grid location,
# and then write the sampling grid to an output file in .img format.
# 
# Author: Rick Reeves
# Date created: 29-Apr-2009    
# Date modified: 30-Sept-2010                                                            
# NCEAS
#
#########################################################################
#
SubsampleImageSpGdal <- function()
{
   library(rgdal) # loads sp package
   library(Hmisc)
 
# SampleFactor: the reduction in image size and resolution: 1 / SampleFactor
# e.g., Sample Factor  = 2 creates an image with 1/2 the resolution 
# and 1/4 the number of pixels (x and y dimension subsampling)   
 
   SampleFactor <- 2

# Step 1) Read input image into a SpatialGridDataFrame object, then display.
 
   setwd("./Img1")
   InputImage <- readGDAL("B30SubNP.tif")      

# Create a sampling grid, using the input image's spatial parameters  
   
   SamplingGridTopology <- GridTopology(InputImage@bbox[,1],
                                        InputImage@grid@cellsize * SampleFactor,
                                        InputImage@grid [at] cells [dot] dim / SampleFactor)

# Create the necessary data structures
                                       
   SamplingGrid <- SpatialGrid(SamplingGridTopology)  
   OutSpatialGridSGDF <- as(SamplingGrid,"SpatialGridDataFrame")
   SamplingPoints <- as(SamplingGrid,"SpatialPoints")   

# Create the subsampled image in a SpatialGridDataFrame object.'
  
   SamplingPoints <- as(as(SamplingGrid,"SpatialPixels"),"SpatialPoints")
   SamplingResultsSPDF <- overlay(InputImage,SamplingPoints)    
   OutSpatialGridSGDF@data <- SamplingResultsSPDF@data 
   OutSpatialGridSGDF@proj4string <- CRS(proj4string(InputImage))   
   gridded(OutSpatialGridSGDF) <- TRUE  # necc for writeGDAL to write image.
 
# use Hmisc library describe() function to compare the incoming and resampled images
 
   Incoming = describe(InputImage@data)
   print(Incoming)
   Outgoing = describe(SamplingResultsSPDF@data)  
   print(Outgoing)
   message("Summary of incoming and outgoing images (hit key to continue):")   
   browser() 
#                                                                                                                             
   writeGDAL(OutSpatialGridSGDF,"B30SubResampled.tif")
   
# Finished
    
}

Example 2: Using raster package.

Satellite Image: MODIS 8-Day (Average) Land Surface Temperature Image, 500 meter resolution

Workflow:

1) Read input files into R Spatial data objects.

2) Display the satellite image, creating custom gray/color scales for the image .

3) If necessary, transform one or more of the input files into a new spatial frame of reference (i.e., map projection) to make it compatible with the other file.

4) Perform Point Overlay operation: Assign to each point location the spatially corresponding satellite image pixel value.

5) Write the point locations, along with assigned MODIS surface temperature values, to a new ESRI Shape File, as well as a separate Comma Separated Value file.

Discussion:

R Script:

  #########################################################################
# SubsampleImageRaster.r
#  
# This function demonstrates resampling of raster images to a new
# spatial resolution using the R raster package.
# 
# Author: Rick Reeves
# Date created: 6-October- 2010
# Date modified:                                                             
# NCEAS
#
#########################################################################
#
SubsampleImageRaster <- function()
{
   library(raster)

   sOutFile <- ""
   resampleFactor <- 4  # For test, subsample incoming image by factor of 10
                 
# Read the mosaic components, stored in a subfolder, into a raster object list.
# Within the same loop, obtain the (geographic) extent of each component.
# Note: these images do not have same spatial extent, so they cant be stored
# in a rasterStack. Instead, use a list of rasterLayers.
                                         
   setwd("./ForUseCase")
   inFiles <- list.files(pattern="*.tif")
   nFiles <-  length(inFiles)
   inputRaster <- raster()
   CoarseResampRaster <- raster()   
   FineResampRaster <- raster()   
   for (iCtr in 1 : nFiles)
   {
      message(sprintf("resampling file: %s",inFiles[iCtr]))
      inputRaster <- raster(inFiles[iCtr])
     
# The aggregate() / disaggregate methods resample rasters to COARSER (bigger cells)
# and FINER (smaller cells) resolutions, respectively
      
      CoarseResampRaster <- aggregate(inputRaster,fact=resampleFactor,fun=mean) 
      sOutFile <- sprintf("CoarseSubsamp%s",inFiles[iCtr]) 
      writeRaster(CoarseResampRaster,filename=sOutFile,format="GTiff",datatype="INT1U",overwrite=TRUE)      
      FineResampRaster <- disaggregate(inputRaster,fact=resampleFactor,fun=mean) 
      sOutFile <- sprintf("FineSubsamp%s",inFiles[iCtr]) 
      writeRaster(FineResampRaster,filename=sOutFile,format="GTiff",datatype="INT1U",overwrite=TRUE)            
   }   
   
   message("resample demo")
   browser()

# second method: use the resample() method from raster package
   
# Simple example:
# This code produces a resampled raster, 's',
# with correct resampled values. e.g.; 
# s[] prints a vector of resampled cell values.

   r <- raster(nrow=3, ncol=3)
   r[] <- 1:ncell(r)
   s <- raster(nrow=10, ncol=10)
   s <- resample(r, s, method='bilinear')
     
# Useful example:
# Resample a satellite image, stored in a GeoTiff file
# into a NEW raster with 2x spatial resolution in 
# both dimensions (four times the number of cells)
# Here is the technique: 
#  1) Create a new raster object with the correct 'resampled' number of cells.
#  2) Set the extent (geographic 'bounding box') of the new raster 
#     to the extent of the original raster
#  3) Generate the resampled raster.
 
   resampleFactor <- .5  # reduce the cell size by 50% and double the number of rows and columns.      
   inputRaster <- raster("TmB50MosaicImg1.tif")      
   inCols <- ncol(inputRaster)
   inRows <- nrow(inputRaster)
   resampledRaster <- raster(ncol=(inCols / resampleFactor), nrow=(inRows / resampleFactor))
   extent(resampledRaster) <- extent(inputRaster)

# The resample method will write the resampled raster image to a NEW disk file..

   resampledRaster <- resample(inputRaster,resampledRaster,datatype="INT1U",method='bilinear',filename="testOutResamp.tif",overwrite=TRUE)

# Or, use writeRaster method to create the output file.

   writeRaster(resampledRaster,filename="ResampleProduct.tif",format="GTiff",datatype="INT1U",overwrite=TRUE)   


# end
}

Learning More:

The MODIS Imagery Web Site - NASA Goddard Space Flight Center

The GIS Primer: The Nature of Spatial Information

scicomp: