Lesson 3 updated by replacing calls to raster and rgdal packages with calls to terra package. #368 #363

Author: albhasan

episodes/03-raster-reproject-in-r.Rmd

@@ -22,8 +22,7 @@ source("setup.R")
 ::::::::::::::::::::::::::::::::::::::::::::::::::
 
 ```{r load-libraries, echo=FALSE, results="hide", message=FALSE, warning=FALSE}
-library(raster)
-library(rgdal)
+library(terra)
 library(ggplot2)
 library(dplyr)
 ```
@@ -33,16 +32,17 @@ library(dplyr)
 ## Things You'll Need To Complete This Episode
 
 See the [lesson homepage](.) for detailed information about the software,
-data, and other prerequisites you will need to work through the examples in this episode.
+data, and other prerequisites you will need to work through the examples in 
+this episode.
 
 
 ::::::::::::::::::::::::::::::::::::::::::::::::::
 
 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 episode explains how to deal with rasters in different, known CRSs. It
-will walk though reprojecting rasters in R using the `projectRaster()`
-function in the `raster` package.
+Reference Systems (CRS). This episode explains how to deal with rasters in 
+different, known CRSs. It will walk though reprojecting rasters in R using 
+the `project()` function in the `terra` package.
 
 ## Raster Projection in R
 
@@ -57,21 +57,23 @@ far in that the digital surface model (DSM) includes the tops of trees, while
 the digital terrain model (DTM) shows the ground level.
 
 We'll be looking at another model (the canopy height model) in
-[a later episode](04-raster-calculations-in-r/) and will see how to calculate the CHM from the
-DSM and DTM. Here, we will create a map of the Harvard Forest Digital
-Terrain Model
-(`DTM_HARV`) draped or layered on top of the hillshade (`DTM_hill_HARV`).
-The hillshade layer maps the terrain using light and shadow to create a 3D-looking image,
-based on a hypothetical illumination of the ground level.
+[a later episode](04-raster-calculations-in-r/) and will see how to calculate 
+the CHM from the DSM and DTM. Here, we will create a map of the Harvard Forest 
+Digital Terrain Model (`DTM_HARV`) draped or layered on top of the hillshade 
+(`DTM_hill_HARV`).
+The hillshade layer maps the terrain using light and shadow to create a 
+3D-looking image, based on a hypothetical illumination of the ground level.
 
 ![](fig/dc-spatial-raster/lidarTree-height.png){alt='Source: National Ecological Observatory Network (NEON).'}
 
 First, we need to import the DTM and DTM hillshade data.
 
 ```{r import-DTM-hillshade}
-DTM_HARV <- raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/DTM/HARV_dtmCrop.tif")
+DTM_HARV <- 
+    rast("data/NEON-DS-Airborne-Remote-Sensing/HARV/DTM/HARV_dtmCrop.tif")
 
-DTM_hill_HARV <- raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/DTM/HARV_DTMhill_WGS84.tif")
+DTM_hill_HARV <- 
+    rast("data/NEON-DS-Airborne-Remote-Sensing/HARV/DTM/HARV_DTMhill_WGS84.tif")
 ```
 
 Next, we will convert each of these datasets to a dataframe for
@@ -99,8 +101,7 @@ ggplot() +
 
 Our results are curious - neither the Digital Terrain Model (`DTM_HARV_df`)
 nor the DTM Hillshade (`DTM_hill_HARV_df`) plotted.
-Let's try to
-plot the DTM on its own to make sure there are data there.
+Let's try to plot the DTM on its own to make sure there are data there.
 
 ```{r plot-DTM}
 ggplot() +
@@ -123,10 +124,12 @@ geom_raster(data = DTM_hill_HARV_df,
     coord_quickmap()
 ```
 
-If we look at the axes, we can see that the projections of the two rasters are different.
-When this is the case, `ggplot` won't render the image. It won't even
-throw an error message to tell you something has gone wrong. We can look at Coordinate Reference Systems (CRSs) of the DTM and
-the hillshade data to see how they differ.
+If we look at the axes, we can see that the projections of the two rasters are 
+different.
+When this is the case, `ggplot` won't render the image. It won't even throw an 
+error message to tell you something has gone wrong. We can look at Coordinate 
+Reference Systems (CRSs) of the DTM and the hillshade data to see how they 
+differ.
 
 :::::::::::::::::::::::::::::::::::::::  challenge
 
@@ -141,10 +144,10 @@ does each use?
 
 ```{r explore-crs}
 # view crs for DTM
-crs(DTM_HARV)
+crs(DTM_HARV, parse = TRUE)
 
 # view crs for hillshade
-crs(DTM_hill_HARV)
+crs(DTM_hill_HARV, parse = TRUE)
 ```
 
 `DTM_HARV` is in the UTM projection, with units of meters.
@@ -158,12 +161,12 @@ crs(DTM_hill_HARV)
 ::::::::::::::::::::::::::::::::::::::::::::::::::
 
 Because the two rasters are in different CRSs, they don't line up when plotted
-in R. We need to reproject (or change the projection of) `DTM_hill_HARV` into the UTM CRS. Alternatively,
-we could reproject `DTM_HARV` into WGS84.
+in R. We need to reproject (or change the projection of) `DTM_hill_HARV` into 
+the UTM CRS. Alternatively, we could reproject `DTM_HARV` into WGS84.
 
 ## Reproject Rasters
 
-We can use the `projectRaster()` function to reproject a raster into a new CRS.
+We can use the `project()` 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.
@@ -172,46 +175,46 @@ 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.
+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:
+To use the `project()` function, we need to define two things:
 
 1. the object we want to reproject and
 2. the CRS that we want to reproject it to.
 
-The syntax is `projectRaster(RasterObject, crs = CRSToReprojectTo)`
+The syntax is `project(RasterObject, crs)`
 
 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 `projectRaster()`
-function as follows: `crs = crs(DTM_HARV)`.
-Note that we are using the `projectRaster()` function on the raster object,
+assign the CRS of our `DTM_HARV` to our hillshade within the `project()`
+function as follows: `crs(DTM_HARV)`.
+Note that we are using the `project()` function on the raster object,
 not the `data.frame()` we use for plotting with `ggplot`.
 
-First we will reproject our `DTM_hill_HARV` raster data to match the `DTM_HARV` raster CRS:
+First we will reproject our `DTM_hill_HARV` raster data to match the `DTM_HARV` 
+raster CRS:
 
 ```{r reproject-raster}
-DTM_hill_UTMZ18N_HARV <- projectRaster(DTM_hill_HARV,
-                                       crs = crs(DTM_HARV))
+DTM_hill_UTMZ18N_HARV <- project(DTM_hill_HARV,
+                                 crs(DTM_HARV))
 ```
 
-Now we can compare the CRS of our original DTM hillshade
-and our new DTM hillshade, to see how they are different.
+Now we can compare the CRS of our original DTM hillshade and our new DTM 
+hillshade, to see how they are different.
 
 ```{r}
-crs(DTM_hill_UTMZ18N_HARV)
-crs(DTM_hill_HARV)
+crs(DTM_hill_UTMZ18N_HARV, parse = TRUE)
+crs(DTM_hill_HARV, parse = TRUE)
 ```
 
 We can also compare the extent of the two objects.
 
 ```{r}
-extent(DTM_hill_UTMZ18N_HARV)
-extent(DTM_hill_HARV)
+ext(DTM_hill_UTMZ18N_HARV)
+ext(DTM_hill_HARV)
 ```
 
 Notice in the output above that the `crs()` of `DTM_hill_UTMZ18N_HARV` is now
@@ -228,8 +231,9 @@ Why do you think the two extents differ?
 
 ## Answers
 
-The extent for DTM\_hill\_UTMZ18N\_HARV is in UTMs so the extent is in meters. The extent for DTM\_hill\_HARV is in lat/long so the extent is expressed
-in decimal degrees.
+The extent for DTM\_hill\_UTMZ18N\_HARV is in UTMs so the extent is in meters. 
+The extent for DTM\_hill\_HARV is in lat/long so the extent is expressed in 
+decimal degrees.
 
 
 
@@ -239,31 +243,36 @@ in decimal degrees.
 
 ## Deal with Raster Resolution
 
-Let's next have a look at the resolution of our reprojected hillshade versus our original data.
+Let's next have a look at the resolution of our reprojected hillshade versus 
+our original data.
 
 ```{r view-resolution}
 res(DTM_hill_UTMZ18N_HARV)
 res(DTM_HARV)
 ```
 
-These two resolutions are different, but they're representing the same data. We can tell R to force our
-newly reprojected raster to be 1m x 1m resolution by adding a line of code
-`res=1` within the `projectRaster()` function. In the example below, we ensure a resolution match by using `res(DTM_HARV)` as a variable.
+These two resolutions are different, but they're representing the same data. We 
+can tell R to force our newly reprojected raster to be 1m x 1m resolution by 
+adding a line of code `res=1` within the `project()` function. In the 
+example below, we ensure a resolution match by using `res(DTM_HARV)` as a 
+variable.
 
 ```{r reproject-assign-resolution}
-  DTM_hill_UTMZ18N_HARV <- projectRaster(DTM_hill_HARV,
-                                         crs = crs(DTM_HARV),
-                                         res = res(DTM_HARV)) 
+  DTM_hill_UTMZ18N_HARV <- project(DTM_hill_HARV, 
+                                   crs(DTM_HARV), 
+                                   res = res(DTM_HARV)) 
 ```
 
-Now both our resolutions and our CRSs match, so we can plot these two data sets together. Let's double-check our resolution to be sure:
+Now both our resolutions and our CRSs match, so we can plot these two data sets 
+together. Let's double-check our resolution to be sure:
 
 ```{r}
 res(DTM_hill_UTMZ18N_HARV)
 res(DTM_HARV)
 ```
 
-For plotting with `ggplot()`, we will need to create a dataframe from our newly reprojected raster.
+For plotting with `ggplot()`, we will need to create a dataframe from our newly 
+reprojected raster.
 
 ```{r make-df-projected-raster}
 DTM_hill_HARV_2_df <- as.data.frame(DTM_hill_UTMZ18N_HARV, xy = TRUE)
@@ -302,15 +311,16 @@ Reproject the data as necessary to make things line up!
 
 ```{r challenge-code-reprojection, echo=TRUE}
 # import DSM
-DSM_SJER <- raster("data/NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_dsmCrop.tif")
+DSM_SJER <- 
+    rast("data/NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_dsmCrop.tif")
 # import DSM hillshade
 DSM_hill_SJER_WGS <-
-raster("data/NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_DSMhill_WGS84.tif")
+    rast("data/NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_DSMhill_WGS84.tif")
 
 # reproject raster
-DTM_hill_UTMZ18N_SJER <- projectRaster(DSM_hill_SJER_WGS,
-                                  crs = crs(DSM_SJER),
-                                  res = 1)
+DTM_hill_UTMZ18N_SJER <- project(DSM_hill_SJER_WGS,
+                                 crs(DSM_SJER),
+                                 res = 1)
 
 # convert to data.frames
 DSM_SJER_df <- as.data.frame(DSM_SJER, xy = TRUE)
@@ -355,7 +365,7 @@ is this one was reprojected from WGS84 to UTM prior to plotting.
 :::::::::::::::::::::::::::::::::::::::: keypoints
 
 - In order to plot two raster data sets together, they must be in the same CRS.
-- Use the `projectRaster()` function to convert between CRSs.
+- Use the `project()` function to convert between CRSs.
 
 ::::::::::::::::::::::::::::::::::::::::::::::::::