Crop (spatially subset) rasters to an area

To subset the extent of a raster to a desired area, we use the 'crop' function in the 'raster' package. This is the most commonly used function to modify the extent of a raster object. The raster function takes two main arguments:

(1) a Raster* object (or an SpatialPolygons*, or SpatialLines*, or SpatialPoints* object)

(2) an Extent object or any object from which an Extent object can be extracted (i.e. RasterLayer, RasterStack, RasterBrick and objects of the Spatial* classes). The extent object will act as a cookie cutter.

For an adequate behaviour, the raster and extent need to be in the same Coordinate Reference System (CRS).

EXAMPLES

Below examples on how to crop raster objects is provided. These examples were taken from the “Effects of Cyclone Yasi on Green Cover at Mission Beach” tutorial. Putting the code snippets in context, by looking at a broader section of the R script, could be benefitial. Boxes with grey background contain code snippets, and boxes with white background containt code (text) outputs.

In the example, a SpatialPolygons object for the area of interest (i.e. the Study Area) had been previously created to prepare a map showing the Region & Study Area. The raster is in EPSG:3577(=‘GDA94’/‘Australian Albers’), an Australian-specific CRS (further details for the EPSG:3577 CRS in this link). The polygon for the Study Area, on the other hand, is in EPSG:4326 CRS, a general world spatial reference (more details for this CRS in this link. Therefore, we need to transform the CRS for the Study Area spatial polygon form EPSG:4326 to EPSG:3577. We could then directly use this re-projected polygon as an argument in the ‘crop’ function, as the ‘crop’ function would internally extract the extent for this object. However, for clarity we first create an Extent object from the re-projected Study Area, and use this Extent object as an argument in the ‘crop’ function.

RasterLayer object Coordinate Reference System (CRS)

.

#==================================================================================
# Subset Rasters to Study Area
#==================================================================================

#----------------------------------------------------------------------------------
# Create Cookie Cutter for the Study Area
#----------------------------------------------------------------------------------
# Change Study Area spatial polygon CRS (EPSG:4326) to the Australian-specific CRS used for raster (EPSG:3577).

# Find Raster CRS
SPGC.2010q3.rl

.

## class       : RasterLayer 
## dimensions  : 135159, 141481, 19122430479  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : -1944645, 2299785, -4910195, -855425  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : C:/Users/uqbblanc/Documents/TERN/04b-DSDP_GitHub/Prep/Landscapes_AusCover-RemoteSensing/YasiEffectsonGCatMB/YasiEffectsonGCatMB_Tutorial/SPGC_2010q3.vrt 
## names       : SPGC_2010q3 
## values      : 0, 255  (min, max)

.

crs(SPGC.2010q3.rl)

.

## CRS arguments:
##  +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0
## +ellps=GRS80 +units=m +no_defs

.

# Double check that this is the EPSG:3577 CRS (= to the raster's CRS)
CRS("+init=epsg:3577")

.

## CRS arguments:
##  +init=epsg:3577 +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0
## +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0
## +units=m +no_defs

.

Reproject SpatialPolygons object and create an Extent object

.

# Transform the Spatial Polygon CRS (from EPSG:4326 to EPSG:3577)
StudyArea.SP.reprj = spTransform(StudyArea.SP, CRS("+init=epsg:3577"))

# Study Area Reprojected Extent
StudyArea.extent.reprj = extent(StudyArea.SP.reprj) 
StudyArea.extent.reprj

.

## class       : Extent 
## xmin        : 1484019 
## xmax        : 1490198 
## ymin        : -1996975 
## ymax        : -1988144

.

Crop (Spatially Subset) RasterLayer objects

.

Example 1: SPGC for the Winter of 2010

.

#----------------------------------------------------------------------------------
# Subset (i.e. crop) raster layer to area of interest
#----------------------------------------------------------------------------------
# Subset both Year-Season rasters
SPGC.StudyArea.2010q3.rl = crop(SPGC.2010q3.rl, StudyArea.extent.reprj)
SPGC.StudyArea.2010q3.rl

.

## class       : RasterLayer 
## dimensions  : 295, 206, 60770  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 1484025, 1490205, -1996985, -1988135  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : SPGC_2010q3 
## values      : 100, 197  (min, max)

.

Example 2: SPGC for the Winter of 2011

.

SPGC.StudyArea.2011q3.rl = crop(SPGC.2011q3.rl, StudyArea.extent.reprj)
SPGC.StudyArea.2011q3.rl

.

## class       : RasterLayer 
## dimensions  : 295, 206, 60770  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 1484025, 1490205, -1996985, -1988135  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : SPGC_2011q3 
## values      : 100, 196  (min, max)

.