Skip to content

Commit

Permalink
Revert "Lesson 3 updated by replacing calls to raster and rgdal packa…
Browse files Browse the repository at this point in the history 
…ges with calls to terra package. datacarpentry#368 datacarpentry#363" because it was included by accident.
  • Loading branch information
@albhasan
albhasan committed Mar 14, 2023
1 parent 2461571 commit f7f6c2c
Showing 1 changed file with 59 additions and 69 deletions.
128 changes: 59 additions & 69 deletions episodes/03-raster-reproject-in-r.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,8 @@ source("setup.R")
::::::::::::::::::::::::::::::::::::::::::::::::::

```{r load-libraries, echo=FALSE, results="hide", message=FALSE, warning=FALSE}
library(terra)
library(raster)
library(rgdal)
library(ggplot2)
library(dplyr)
```
Expand All @@ -32,17 +33,16 @@ 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 `project()` function in the `terra` 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 `projectRaster()`
function in the `raster` package.

## Raster Projection in R

Expand All @@ -57,23 +57,21 @@ 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 <-
rast("data/NEON-DS-Airborne-Remote-Sensing/HARV/DTM/HARV_dtmCrop.tif")
DTM_HARV <- raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/DTM/HARV_dtmCrop.tif")
DTM_hill_HARV <-
rast("data/NEON-DS-Airborne-Remote-Sensing/HARV/DTM/HARV_DTMhill_WGS84.tif")
DTM_hill_HARV <- raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/DTM/HARV_DTMhill_WGS84.tif")
```

Next, we will convert each of these datasets to a dataframe for
Expand Down Expand Up @@ -101,7 +99,8 @@ 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() +
Expand All @@ -124,12 +123,10 @@ 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

Expand All @@ -144,10 +141,10 @@ does each use?

```{r explore-crs}
# view crs for DTM
crs(DTM_HARV, parse = TRUE)
crs(DTM_HARV)
# view crs for hillshade
crs(DTM_hill_HARV, parse = TRUE)
crs(DTM_hill_HARV)
```

`DTM_HARV` is in the UTM projection, with units of meters.
Expand All @@ -161,12 +158,12 @@ crs(DTM_hill_HARV, parse = TRUE)
::::::::::::::::::::::::::::::::::::::::::::::::::

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 `project()` function to reproject a raster into a new CRS.
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.
Expand All @@ -175,46 +172,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 `project()` function, we need to define two things:
To use the `projectRaster()` 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 `project(RasterObject, crs)`
The syntax is `projectRaster(RasterObject, crs = CRSToReprojectTo)`

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 `project()`
function as follows: `crs(DTM_HARV)`.
Note that we are using the `project()` function on the raster object,
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,
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 <- project(DTM_hill_HARV,
crs(DTM_HARV))
DTM_hill_UTMZ18N_HARV <- projectRaster(DTM_hill_HARV,
crs = 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, parse = TRUE)
crs(DTM_hill_HARV, parse = TRUE)
crs(DTM_hill_UTMZ18N_HARV)
crs(DTM_hill_HARV)
```

We can also compare the extent of the two objects.

```{r}
ext(DTM_hill_UTMZ18N_HARV)
ext(DTM_hill_HARV)
extent(DTM_hill_UTMZ18N_HARV)
extent(DTM_hill_HARV)
```

Notice in the output above that the `crs()` of `DTM_hill_UTMZ18N_HARV` is now
Expand All @@ -231,9 +228,8 @@ 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.



Expand All @@ -243,36 +239,31 @@ 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 `project()` 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 `projectRaster()` 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 <- project(DTM_hill_HARV,
crs(DTM_HARV),
res = res(DTM_HARV))
DTM_hill_UTMZ18N_HARV <- projectRaster(DTM_hill_HARV,
crs = 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)
Expand Down Expand Up @@ -311,16 +302,15 @@ Reproject the data as necessary to make things line up!

```{r challenge-code-reprojection, echo=TRUE}
# import DSM
DSM_SJER <-
rast("data/NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_dsmCrop.tif")
DSM_SJER <- raster("data/NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_dsmCrop.tif")
# import DSM hillshade
DSM_hill_SJER_WGS <-
rast("data/NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_DSMhill_WGS84.tif")
raster("data/NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_DSMhill_WGS84.tif")
# reproject raster
DTM_hill_UTMZ18N_SJER <- project(DSM_hill_SJER_WGS,
crs(DSM_SJER),
res = 1)
DTM_hill_UTMZ18N_SJER <- projectRaster(DSM_hill_SJER_WGS,
crs = crs(DSM_SJER),
res = 1)
# convert to data.frames
DSM_SJER_df <- as.data.frame(DSM_SJER, xy = TRUE)
Expand Down Expand Up @@ -365,7 +355,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 `project()` function to convert between CRSs.
- Use the `projectRaster()` function to convert between CRSs.

::::::::::::::::::::::::::::::::::::::::::::::::::

Expand Down

0 comments on commit f7f6c2c

Please sign in to comment.