Coordinate Reference Systems (CRSs) in Rasters


A Coordinate Reference Systems (CRS) is a coordinate-based scheme use to locate entities in geographical space. It is an essential aspect of spatial data.


CRS can be stored in many different formats. The most common formats are: 

  • 4:the default output from many spatial data R packages (e.g. ‘sp’, ‘rgdal’, ‘raster’)
  • EPSG codes: a more concise format, found for example in the package ‘sf’
  • Well Known Text (WKT): format that succinctly describes the elements of a CRS.


Set a CRS

Functions of R spatial packages can be used to set (and retrieve) the CRS of spatial object. Note that these functions add CRS information to spatial objects, but they do not change their projections. Thus, you need to be sure that you are assigning the correct CRS to your spatial object.

  • For Spatial* objects: ‘proj4string’ function in the ‘sp’ package
  • For Raster* objects:
    • ‘crs’ function in the ‘raster’ package. It has replaced the earlier function ‘projection’.
    • ‘projection’ function in the ‘raster’ package: Earlier function, still working. However, the use of the ‘crs’ function is preferred.
    • ‘proj4string’: For compatibility with ‘sp’, ‘proj4string’ can also be used instead of ‘crs’.


Change a CRS

The ‘projectRaster’ function in the ‘raster’ package is used to project the values of a Raster* object to a new Raster* object with another projection (i.e. CRS). This can be done by:

  • Providing the new projection as a single argument: The function sets the extent and resolution of the new object.
  • Providing a Raster* object with the properties that the input data should be projected to: This allows more control over the transformation (e.g. assuring that the new object lines up with other datasets).




EXAMPLES


Examples on how to set and change the CRS of raster/polygon objects are provided below. The examples are 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. For more details on CRS and CRS operations, see TERN.s DSDP "Using Raster Data in R" tutorial.



Define a CRS

Below we first create a SpatialPolygons* object for the Study Area in the “Effects of Cyclone Yasi on Green Cover at Mission Beach” tutorial. This polygon does not have a CRS, so we provide an appropriate one. We use the EPSG:4326, which is the EPSG code for geographic coordinates on the WGS84 datum. This is a general spatial reference (i.e. area: World; scope: Horizontal component of 3D system). It is used by GPS satellite navigation systems and NATO military geodetic surveying). For more details see this spatialreferece.org link.


#==================================================================================
# Study Area
#==================================================================================

# Define Spatial 'Extent' of Study Area: Mission Beach backdrop
# --------------------------------------------------------------
StudyArea.extent = extent(146.0333, 146.0833, -17.9333, -17.8583)

# Create a working spatial structure from the spatial 'Extent'
# ------------------------------------------------------------
# Create a 'Spatial Polygon' from Spatial 'Extent'
StudyArea.SP = as(StudyArea.extent, "SpatialPolygons")
StudyArea.SP # Has not CRS


## class       : SpatialPolygons 
## features    : 1 
## extent      : 146.0333, 146.0833, -17.9333, -17.8583  (xmin, xmax, ymin, ymax)
## coord. ref. : NA



# Add a CRS. Use: EPSG:4326 = 'WGS84' - General Spatial Reference (https://spatialreference.org/ref/epsg/wgs-84/).
proj4string(StudyArea.SP) = CRS("+init=epsg:4326") # Add CRS
StudyArea.SP # Has a EPSG:4326 CRS


## class       : SpatialPolygons 
## features    : 1 
## extent      : 146.0333, 146.0833, -17.9333, -17.8583  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0



Change a CRS

Below we reproject a set of raster objects from their Australian-specific CRS (EPSG:3577= ‘GDA94’ = ‘Australian Albers’) to a general world CRS (EPSG:4326 = ‘WGS 84’), so that it can better compared to 

downloaded Stamen Maps for the region.


# Reproject Rasters to EPSG:4326 (general spatial CRS with coordinates on the WGS84
# datum on the Stamen Maps)
# -----------------------------------------------------------------------------------
# SPGC Study Area 2010 Winter
SPGC.StudyArea.2010q3.reprj.rl = projectRaster(SPGC.StudyArea.2010q3.rl, crs=CRS("+init=epsg:4326"))
# SPGC Study Area 20101Winter
SPGC.StudyArea.2011q3.reprj.rl = projectRaster(SPGC.StudyArea.2011q3.rl, crs=CRS("+init=epsg:4326"))
# SGGC Study Area 2010 Winter
SGGC.StudyArea.2010q3.reprj.rl = projectRaster(SGGC.StudyArea.2010q3.rl, crs=CRS("+init=epsg:4326"))
# SGGC Study Area 2011 Winter
SGGC.StudyArea.2011q3.reprj.rl = projectRaster(SGGC.StudyArea.2011q3.rl, crs=CRS("+init=epsg:4326"))
# SPGC Study Area Difference between Winters of 2011 and 2010
SPGC.StudyArea.Diffq3.reprj.rl = projectRaster(SPGC.StudyArea.Diffq3.rl, crs=CRS("+init=epsg:4326"))
# SGGC Study Area Difference between Winters of 2011 and 2010
SGGC.StudyArea.Diffq3.reprj.rl = projectRaster(SGGC.StudyArea.Diffq3.rl, crs=CRS("+init=epsg:4326"))