Raster Data

Library Loading

Before we can dive into “actual” spatial work we’ll need to load some libraries.

Load terra to work with raster data.

# Load needed libraries
library(terra)

Load the rioxarray and rasterio libraries to work with raster data. We’ll also load the os library to deal with file path issues.

# Load needed libraries
import os
import rasterio as rio
import rioxarray as rxr

To demonstrate raster data operations we’ll use NASA’s Shuttle Radar Topography Mission Global 3 arc second elevation data. Note that some minor preparatory work was necessary to get the data ready for our purposes here and is preserved in this folder of the website’s GitHub repository.

Loading Rasters

We can begin by reading in the data.

There are several R packages for working with raster data but we’ll focus on terra.

# Load a raster of elevation data
rast_r <- terra::rast(x = file.path("data", "elevation.tif"))

# Check the class of this object
class(rast_r)
[1] "SpatRaster"
attr(,"package")
[1] "terra"

The rioxarray library will be our focal library for raster operations with Python.

# Load a raster of elevation data
rast_py = rxr.open_rasterio(os.path.join("data", "elevation.tif"), masked = True).squeeze()

# Check the type of this variable
type(rast_py.rio.crs)
<class 'rasterio.crs.CRS'>

Checking Raster CRS

As noted earlier, the very first thing we should do after reading in spatial data is check the coordinate reference system!

We can check the CRS with the crs function (from terra).

# Check the CRS of the elevation data
terra::crs(x = rast_r)
[1] "GEOGCRS[\"WGS 84\",\n    ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n        MEMBER[\"World Geodetic System 1984 (Transit)\"],\n        MEMBER[\"World Geodetic System 1984 (G730)\"],\n        MEMBER[\"World Geodetic System 1984 (G873)\"],\n        MEMBER[\"World Geodetic System 1984 (G1150)\"],\n        MEMBER[\"World Geodetic System 1984 (G1674)\"],\n        MEMBER[\"World Geodetic System 1984 (G1762)\"],\n        MEMBER[\"World Geodetic System 1984 (G2139)\"],\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ENSEMBLEACCURACY[2.0]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Horizontal component of 3D system.\"],\n        AREA[\"World.\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"

Python stores raster CRS as an attribute.

# Check the CRS of the elevation data
rast_py.rio.crs
CRS.from_wkt('GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]')

You may note in the above operations that Python and R seem to be returning different information for the same data. In fact they are saying the same thing just in slightly different ways. R is telling us the name of the coordinate reference system (World Geodetic System 1984) while Python is giving us the four-digit “EPSG” code. For context, the European Petroleum Survey Group (EPSG) is a major authority on accepted coordinate reference systems. Each CRS is given a unique, four-digit code. As you may now be able to guess, “4326” is the EPSG code for the World Geodetic System 1984 CRS!

This CRS (WGS84) is an especially common one so you may eventually commit its EPSG code to memory but often you’ll encounter either CRS names or EPSG codes with which you are not familiar. It may not sound like sage advice but it can be simplest to just use whichever piece of information your coding language returns to Google the “missing” information.

Transforming Raster CRS

Once we know the current CRS, we can transform into a different coordinate reference system as needed. Let’s transform from WGS84 into another commonly-used CRS: Albers Equal Area (EPSG code 3083).

Quick warning, transforming rasters can take a long time and is very computationally intensive. It’ll work for a relatively small raster like the one we’re working with here but in general you should be careful with attempting to re-project a raster into a different CRS.

For CRS transformations in R we can use the project function (pronounced like the verb not the noun).

# Transform the raster CRS
rast_alber_r <- terra::project(x = rast_r, y = "epsg:3083")

# Re-check the CRS to confirm it worked
terra::crs(rast_alber_r)
[1] "PROJCRS[\"NAD83 / Texas Centric Albers Equal Area\",\n    BASEGEOGCRS[\"NAD83\",\n        DATUM[\"North American Datum 1983\",\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4269]],\n    CONVERSION[\"Texas Centric Albers Equal Area\",\n        METHOD[\"Albers Equal Area\",\n            ID[\"EPSG\",9822]],\n        PARAMETER[\"Latitude of false origin\",18,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8821]],\n        PARAMETER[\"Longitude of false origin\",-100,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8822]],\n        PARAMETER[\"Latitude of 1st standard parallel\",27.5,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8823]],\n        PARAMETER[\"Latitude of 2nd standard parallel\",35,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8824]],\n        PARAMETER[\"Easting at false origin\",1500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8826]],\n        PARAMETER[\"Northing at false origin\",6000000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8827]]],\n    CS[Cartesian,2],\n        AXIS[\"easting (X)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"northing (Y)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"State-wide spatial data presentation requiring true area measurements.\"],\n        AREA[\"United States (USA) - Texas.\"],\n        BBOX[25.83,-106.66,36.5,-93.5]],\n    ID[\"EPSG\",3083]]"

Python CRS re-projections have slightly more steps to work through.

# Create a CRS variable for the desired ending CRS
wgs84_crs = rio.crs.CRS.from_string("EPSG:3083")

# Reproject the raster we had
rast_alber_py = rast_py.rio.reproject(wgs84_crs)

# Re-check the CRS to confirm it worked
rast_alber_py.rio.crs
CRS.from_wkt('PROJCS["NAD83 / Texas Centric Albers Equal Area",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["latitude_of_center",18],PARAMETER["longitude_of_center",-100],PARAMETER["standard_parallel_1",27.5],PARAMETER["standard_parallel_2",35],PARAMETER["false_easting",1500000],PARAMETER["false_northing",6000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","3083"]]')