How to change the resolution of a raster layer in R

Question

I have several high resolution raster layers in R that I am working with. The level of detail is excessive for some of the analyses I am running, so I would like to speed things up by reducing the resolution.

The coordinate system is UTM so the units are meters. The resolution says it is 30, 30 (x, y). So it seems that the resolution here is 30m.

Could someone please advise me on how to change the resolution to 120m instead? I have read the help for the resample() and projectRaster() functions but they seem to require a template raster with the desired resolution, which I don't have.

Here is an example of one of my raster layers:

alt.utm
class : RasterLayer
dimensions : 4572, 2495, 11407140 (nrow, ncol, ncell)
resolution : 30, 30 (x, y)
extent : 421661, 496511, 4402939, 4540099 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=13 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
data source : in memory
names : layer
values : 1485.127, 4275.202 (min, max)

Terce

Answer 1

You can use aggregate or disaggregate.

library(raster)

#get some sample data
data(meuse.grid)
gridded(meuse.grid) <- ~x+y
meuse.raster <- raster(meuse.grid)
res(meuse.raster)
#[1] 40 40

#aggregate from 40x40 resolution to 120x120 (factor = 3)
meuse.raster.aggregate <- aggregate(meuse.raster, fact=3)
res(meuse.raster.aggregate)
#[1] 120 120

#disaggregate from 40x40 resolution to 10x10 (factor = 4)
meuse.raster.disaggregate <- disaggregate(meuse.raster, fact=4)
res(meuse.raster.disaggregate)
#[1] 10 10

Noeminoesis

Thank you for this simple solution! This worked perfectly. – Terce

doesn't work if x and y resolution different (e.g. rectangular pixels) – Kalil

Service message: Even though using raster functions should be fine most of the time, keep in mind, future reader, that rgdal, rgeos, maptools, and, by extension, many of their dependencies, are being retired by the end of 2023. "Plan transition to sf/stars/terra functions [...] at your earliest convenience." Read more about it here. – Corinacorine

Answer 2

These days we can use terra, the replacement for the raster package

Example data

library(terra)
f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)

Aggregate raster cells. To go from 30 to 120m is a factor of 4

a1 <- aggregate(r, 4)

You can use a different factor for the rows and columns, and also different aggregation functions (the default is "mean")

a2 <- aggregate(r, c(2,3), fun=sum)

You can also go the other way and disaggregate:

d <- disagg(a2, 2)

Aggregate can only combine entire cells. But to combine raster data from different sources, you may need to match the geometry of a raster that does not align. In that case you can use resample

Example of a non-aligned SpatRaster:

x <- rast(r)
res(x) <- 0.01

Solution

x <- resample(r, x)

It could also be that you want to transform raster data to a geometry with another coordinate reference system ("map projection"). You can do that like this:

u1 <- project(r, "+proj=utm +zone=32")

But, unlike for vector data, the transformation of raster data is not well-defined. Therefore, the preferred approach is to provide a template to which you want the output to align with. This would normally be a SpatRaster that you already have from another data source. But here I create one for the example:

temp <- rast(xmin=264000, xmax=324000, ymin=5479000, ymax=5565000, res=100, crs="+proj=utm +zone=32")

And now use it

u2 <- project(r, temp)

Further reading

Anaphrodisiac

Answer 3

I've tried 3 different options for upscaling a DEM file. First I used gdal_translate like this:

from CMD:
gdal_translate -tr 0.1 0.1 "C:\dem.tif" "C:\dem_0.1.tif"

In R:

res(raster('C:\dem_0.1.tif')
[1] 0.1 0.1

Then I tried both aggregate and resample commands from raster package in R:

#
r <- raster("dem.tif")
> res(r)
[1] 0.0002777778 0.0002777778 # original dem resolution
#
r_up_0.1 <- aggregate(r, fact = 0.1/res(r)) # aggregate output
> res(r_up_0.1)
[1] 0.1 0.1
#
s <- raster(nrow = 10, ncol = 10)
extent(s) <- extent(r)
s <- resample(r, s, method = 'bilinear') # resample output
> res(s)
[1] 0.1000278 0.1000278

These are the outputs: The 3 outputs were upscaled to a res = 0.1x0.1 and there are some differences between them but I intend to use the gdal output.

Hope this helps.

Wormeaten

how did the pre and post totals compare for all 3 methods? i do these all the time, and the cell values (and raster total) can vary quite a bit – Steinke

Answer 4

Here is an example of how to do this. (link to original)

  #########################################################################
# 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
}

Songer

Thanks scribbles - this is useful to know for manipulating resolution in either direction. – Terce

Answer 5

This is a late answer, but creating a template raster and using projectRaster is very straightforward. You can set the template raster resolution to any value you would like, including two values if you want rectangular rather than square cells:

# Create the template raster, with 120m cells
TEMPLATE.RASTER <- raster(extent(ORIGINAL.RASTER), resolution = 120, 
                          crs = st_crs(ORIGINAL.RASTER)$proj4string)

# Project the original raster to the new resolution
PROJECTED.RASTER <- projectRaster(from = ORIGINAL.RASTER, to = TEMPLATE.RASTER)

I like to create a template, as it copies over the target CRS as well. However, you can actually skip creating a template raster and set the resolution and CRS directly in projectRaster. Do check the help files, as there are options for calculation method etc which might be important to consider.

Uboat

Answer 6

I recently had a requirement to reduce the resolution of a ggmap object. This involves extracting and transforming the ggmap raster (using Robin Lovelace's ggmap_rast()), aggregating the raster as discussed in this thread, and then replacing the ggmap high resolution raster with the lower resolution raster below. Hope this is useful:

original_map <- get_map("New York", scale = 1) #download a map at lowest resolution
#ggmap_rast function courtesy of Robin Lovelace
original_map.rast1 <- ggmap_rast(original_map) #extract RasterStack from ggmap object

original_map.rast2 <- aggregate(original_map.rast, 2) #compress raster

rast2_length <- sqrt(length(original_map.rast2)/3)  ## find the number of cells in compressed raster

map <- ggmap::ggmap(original_map) 
# from https://mcmap.net/q/427367/-plot-ggmap-image-over-raster
r <- map$layers[[2]]$geom_params$raster #pull hex raster pixel values from ggmap object
#x <- r[,] #assign values to a variable (probably unnecessary)

xx <- r[1:rast2_length,1:rast2_length] #filler values for raster to be created

rgb_map <- original_map.rast2

for(i in 1:rast2_length){
  for(j in 1:rast2_length){
    k=i*rast2_length+j
    #many rgv values are non-integers; rgb2hex requires integer
    red = as.integer(round(rgb_map$layer.1[k],0))
    green = as.integer(round(rgb_map$layer.2[k],0))
    blue = as.integer(round(rgb_map$layer.3[k],0))
    #rgb2hex from ggtern package
    xx[i,j] <- rgb2hex(red,green,blue)
    #rgb2hex is slow; export, calculate in excel, import is faster, sadly
    #xx[i,j] <- as.character(rgb_from_excel$V1[k])
  }
}

map$layers[[2]]$geom_params$raster <- xx
map

Mononucleosis