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

Author: albhasan

episodes/13-plot-time-series-rasters-in-r.Rmd

@@ -23,8 +23,7 @@ source("setup.R")
 ::::::::::::::::::::::::::::::::::::::::::::::::::
 
 ```{r load-libraries, echo=FALSE, results="hide", message=FALSE}
-library(raster)
-library(rgdal)
+library(terra)
 library(ggplot2)
 library(dplyr)
 library(reshape)
@@ -35,10 +34,13 @@ library(scales)
 ```{r load-data, echo=FALSE, results="hide"}
 # learners will have this data loaded from the previous episode
 
-all_NDVI_HARV <- list.files("data/NEON-DS-Landsat-NDVI/HARV/2011/NDVI", full.names = TRUE, pattern = ".tif$")
+all_NDVI_HARV <- list.files("data/NEON-DS-Landsat-NDVI/HARV/2011/NDVI", 
+                            full.names = TRUE, pattern = ".tif$")
 
 # Create a time series raster stack
-NDVI_HARV_stack <- stack(all_NDVI_HARV)
+NDVI_HARV_stack <- rast(all_NDVI_HARV)
+# NOTE: Fix the bands' names so they don't start with a number!
+names(NDVI_HARV_stack) <- paste0("X", names(NDVI_HARV_stack))
 
 # apply scale factor
 NDVI_HARV_stack <- NDVI_HARV_stack/10000
@@ -54,20 +56,22 @@ NDVI_HARV_stack_df <- as.data.frame(NDVI_HARV_stack, xy = TRUE) %>%
 
 ## 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.
+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.
 
 
 ::::::::::::::::::::::::::::::::::::::::::::::::::
 
-This episode covers how to customize your raster plots using the `ggplot2` package
-in R to create publication-quality plots.
+This episode covers how to customize your raster plots using the `ggplot2` 
+package in R to create publication-quality plots.
 
 ## Before and After
 
-In [the previous episode](12-time-series-raster/), we learned how to plot multi-band
-raster data in R using the `facet_wrap()` function. This created a separate panel in our plot
-for each raster band. The plot we created together is shown below:
+In [the previous episode](12-time-series-raster/), we learned how to plot 
+multi-band raster data in R using the `facet_wrap()` function. This created a 
+separate panel in our plot for each raster band. The plot we created together 
+is shown below:
 
 ```{r levelplot-time-series-before, echo=FALSE}
 # code not shown, demonstration only
@@ -77,10 +81,13 @@ ggplot() +
   ggtitle("Landsat NDVI", subtitle = "NEON Harvard Forest")
 ```
 
-Although this plot is informative, it isn't something we would expect to see in a journal publication. The x
-and y-axis labels aren't informative. There is a lot of unnecessary gray background and the titles of
-each panel don't clearly state that the number refers to the Julian day the data was collected. In this
-episode, we will customize this plot above to produce a publication quality graphic. We will go through these steps iteratively. When we're done, we will have created the plot shown below.
+Although this plot is informative, it isn't something we would expect to see in 
+a journal publication. The x and y-axis labels aren't informative. There is a 
+lot of unnecessary gray background and the titles of each panel don't clearly 
+state that the number refers to the Julian day the data was collected. In this
+episode, we will customize this plot above to produce a publication quality 
+graphic. We will go through these steps iteratively. When we're done, we will 
+have created the plot shown below.
 
 ```{r levelplot-time-series-after, echo=FALSE}
 # code not shown, demonstration only
@@ -94,10 +101,12 @@ green_colors <- brewer.pal(9, "YlGn") %>%
 
 ggplot() +
   geom_raster(data = NDVI_HARV_stack_df , aes(x = x, y = y, fill = value)) +
-  facet_wrap(~variable, nrow = 3, ncol = 5, labeller = labeller(variable = labels_names)) +
+  facet_wrap(~variable, nrow = 3, ncol = 5, 
+             labeller = labeller(variable = labels_names)) +
   ggtitle("Landsat NDVI - Julian Days", subtitle = "Harvard Forest 2011") +
   theme_void() +
-  theme(plot.title = element_text(hjust = 0.5, face = "bold"), plot.subtitle = element_text(hjust = 0.5)) + 
+  theme(plot.title = element_text(hjust = 0.5, face = "bold"), 
+        plot.subtitle = element_text(hjust = 0.5)) + 
   scale_fill_gradientn(name = "NDVI", colours = green_colors(20))
 
 # cleanup
@@ -106,8 +115,9 @@ rm(raster_names, labels_names, green_colors)
 
 ## Adjust the Plot Theme
 
-The first thing we will do to our plot remove the x and y-axis labels and axis ticks, as these are
-unnecessary and make our plot look messy. We can do this by setting the plot theme to `void`.
+The first thing we will do to our plot remove the x and y-axis labels and axis 
+ticks, as these are unnecessary and make our plot look messy. We can do this by 
+setting the plot theme to `void`.
 
 ```{r adjust-theme}
 ggplot() +
@@ -117,13 +127,16 @@ ggplot() +
   theme_void()
 ```
 
-Next we will center our plot title and subtitle. We need to do this **after** the `theme_void()` layer,
-because R interprets the `ggplot` layers in order. If we first tell R to center our plot title,
-and then set the theme to `void`, any adjustments we've made to the plot theme will be over-written
-by the `theme_void()` function. So first we make the theme `void` and then we center the title.
-We center both the title and subtitle by using the `theme()` function and setting the `hjust`
-parameter to 0.5. The `hjust` parameter stands for "horizontal justification" and takes any value between
-0 and 1. A setting of 0 indicates left justification and a setting of 1 indicates right justification.
+Next we will center our plot title and subtitle. We need to do this **after** 
+the `theme_void()` layer, because R interprets the `ggplot` layers in order. If 
+we first tell R to center our plot title, and then set the theme to `void`, any 
+adjustments we've made to the plot theme will be over-written by the 
+`theme_void()` function. So first we make the theme `void` and then we center 
+the title. We center both the title and subtitle by using the `theme()` 
+function and setting the `hjust` parameter to 0.5. The `hjust` parameter stands 
+for "horizontal justification" and takes any value between 0 and 1. A setting 
+of 0 indicates left justification and a setting of 1 indicates right 
+justification.
 
 ```{r adjust-theme-2}
 ggplot() +
@@ -139,9 +152,9 @@ ggplot() +
 
 ## Challenge
 
-Change the plot title (but not the subtitle) to bold font. You can
-(and should!) use the help menu in RStudio or any internet resources
-to figure out how to change this setting.
+Change the plot title (but not the subtitle) to bold font. You can (and 
+should!) use the help menu in RStudio or any internet resources to figure out 
+how to change this setting.
 
 :::::::::::::::  solution
 
@@ -170,12 +183,13 @@ ggplot() +
 Next, let's adjust the color ramp used to render the rasters. First, we can
 change the blue color ramp to a green one that is more visually suited to our
 NDVI (greenness) data using the `colorRampPalette()` function in combination
-with `colorBrewer` which requires loading the `RColorBrewer` library. Then we use `scale_fill_gradientn` to pass the list of
-colours (here 20 different colours) to ggplot.
+with `colorBrewer` which requires loading the `RColorBrewer` library. Then we 
+use `scale_fill_gradientn` to pass the list of colours (here 20 different 
+colours) to ggplot.
 
-First we need to create a set of colors to use. We will select a set of
-nine colors from the "YlGn" (yellow-green) color palette. This returns a
-set of hex color codes:
+First we need to create a set of colors to use. We will select a set of nine 
+colors from the "YlGn" (yellow-green) color palette. This returns a set of hex 
+color codes:
 
 ```{r}
 library(RColorBrewer)
@@ -190,9 +204,9 @@ green_colors <- brewer.pal(9, "YlGn") %>%
   colorRampPalette()
 ```
 
-We can
-tell the `colorRampPalette()` function how many discrete colors within this color range to
-create. In our case, we will use 20 colors when we plot our graphic.
+We can tell the `colorRampPalette()` function how many discrete colors within 
+this color range to create. In our case, we will use 20 colors when we plot our 
+graphic.
 
 ```{r change-color-ramp}
 ggplot() +
@@ -213,7 +227,8 @@ pixels that are more green have a higher NDVI value.
 
 ## Data Tip
 
-For all of the `brewer.pal` ramp names see the [brewerpal page](https://www.datavis.ca/sasmac/brewerpal.html).
+For all of the `brewer.pal` ramp names see the 
+[brewerpal page](https://www.datavis.ca/sasmac/brewerpal.html).
 
 
 ::::::::::::::::::::::::::::::::::::::::::::::::::
@@ -222,19 +237,19 @@ For all of the `brewer.pal` ramp names see the [brewerpal page](https://www.data
 
 ## Data Tip
 
-Cynthia Brewer, the creator of
-ColorBrewer, offers an online tool to help choose suitable color ramps, or to
-create your own. [ColorBrewer 2.0; Color Advise for Cartography](https://colorbrewer2.org/)
+Cynthia Brewer, the creator of ColorBrewer, offers an online tool to help 
+choose suitable color ramps, or to create your own. 
+[ColorBrewer 2.0; Color Advise for Cartography](https://colorbrewer2.org/)
 
 
 ::::::::::::::::::::::::::::::::::::::::::::::::::
 
 ## Refine Plot \& Tile Labels
 
-Next, let's label each panel in our plot with the Julian day that the
-raster data for that panel was collected. The current names come from the band
-"layer names"" stored in the
-`RasterStack` and the first part of each name is the Julian day.
+Next, let's label each panel in our plot with the Julian day that the raster 
+data for that panel was collected. The current names come from the band "layer 
+names"" stored in the `RasterStack` and the first part of each name is the 
+Julian day.
 
 To create a more meaningful label we can remove the "x" and replace it with
 "day" using the `gsub()` function in R. The syntax is as follows:
@@ -249,9 +264,9 @@ names(NDVI_HARV_stack)
 ```
 
 Now we will use the `gsub()` function to find the character string
-"\_HARV\_ndvi\_crop" and replace it with a blank string (""). We will
-assign this output to a new object (`raster_names`) and look
-at that object to make sure our code is doing what we want it to.
+"\_HARV\_ndvi\_crop" and replace it with a blank string (""). We will assign 
+this output to a new object (`raster_names`) and look at that object to make 
+sure our code is doing what we want it to.
 
 ```{r}
 raster_names <- names(NDVI_HARV_stack)
@@ -291,13 +306,14 @@ ggplot() +
 ## Change Layout of Panels
 
 We can adjust the columns of our plot by setting the number of columns `ncol`
-and the number of rows `nrow` in `facet_wrap`. Let's make our plot so that
-it has a width of five panels.
+and the number of rows `nrow` in `facet_wrap`. Let's make our plot so that it 
+has a width of five panels.
 
 ```{r adjust-layout}
 ggplot() +
   geom_raster(data = NDVI_HARV_stack_df , aes(x = x, y = y, fill = value)) +
-  facet_wrap(~variable, ncol = 5, labeller = labeller(variable = labels_names)) +
+  facet_wrap(~variable, ncol = 5, 
+             labeller = labeller(variable = labels_names)) +
   ggtitle("Landsat NDVI", subtitle = "NEON Harvard Forest") + 
   theme_void() + 
   theme(plot.title = element_text(hjust = 0.5, face = "bold"), 
@@ -311,18 +327,17 @@ Now we have a beautiful, publication quality plot!
 
 ## Challenge: Divergent Color Ramps
 
-When we used the `gsub()` function to modify the tile labels we replaced the beginning of each
-tile title with "Day". A more descriptive name could be "Julian Day". Update the
-plot above with the following changes:
+When we used the `gsub()` function to modify the tile labels we replaced the 
+beginning of each tile title with "Day". A more descriptive name could be 
+"Julian Day". Update the plot above with the following changes:
 
-1. Label each tile "Julian Day" with the julian day value
-  following.
+1. Label each tile "Julian Day" with the julian day value following.
 2. Change the color ramp to a divergent brown to green color ramp.
 
 **Questions:**
-Does having a divergent color ramp represent the data
-better than a sequential color ramp (like "YlGn")? Can you think of other data
-sets where a divergent color ramp may be best?
+Does having a divergent color ramp represent the data better than a sequential 
+color ramp (like "YlGn")? Can you think of other data sets where a divergent 
+color ramp may be best?
 
 :::::::::::::::  solution
 
@@ -345,8 +360,8 @@ ggplot() +
 ```
 
 For NDVI data, the sequential color ramp is better than the divergent as it is
-more akin to the process
-of greening up, which starts off at one end and just keeps increasing.
+more akin to the process of greening up, which starts off at one end and just 
+keeps increasing.