Sometimes we encounter raster datasets that do not "line up" when plotted or analyzed. Rasters that don't line up are most often in different Coordinate Reference Systems (CRS).
This tutorial explains how to deal with rasters in different, known CRSs. It
will walk though reprojecting rasters in R using the
function in the
After completing this tutorial, you will be able to:
- Be able to reproject a raster in R.
Things You’ll Need To Complete This Tutorial
You will need the most current version of R and, preferably,
on your computer to complete this tutorial.
Install R Packages
More on Packages in R – Adapted from Software Carpentry.
Data to Download
The LiDAR and imagery data used to create this raster teaching data subset were collected over the National Ecological Observatory Network's Harvard Forest and San Joaquin Experimental Range field sites and processed at NEON headquarters. The entire dataset can be accessed by request from the NEON Data Portal.
Set Working Directory: This lesson assumes that you have set your working directory to the location of the downloaded and unzipped data subsets.
R Script & Challenge Code: NEON data lessons often contain challenges that reinforce learned skills. If available, the code for challenge solutions is found in the downloadable R script of the entire lesson, available in the footer of each lesson page.
Raster Projection in R
In the Plot Raster Data in R tutorial, we learned how to layer a raster file on top of a hillshade for a nice looking basemap. In this tutorial, all of our data were in the same CRS. What happens when things don't line up?
We will use the
rgdal packages in this tutorial.
# load raster package library(raster) library(rgdal) # set working directory to ensure R can find the file we wish to import wd <- "~/Git/data/" # this will depend on your local environment # be sure that the downloaded file is in this directory setwd(wd)
Let's create a map of the Harvard Forest Digital Terrain Model
DTM_HARV) draped or layered on top of the hillshade (
# import DTM DTM_HARV <- raster(paste0(wd,"NEON-DS-Airborne-Remote-Sensing/HARV/DTM/HARV_dtmCrop.tif")) # import DTM hillshade DTM_hill_HARV <- raster(paste0(wd,"NEON-DS-Airborne-Remote-Sensing/HARV/DTM/HARV_DTMhill_WGS84.tif")) # plot hillshade using a grayscale color ramp plot(DTM_hill_HARV, col=grey(1:100/100), legend=FALSE, main="DTM Hillshade\n NEON Harvard Forest Field Site") # overlay the DTM on top of the hillshade plot(DTM_HARV, col=terrain.colors(10), alpha=0.4, add=TRUE, legend=FALSE)
Our results are curious - the Digital Terrain Model (
DTM_HARV) did not plot on
top of our hillshade. The hillshade plotted just fine on it's own. Let's try to
plot the DTM on it's own to make sure there are data there.
Code Tip: For boolean R elements, such as
add=TRUE, you can use
F in place of
# Plot DTM plot(DTM_HARV, col=terrain.colors(10), alpha=1, legend=F, main="Digital Terrain Model\n NEON Harvard Forest Field Site")
Our DTM seems to contain data and plots just fine. Let's next check the Coordinate Reference System (CRS) and compare it to our hillshade.
# view crs for DTM crs(DTM_HARV) ## CRS arguments: ## +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs # view crs for hillshade crs(DTM_hill_HARV) ## CRS arguments: +proj=longlat +datum=WGS84 +no_defs
DTM_HARV is in the UTM projection.
DTM_hill_HARV is in
Geographic WGS84 - which is represented by latitude and longitude values.
Because the two rasters are in different CRSs, they don't line up when plotted
in R. We need to reproject
DTM_hill_HARV into the UTM CRS. Alternatively,
we could project
DTM_HARV into WGS84.
We can use the
projectRaster function to reproject a raster into a new CRS.
Keep in mind that reprojection only works when you first have a defined CRS
for the raster object that you want to reproject. It cannot be used if no
CRS is defined. Lucky for us, the
DTM_hill_HARV has a defined CRS.
Data Tip: When we reproject a raster, we move it from one "grid" to another. Thus, we are modifying the data! Keep this in mind as we work with raster data.
To use the
projectRaster function, we need to define two things:
- the object we want to reproject and
- the CRS that we want to reproject it to.
The syntax is
We want the CRS of our hillshade to match the
DTM_HARV raster. We can thus
assign the CRS of our
DTM_HARV to our hillshade within the
function as follows:
# reproject to UTM DTM_hill_UTMZ18N_HARV <- projectRaster(DTM_hill_HARV, crs=crs(DTM_HARV)) # compare attributes of DTM_hill_UTMZ18N to DTM_hill crs(DTM_hill_UTMZ18N_HARV) ## CRS arguments: ## +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs crs(DTM_hill_HARV) ## CRS arguments: +proj=longlat +datum=WGS84 +no_defs # compare attributes of DTM_hill_UTMZ18N to DTM_hill extent(DTM_hill_UTMZ18N_HARV) ## class : Extent ## xmin : 731397.3 ## xmax : 733205.3 ## ymin : 4712403 ## ymax : 4713907 extent(DTM_hill_HARV) ## class : Extent ## xmin : -72.18192 ## xmax : -72.16061 ## ymin : 42.52941 ## ymax : 42.54234
Notice in the output above that the
DTM_hill_UTMZ18N_HARV is now
UTM. However, the extent values of
DTM_hillUTMZ18N_HARV are different from
Challenge: Extent Change with CRS Change
Why do you think the two extents differ?
Deal with Raster Resolution
Let's next have a look at the resolution of our reprojected hillshade.
# compare resolution res(DTM_hill_UTMZ18N_HARV) ##  1.000 0.998
The output resolution of
DTM_hill_UTMZ18N_HARV is 1 x 0.998. Yet, we know that
the resolution for the data should be 1m x 1m. We can tell R to force our
newly reprojected raster to be 1m x 1m resolution by adding a line of code
# adjust the resolution DTM_hill_UTMZ18N_HARV <- projectRaster(DTM_hill_HARV, crs=crs(DTM_HARV), res=1) # view resolution res(DTM_hill_UTMZ18N_HARV) ##  1 1
Let's plot our newly reprojected raster.
# plot newly reprojected hillshade plot(DTM_hill_UTMZ18N_HARV, col=grey(1:100/100), legend=F, main="DTM with Hillshade\n NEON Harvard Forest Field Site") # overlay the DTM on top of the hillshade plot(DTM_HARV, col=terrain.colors(100), alpha=0.4, add=T, legend=F)
We have now successfully draped the Digital Terrain Model on top of our hillshade to produce a nice looking, textured map!
Challenge: Reproject, then Plot a Digital Terrain Model
Create a map of the
San Joaquin Experimental Range
field site using the
Reproject the data as necessary to make things line up!
## Warning in projectRaster(DTM_hill_SJER, crs = crs(DTM_SJER), res = 1): ## input and ouput crs are the same
If you completed the San Joaquin plotting challenge in the Plot Raster Data in R tutorial, how does the map you just created compare to that map?
Get Lesson Code
If you have questions or comments on this content, please contact us.Contact Us