R for Spatial Statistics 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.
Note that rgdal is supposed to go away in 2023. terra is a potential replacement (https://www.r-bloggers.com/2021/05/a-comparison-of-terra-and-raster-packages/)
Required Packages
Install the "raster" and "rgdal" packages. "Raster" contains a set of classes for working with raster data and "rgal" 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(raster)
library(rgdal)
Reading Raster Files
Once you have the libraries loaded it is easy to read an existing raster file.
AnnualMeanTemp=raster("C:/ProjectsModeling/Cali_5minute/bio_1_AnnualMeanTemp_2.img")
AnnualMeanTemp=raster("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=raster(FilePath)
Again, R Studio will open the file chooser dialog behind its main window so you may need to minimize it.
Saving Raster Files
It's also easy to write a raster to a new file.
FilePath=file.choose() writeRaster(TheRaster, FilePath)
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.
library(raster)
# Load a raster file
FilePath="C:\\Temp3\\LandSat_WGS83.tif"
raster1=raster(FilePath)
# Load a CSV that contains colmns of X and Y coordinate values
pointCoordinates=read.csv("C:\\Temp3\\LandSatCoordinates.csv")
coordinates(pointCoordinates)= ~ X+ Y # Convert to a coordinates object
# Extract the values using the coordinates
rasValue=raster::extract(raster1, pointCoordinates)
# Add the pixel values back into the data frame
pointCoordinates=data.frame(pointCoordinates,rasValue)
Converting rasters to data frames
The code below converts a raster that has been loaded into a data frame. If you wish to have coordinates attached to the data frame, you can use the coordinates(...) call.
######################################################### # Convert a raster to a data frame and a SpatialPoints ThePoints=as(TheRaster, "SpatialPoints") # 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
Creating rasters from data frames
The raster() 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=raster(ncols=180, nrows=180) # Create a bounding box in the same units as the raster data BoundingBox=extent(-125,-110,30,45) # XMin, XMax, YMin, YMax
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 the values from a column in a data frame are used to add the pixel values to the raseter. Note that this requires the rasterize library.
library(rasterize) # 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) # 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 TheRaster=raster(ncols=WidthInPixels, nrows=HeightInPixels) # Create a bounding box in the same units as the raster data BoundingBox=extent(XMin,XMax,YMin,YMax) # Add the bounding box to the raster as its extent extent(TheRaster)=BoundingBox # Specify the columns for the coordinate values coordinates(NewData)=~X+Y # Render the points into the raster using the prediction as the values for the pixels TheRaster=rasterize(NewData,TheRaster,NewData$ThePrediction) # 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.
VectorsToRaster=function(X,Y,Pixels,WidthInPixels=200)
{
library(rasterize)
# Find the spatial extent of the points in reference units (i.e. degrees, meters or feet)
XMin=min(X)
XMax=max(X)
YMin=min(Y)
YMax=max(Y)
# 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
TheRaster=raster(ncols=WidthInPixels, nrows=HeightInPixels)
# Create a bounding box in the same units as the raster data
BoundingBox=extent(XMin,XMax,YMin,YMax)
# Add the bounding box to the raster as its extent
extent(TheRaster)=BoundingBox
TheData=data.frame(X,Y,Pixels)
# Specify the columns for the coordinate values
coordinates(TheData)=~X+Y
# Render the points into the raster using the prediction as the values for the pixels
TheRaster=rasterize(TheData,TheRaster,TheData$Pixels)
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)