Skip to content

Commit 2cc2f7b

Browse files
authored
Merge pull request #384 from albhasan/m2022q1
Lesson 1 updated by replacing calls to raster and rgdal packages
2 parents 1be4ef4 + 763d842 commit 2cc2f7b

File tree

1 file changed

+65
-48
lines changed

1 file changed

+65
-48
lines changed

episodes/01-raster-structure.Rmd

Lines changed: 65 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -14,7 +14,7 @@ knitr::opts_chunk$set(fig.height = 6)
1414

1515
- Describe the fundamental attributes of a raster dataset.
1616
- Explore raster attributes and metadata using R.
17-
- Import rasters into R using the `raster` package.
17+
- Import rasters into R using the `terra` package.
1818
- Plot a raster file in R using the `ggplot2` package.
1919
- Describe the difference between single- and multi-band rasters.
2020

@@ -29,8 +29,7 @@ knitr::opts_chunk$set(fig.height = 6)
2929
::::::::::::::::::::::::::::::::::::::::::::::::::
3030

3131
```{r load-libraries, echo=FALSE, results="hide", message=FALSE, warning=FALSE}
32-
library(raster)
33-
library(rgdal)
32+
library(terra)
3433
library(ggplot2)
3534
library(dplyr)
3635
```
@@ -52,12 +51,13 @@ rasters in R, including CRS and resolution. We will also explore missing and bad
5251
data values as stored in a raster and how R handles these elements.
5352

5453
We will continue to work with the `dplyr` and `ggplot2` packages that were introduced
55-
in the [Introduction to R for Geospatial Data](https://datacarpentry.org/r-intro-geospatial/) lesson. We will use two additional packages in this episode to work with raster data - the
56-
`raster` and `rgdal` packages. Make sure that you have these packages loaded.
54+
in the [Introduction to R for Geospatial Data](https://datacarpentry.org/r-intro-geospatial/)
55+
lesson. We will use two additional packages in this episode to work with raster
56+
data - the `terra` and `sf` packages. Make sure that you have these packages
57+
loaded.
5758

5859
```{r load-libraries-2, eval=FALSE}
59-
library(raster)
60-
library(rgdal)
60+
library(terra)
6161
library(ggplot2)
6262
library(dplyr)
6363
```
@@ -81,19 +81,19 @@ page](https://datacarpentry.org/geospatial-workshop/data/).
8181

8282
We will be working with a series of GeoTIFF files in this lesson. The
8383
GeoTIFF format contains a set of embedded tags with metadata about the raster
84-
data. We can use the function `GDALinfo()` to get information about our raster
84+
data. We can use the function `describe()` to get information about our raster
8585
data before we read that data into R. It is ideal to do this before importing
8686
your data.
8787

8888
```{r view-attributes-gdal}
89-
GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")
89+
describe("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")
9090
```
9191

9292
If you wish to store this information in R, you can do the following:
9393

9494
```{r}
9595
HARV_dsmCrop_info <- capture.output(
96-
GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")
96+
describe("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")
9797
)
9898
```
9999

@@ -104,7 +104,7 @@ episode. By the end of this episode, you will be able to explain and understand
104104
## Open a Raster in R
105105

106106
Now that we've previewed the metadata for our GeoTIFF, let's import this
107-
raster dataset into R and explore its metadata more closely. We can use the `raster()`
107+
raster dataset into R and explore its metadata more closely. We can use the `rast()`
108108
function to open a raster in R.
109109

110110
::::::::::::::::::::::::::::::::::::::::: callout
@@ -123,7 +123,7 @@ First we will load our raster file into R and view the data structure.
123123

124124
```{r}
125125
DSM_HARV <-
126-
raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")
126+
rast("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")
127127
128128
DSM_HARV
129129
```
@@ -138,20 +138,17 @@ summary(DSM_HARV)
138138

139139
but note the warning - unless you force R to calculate these statistics using
140140
every cell in the raster, it will take a random sample of 100,000 cells and
141-
calculate from that instead. To force calculation on more, or even all values,
142-
you can use the parameter `maxsamp`:
141+
calculate from that instead. To force calculation all the values, you can use
142+
the function `values`:
143143

144144
```{r}
145-
summary(DSM_HARV, maxsamp = ncell(DSM_HARV))
145+
summary(values(DSM_HARV))
146146
```
147147

148-
You may not see major differences in summary stats as `maxsamp` increases,
149-
except with very large rasters.
150-
151148
To visualise this data in R using `ggplot2`, we need to convert it to a
152149
dataframe. We learned about dataframes in [an earlier
153150
lesson](https://datacarpentry.org/r-intro-geospatial/04-data-structures-part2/index.html).
154-
The `raster` package has an built-in function for conversion to a plotable dataframe.
151+
The `terra` package has an built-in function for conversion to a plotable dataframe.
155152

156153
```{r}
157154
DSM_HARV_df <- as.data.frame(DSM_HARV, xy = TRUE)
@@ -164,8 +161,12 @@ dataframe format.
164161
str(DSM_HARV_df)
165162
```
166163

167-
We can use `ggplot()` to plot this data. We will set the color scale to `scale_fill_viridis_c`
168-
which is a color-blindness friendly color scale. We will also use the `coord_quickmap()` function to use an approximate Mercator projection for our plots. This approximation is suitable for small areas that are not too close to the poles. Other coordinate systems are available in ggplot2 if needed, you can learn about them at their help page `?coord_map`.
164+
We can use `ggplot()` to plot this data. We will set the color scale to
165+
`scale_fill_viridis_c` which is a color-blindness friendly color scale. We will
166+
also use the `coord_quickmap()` function to use an approximate Mercator
167+
projection for our plots. This approximation is suitable for small areas that
168+
are not too close to the poles. Other coordinate systems are available in
169+
ggplot2 if needed, you can learn about them at their help page `?coord_map`.
169170

170171
```{r ggplot-raster, fig.cap="Raster plot with ggplot2 using the viridis color scale"}
171172
@@ -189,7 +190,7 @@ More information about the Viridis palette used above at
189190

190191
## Plotting Tip
191192

192-
For faster, simpler plots, you can use the `plot` function from the `raster` package.
193+
For faster, simpler plots, you can use the `plot` function from the `terra` package.
193194

194195
::::::::::::::: solution
195196

@@ -222,7 +223,7 @@ We can view the CRS string associated with our R object using the`crs()`
222223
function.
223224

224225
```{r view-resolution-units}
225-
crs(DSM_HARV)
226+
crs(DSM_HARV, proj = TRUE)
226227
```
227228

228229
::::::::::::::::::::::::::::::::::::::: challenge
@@ -254,7 +255,8 @@ and datum (`datum=`).
254255

255256
### UTM Proj4 String
256257

257-
Our projection string for `DSM_HARV` specifies the UTM projection as follows:
258+
A projection string (like the one of `DSM_HARV`) specifies the UTM projection
259+
as follows:
258260

259261
`+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0`
260262

@@ -269,7 +271,8 @@ Our projection string for `DSM_HARV` specifies the UTM projection as follows:
269271
Note that the zone is unique to the UTM projection. Not all CRSs will have a
270272
zone. Image source: Chrismurf at English Wikipedia, via [Wikimedia Commons](https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system#/media/File:Utm-zones-USA.svg) (CC-BY).
271273

272-
![The UTM zones across the continental United States. From: https://upload.wikimedia.org/wikipedia/commons/8/8d/Utm-zones-USA.svg](fig/Utm-zones-USA.svg)
274+
275+
![The UTM zones across the continental United States. From: https://upload.wikimedia.org/wikipedia/commons/8/8d/Utm-zones-USA.svg](fig/Utm-zones-USA.svg){alt='UTM zones in the USA.'}
273276

274277
## Calculate Raster Min and Max Values
275278

@@ -281,9 +284,11 @@ Raster statistics are often calculated and embedded in a GeoTIFF for us. We
281284
can view these values:
282285

283286
```{r view-min-max}
284-
minValue(DSM_HARV)
287+
minmax(DSM_HARV)
285288
286-
maxValue(DSM_HARV)
289+
min(values(DSM_HARV))
290+
291+
max(values(DSM_HARV))
287292
```
288293

289294
::::::::::::::::::::::::::::::::::::::::: callout
@@ -300,8 +305,8 @@ DSM_HARV <- setMinMax(DSM_HARV)
300305

301306
::::::::::::::::::::::::::::::::::::::::::::::::::
302307

303-
We can see that the elevation at our site ranges from `r raster::minValue(DSM_HARV)`m to
304-
`r raster::maxValue(DSM_HARV)`m.
308+
We can see that the elevation at our site ranges from `r min(terra::values(DSM_HARV))`m to
309+
`r max(terra::values(DSM_HARV))`m.
305310

306311
## Raster Bands
307312

@@ -311,12 +316,12 @@ raster: surface elevation in meters for one time period.
311316

312317
![](fig/dc-spatial-raster/single_multi_raster.png){alt='Multi-band raster image'}
313318

314-
A raster dataset can contain one or more bands. We can use the `raster()`
319+
A raster dataset can contain one or more bands. We can use the `rast()`
315320
function to import one single band from a single or multi-band raster. We can
316-
view the number of bands in a raster using the `nlayers()` function.
321+
view the number of bands in a raster using the `nly()` function.
317322

318323
```{r view-raster-bands}
319-
nlayers(DSM_HARV)
324+
nlyr(DSM_HARV)
320325
```
321326

322327
However, raster data can also be multi-band, meaning that one raster file
@@ -343,7 +348,7 @@ did not collect data in these areas.
343348
# no data demonstration code - not being taught
344349
# Use stack function to read in all bands
345350
RGB_stack <-
346-
stack("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif")
351+
rast("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif")
347352
348353
# aggregate cells from 0.25m to 2m for plotting to speed up the lesson and
349354
# save memory
@@ -389,19 +394,25 @@ ggplot() +
389394
390395
```
391396

392-
If your raster already has `NA` values set correctly but you aren't sure where they are, you can deliberately plot them in a particular colour. This can be useful when checking a dataset's coverage. For instance, sometimes data can be missing where a sensor could not 'see' its target data, and you may wish to locate that missing data and fill it in.
397+
If your raster already has `NA` values set correctly but you aren't sure where
398+
they are, you can deliberately plot them in a particular colour. This can be
399+
useful when checking a dataset's coverage. For instance, sometimes data can be
400+
missing where a sensor could not 'see' its target data, and you may wish to
401+
locate that missing data and fill it in.
393402

394-
To highlight `NA` values in ggplot, alter the `scale_fill_*()` layer to contain a colour instruction for `NA` values, like `scale_fill_viridis_c(na.value = 'deeppink')`
403+
To highlight `NA` values in ggplot, alter the `scale_fill_*()` layer to contain
404+
a colour instruction for `NA` values, like `scale_fill_viridis_c(na.value = 'deeppink')`
395405

396406
```{r, echo=FALSE}
397407
# demonstration code
398408
# function to replace 0 with NA where all three values are 0 only
399-
RGB_2m_nas <- calc(RGB_2m,
400-
fun = function(x) {
401-
x[rowSums(x == 0) == 3, ] <- rep(NA, nlayers(RGB_2m))
402-
x
403-
})
404-
RGB_2m_nas <- as.data.frame(RGB_2m_nas, xy = TRUE)
409+
RGB_2m_nas <- app(RGB_2m,
410+
fun = function(x) {
411+
if (sum(x == 0, na.rm = TRUE) == length(x))
412+
return(rep(NA, times = length(x)))
413+
x
414+
})
415+
RGB_2m_nas <- as.data.frame(RGB_2m_nas, xy = TRUE, na.rm = FALSE)
405416
406417
ggplot() +
407418
geom_raster(data = RGB_2m_nas, aes(x = x, y = y, fill = HARV_RGB_Ortho_3)) +
@@ -436,14 +447,15 @@ of `NA` will be ignored by R as demonstrated above.
436447

437448
## Challenge
438449

439-
Use the output from the `GDALinfo()` function to find out what `NoDataValue` is used for our `DSM_HARV` dataset.
450+
Use the output from the `describe()` and `sources()` functions to find out what
451+
`NoDataValue` is used for our `DSM_HARV` dataset.
440452

441453
::::::::::::::: solution
442454

443455
## Answers
444456

445457
```{r}
446-
GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")
458+
describe(sources(DSM_HARV))
447459
```
448460

449461
`NoDataValue` are encoded as -9999.
@@ -482,15 +494,20 @@ elevation values over 400m with a contrasting colour.
482494

483495
```{r demo-bad-data-highlighting, echo=FALSE, message=FALSE, warning=FALSE}
484496
# reclassify raster to ok/not ok
485-
DSM_highvals <- reclassify(DSM_HARV, rcl = c(0, 400, NA_integer_, 400, 420, 1L), include.lowest = TRUE)
497+
DSM_highvals <- classify(DSM_HARV,
498+
rcl = matrix(c(0, 400, NA_integer_, 400, 420, 1L),
499+
ncol = 3, byrow = TRUE),
500+
include.lowest = TRUE)
486501
DSM_highvals <- as.data.frame(DSM_highvals, xy = TRUE)
502+
487503
DSM_highvals <- DSM_highvals[!is.na(DSM_highvals$HARV_dsmCrop), ]
488504
489505
ggplot() +
490506
geom_raster(data = DSM_HARV_df, aes(x = x, y = y, fill = HARV_dsmCrop)) +
491507
scale_fill_viridis_c() +
492508
# use reclassified raster data as an annotation
493-
annotate(geom = 'raster', x = DSM_highvals$x, y = DSM_highvals$y, fill = scales::colour_ramp('deeppink')(DSM_highvals$HARV_dsmCrop)) +
509+
annotate(geom = 'raster', x = DSM_highvals$x, y = DSM_highvals$y,
510+
fill = scales::colour_ramp('deeppink')(DSM_highvals$HARV_dsmCrop)) +
494511
ggtitle("Elevation Data", subtitle = "Highlighting values > 400m") +
495512
coord_quickmap()
496513
@@ -532,7 +549,7 @@ no bad data values in this particular raster.
532549

533550
> ## Challenge: Explore Raster Metadata
534551
>
535-
> Use `GDALinfo()` to determine the following about the `NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_DSMhill.tif` file:
552+
> Use `describe()` to determine the following about the `NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_DSMhill.tif` file:
536553
>
537554
> 1. Does this file have the same CRS as `DSM_HARV`?
538555
> 2. What is the `NoDataValue`?
@@ -547,7 +564,7 @@ no bad data values in this particular raster.
547564
> > ```{r challenge-code-attributes}
548565
> > ```
549566
550-
GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV\_DSMhill.tif")
567+
describe("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV\_DSMhill.tif")
551568
552569
::::::::::::::::::::::::::::::::::::::: challenge
553570
@@ -569,7 +586,7 @@ GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV\_DSMhill.tif")
569586
570587
## More Resources
571588
572-
- [Read more about the `raster` package in R.](https://cran.r-project.org/package=raster)
589+
- [Read more about the `terra` package in R.](https://cran.r-project.org/package=terra)
573590
574591
575592
::::::::::::::::::::::::::::::::::::::::::::::::::

0 commit comments

Comments
 (0)