Working with rasters

R has the ability to read raster file formats but you'll need to install some additional packages. Once you read or create a raster, you can use it as you would other vectors and matricies.

The rgal library is no longer supported in the new versions of R and is requred by raster so we'll be using terra below. rgdal is no longer supported either. I have been able to load TIFF (tif) and IMAGE (img) files but I'm not sure if all types are support or other file formats.

Terra Classes

Terra provides some special classes for working with spatial data.

SpatialPointsDataFrame - a dataframe that has and x and y column specified
SpaRaster - a raster with an extent and spatial reference system

Required Packages

Install the "terra" packages. "terra" contains a set of classes for working with raster data and includes the GDAL OpenSource library for reading and writing Geospatial file formats. This includes TIFF, GeoTIFF, IMG, and over 100 others. Installing these packages may also load a series of other libraries and it will take a while. Progress bars and error messages may appear behind R Studio so you may want to minimize it or move it out of the way.

After installing the packages, load the libraries.

library(terra)
library(sp)

Reading Raster Files

Once you have the libraries loaded it is easy to read an existing raster file.

AnnualMeanTemp=rast("C:/ProjectsModeling/Cali_5minute/bio_1_AnnualMeanTemp_2.img")
AnnualMeanTemp=rast("C:/ProjectsModeling/Cali_5minute/DEM.tif")

You can also use a "File Chooser" dialog to open an existing file as follows.

FilePath=file.choose()
TheRaster=rast(FilePath)

Again, R Studio will open the file chooser dialog behind its main window so you may need to minimize it.

Note that the error messages in terra can be far from usful. When I passed an invalid raster I recieved: " NULL value passed as symbol address". Passing a file path that included a folder that did not exist gave me "Error: external pointer is not valid".

Saving Raster Files

It's also easy to write a raster to a new file.

FilePath=file.choose()
writeRaster(TheRaster, FilePath)

Cropping A Raster

Below is the code to get the existing extent from a raster and then update it and use it to crop the raster. Note that the extent is just a vector.

# create an extent object
TheExtent=ext(TheRaster)
TheExtent # print the extent

# constrain the extent to the new bounds of the raster
TheExtent[1]=-124.75 # west 
TheExtent[2]=-114 # east
TheExtent[3]=30 # south
TheExtent[4]=45 # north

#  you can also just create an extent in one line
TheExtent=ext(-124,-120,40,44) # extent to crop to in spatial reference units (e.g. degrees)

# Crop and plot the new raster
CroppedRaster=crop(TheRaster, TheExtent)
plot(CroppedRaster)

Projecting a Raster

Below is the code to project a raster to a new spatial reference. Note that we just have to specify a "Proj string". See the Proj website for more information.

ProjectedRaster=project(CroppedRaster, "+proj=utm +zone=10")
plot(ProjectedRaster)

Extracting Values from Raster Files

The code below will extract the pixel values from a raster based on the coordinates in a data frame and then add the pixel values back into the data frame. Note that you can also pass a shapefile object directly into extract().

library(raster)

# Load a raster file
FilePath="C:\\Temp3\\LandSat_WGS83.tif"
raster1=rast(FilePath)

# Load a CSV that contains colmns of X and Y coordinate values
pointCoordinates=read.csv("C:\\Temp3\\LandSatCoordinates.csv")

# create a matrix of two columns, first with the x(lon) values and then the y(lat) values
xy=c(TheData$LON,TheData$LAT) 
xy=matrix(xy, length(TheData$LAT), 2)  

# extract the raster values using the matrix with the coordinates       
rasValue=terra::extract(AnnualMeanTemp,xy) 

# Add the pixel values back into the data frame using
pointCoordinates=data.frame(pointCoordinates,rasValue)

# This also works and may be a little more clear
pointCoordinates$Raster1=rasValue$Layer_1

Converting rasters to data frames

The code below converts a raster that has been loaded into a data frame. You can have the x,y coordinate values added by passing in "xy=TRUE".

#########################################################
# Convert a raster to a SpatVector object
ThePoints=as.points(TheRaster)

# Convert the raster to a data frame with x and y coordinate values in columns
TheDataFrame = as.data.frame(TheRaster,xy=TRUE)

# Convert data frame to a SpatialPointsDataFrame
coordinates(TheDataFrame)=~x+y

Once you have the raster in a data frame, you can access it's data as a column in the data frame. Below is an example on how to extract a vector with the pixel values from a raster i a file.

# Load the raster data
AnnualPrecipFilePath=bio_12_AnnualPrecip_CONUS_1.img"
AnnualPrecipRaster = rast(AnnualPrecipFilePath) # loads data as a spatraster (terra)

# Convert the raster to a data frame
AnnualPrecipDataFrame=as.data.frame(AnnualPrecipRaster,xy=TRUE)

# Get the raster data as the first 
AnnualPrecipVector=AnnualPrecipDataFrame$Layer_1 # get the vector of pixel values

Creating rasters from data frames

If you created a dataset from a raster, you can convert the dataset back to a rater by just using the following.

# create raster from dataframe
TheRaster2=rast(RasterDataFrame)
plot(TheRaster2)

If you are going to create a new raster from a dataset, things get more comlicated.

The rast() function creates a new raster and the extent() function sets the bounds of a raster in spatial reference units (e.g. degrees, meters, or feet). The code below shows how these functions are called.

# Create the raster with a specified width and height in pixels
TheRaster=rast(ncols=180, nrows=180)

The code below creates a new raster from a data frame. The spatial bounds of the raster is found from the X and Y values in the data frame. Then, the width of the raster is specified and height of the raster is computed from the width to keep the aspect ratio correct (i.e. pixels are the same width and height). Increasing the width of the raster in pixels will increase the resolution (and data size) of the raster. Then a matrix is created based on the width and height.

Note that this only works if the matrix includes all the pixel values. If you create a data frame or vector from a raster with no data and then remove the no data values, this approach will not work.


# Find the spatial extent of the points in reference units (i.e. degrees, meters or feet)
XMin=min(TheData$X)
XMax=max(TheData$X)
YMin=min(TheData$Y)
YMax=max(TheData$Y)

TheExtent=c(XMin,XMax,YMin,YMax)

# Set the width in pixels of the raster
WidthInPixels=200

# Find the height in pixels of the raster from the width and extent
HeightInPixels=WidthInPixels/(XMax-XMin)*(YMax-YMin)

# Create the raster with a specified width and height in pixels
TheMatrix=matrix(VectorOfPixelValues, nrow=HeightInPixels, ncol=WidthInPixels)

# Render the points into the raster using the prediction as the values for the pixels
TheRaster=rast(TheMatrix,TheExtent)

# Plot the raster
plot(TheRaster)

# Write the raster to a file
writeRaster(TheRaster, filename= "C:/Temp/Test.tif",overwrite=TRUE)

Below is a function that encapsolates the code above. The function takes vectors for the X coordinates, Y coordinates, and the Pixel values as parameters. An optional WidthInPixels parameter allows control over the resolution of the raster.

Note: This has not been tested yet

VectorsToRaster=function(X,Y,Pixels,WidthInPixels=200)
{
	# Find the spatial extent of the points in reference units (i.e. degrees, meters or feet)
	XMin=min(TheData$X)
	XMax=max(TheData$X)
	YMin=min(TheData$Y)
	YMax=max(TheData$Y)

	TheExtent=c(XMin,XMax,YMin,YMax)

	# Set the width in pixels of the raster
	WidthInPixels=200

	# Find the height in pixels of the raster from the width and extent
	HeightInPixels=WidthInPixels/(XMax-XMin)*(YMax-YMin)

	# Create the raster with a specified width and height in pixels
	TheMatrix=matrix(VectorOfPixelValues, nrow=HeightInPixels, ncol=WidthInPixels)

	# Render the points into the raster using the prediction as the values for the pixels
	TheRaster=rast(TheMatrix,TheExtent)

	return(TheRaster)
 }

The code below will call the function above to create a new raster from vectors in a data frame and then plot it and save it to a file.

# Create a raster that is 300 pixels across
TheRaster=VectorsToRaster(TheData$X,TheData$Y,TheData$MaxTemp,300)

# Plot the raster
plot(TheRaster)
		
# Write the raster to a file
writeRaster(TheRaster, filename= "C:/Temp/Test.tif",overwrite=TRUE)

Resources

Terra

R for Spatial Statistics