diff --git a/episodes/05-raster-multi-band-in-r.Rmd b/episodes/05-raster-multi-band-in-r.Rmd index 9d6e0fac9..775042cfa 100644 --- a/episodes/05-raster-multi-band-in-r.Rmd +++ b/episodes/05-raster-multi-band-in-r.Rmd @@ -24,8 +24,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) ``` @@ -34,32 +33,35 @@ 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. +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. :::::::::::::::::::::::::::::::::::::::::::::::::: We introduced multi-band raster data in -[an earlier lesson](https://datacarpentry.org/organization-geospatial/01-intro-raster-data). This episode explores how to import and plot -a multi-band raster in -R. +[an earlier lesson](https://datacarpentry.org/organization-geospatial/01-intro-raster-data). +This episode explores how to import and plot a multi-band raster in R. ## Getting Started with Multi-Band Data in R -In this episode, the multi-band data that we are working with is imagery +In this episode, the multi-band data that we are working with is imagery collected using the [NEON Airborne Observation Platform](https://www.neonscience.org/data-collection/airborne-remote-sensing) high resolution camera over the [NEON Harvard Forest field site](https://www.neonscience.org/field-sites/field-sites-map/HARV). -Each RGB image is a 3-band raster. The same steps would apply to -working with a multi-spectral image with 4 or more bands - like Landsat imagery. +Each RGB image is a 3-band raster. The same steps would apply to working with a +multi-spectral image with 4 or more bands - like Landsat imagery. -If we read a RasterStack object into R using the `raster()` function, it only reads -in the first band. +By using the `rast()` function along with the `lyrs` parameter, we can read +specific raster bands (i.e. the first one); omitting this parater would read +instead all bands. ```{r read-single-band} -RGB_band1_HARV <- raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif") +RGB_band1_HARV <- + rast("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif", + lyrs = 1) ``` We need to convert this data to a data frame in order to plot it with `ggplot`. @@ -71,7 +73,7 @@ RGB_band1_HARV_df <- as.data.frame(RGB_band1_HARV, xy = TRUE) ```{r harv-rgb-band1} ggplot() + geom_raster(data = RGB_band1_HARV_df, - aes(x = x, y = y, alpha = layer)) + + aes(x = x, y = y, alpha = HARV_RGB_Ortho_1)) + coord_quickmap() ``` @@ -79,8 +81,8 @@ ggplot() + ## Challenge -View the attributes of this band. What are its dimensions, CRS, resolution, min and max values, -and band number? +View the attributes of this band. What are its dimensions, CRS, resolution, min +and max values, and band number? ::::::::::::::: solution @@ -91,10 +93,9 @@ RGB_band1_HARV ``` Notice that when we look at the attributes of this band, we see: -`band: 1 (of 3 bands)` +`dimensions : 2317, 3073, 1 (nrow, ncol, nlyr)` -This is R telling us that this particular raster object has more bands (3) -associated with it. +This is R telling us that we read only one its bands. @@ -106,32 +107,33 @@ associated with it. ## Data Tip -The number of bands associated with a -raster object can also be determined using the `nbands()` function: syntax is -`nbands(RGB_band1_HARV)`. +The number of bands associated with a raster's file can also be determined +using the `describe()` function: syntax is `describe(sources(RGB_band1_HARV))`. :::::::::::::::::::::::::::::::::::::::::::::::::: ### Image Raster Data Values -As we saw in the previous exercise, this raster contains values between 0 and 255. These values -represent degrees of brightness associated with the image band. In -the case of a RGB image (red, green and blue), band 1 is the red band. When -we plot the red band, larger numbers (towards 255) represent pixels with more -red in them (a strong red reflection). Smaller numbers (towards 0) represent -pixels with less red in them (less red was reflected). To -plot an RGB image, we mix red + green + blue values into one single color to -create a full color image - similar to the color image a digital camera creates. +As we saw in the previous exercise, this raster contains values between 0 and +255. These values represent degrees of brightness associated with the image +band. In the case of a RGB image (red, green and blue), band 1 is the red band. +When we plot the red band, larger numbers (towards 255) represent pixels with +more red in them (a strong red reflection). Smaller numbers (towards 0) +represent pixels with less red in them (less red was reflected). To plot an RGB +image, we mix red + green + blue values into one single color to create a full +color image - similar to the color image a digital camera creates. ### Import A Specific Band -We can use the `raster()` function to import specific bands in our raster object -by specifying which band we want with `band = N` (N represents the band number we -want to work with). To import the green band, we would use `band = 2`. +We can use the `rast()` function to import specific bands in our raster object +by specifying which band we want with `lyrs = N` (N represents the band number we +want to work with). To import the green band, we would use `lyrs = 2`. ```{r read-specific-band} -RGB_band2_HARV <- raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif", band = 2) +RGB_band2_HARV <- + rast("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif", + lyrs = 2) ``` We can convert this data to a data frame and plot the same way we plotted the red band: @@ -143,7 +145,7 @@ RGB_band2_HARV_df <- as.data.frame(RGB_band2_HARV, xy = TRUE) ```{r rgb-harv-band2} ggplot() + geom_raster(data = RGB_band2_HARV_df, - aes(x = x, y = y, alpha = layer)) + + aes(x = x, y = y, alpha = HARV_RGB_Ortho_2)) + coord_equal() ``` @@ -151,15 +153,16 @@ ggplot() + ## Challenge: Making Sense of Single Band Images -Compare the plots of band 1 (red) and band 2 (green). Is the forested area darker or lighter in band 2 (the green band) compared to band 1 (the red band)? +Compare the plots of band 1 (red) and band 2 (green). Is the forested area +darker or lighter in band 2 (the green band) compared to band 1 (the red band)? ::::::::::::::: solution ## Solution -We'd expect a *brighter* value for the forest in band 2 (green) than in -band 1 (red) because the leaves on trees of most often appear "green" - -healthy leaves reflect MORE green light than red light. +We'd expect a *brighter* value for the forest in band 2 (green) than in band 1 +(red) because the leaves on trees of most often appear "green" - healthy leaves +reflect MORE green light than red light. @@ -169,14 +172,14 @@ healthy leaves reflect MORE green light than red light. ## Raster Stacks in R -Next, we will work with all three image bands (red, green and blue) as an R -RasterStack object. We will then plot a 3-band composite, or full color, -image. +Next, we will work with all three image bands (red, green and blue) as an R +raster object. We will then plot a 3-band composite, or full color, image. -To bring in all bands of a multi-band raster, we use the`stack()` function. +To bring in all bands of a multi-band raster, we use the`rast()` function. ```{r intro-to-raster-stacks} -RGB_stack_HARV <- stack("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif") +RGB_stack_HARV <- + rast("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif") ``` Let's preview the attributes of our stack object: @@ -185,22 +188,16 @@ Let's preview the attributes of our stack object: RGB_stack_HARV ``` -We can view the attributes of each band in the stack in a single output: - -```{r} -RGB_stack_HARV@layers -``` - -If we have hundreds of bands, we can specify which band we'd like to view -attributes for using an index value: +We can view the attributes of each band in the stack in a single output. For +example, if we had hundreds of bands, we could specify which band we'd like to +view attributes for using an index value: ```{r} RGB_stack_HARV[[2]] ``` -We can also use the `ggplot` functions to plot the data in any layer -of our RasterStack object. Remember, we need to convert to a data -frame first. +We can also use the `ggplot` functions to plot the data in any layer of our +raster object. Remember, we need to convert to a data frame first. ```{r} RGB_stack_HARV_df <- as.data.frame(RGB_stack_HARV, xy = TRUE) @@ -236,16 +233,16 @@ To render a final three band, colored image in R, we use the `plotRGB()` functio This function allows us to: -1. Identify what bands we want to render in the red, green and blue regions. The - `plotRGB()` function defaults to a 1=red, 2=green, and 3=blue band order. However, - you can define what bands you'd like to plot manually. Manual definition of - bands is useful if you have, for example a near-infrared band and want to create - a color infrared image. +1. Identify what bands we want to render in the red, green and blue regions. + The `plotRGB()` function defaults to a 1=red, 2=green, and 3=blue band + order. However, you can define what bands you'd like to plot manually. + Manual definition of bands is useful if you have, for example a + near-infrared band and want to create a color infrared image. 2. Adjust the `stretch` of the image to increase or decrease contrast. -Let's plot our 3-band image. Note that we can use the `plotRGB()` -function directly with our RasterStack object (we don't need a -dataframe as this function isn't part of the `ggplot2` package). +Let's plot our 3-band image. Note that we can use the `plotRGB()` function +directly with our RasterStack object (we don't need a dataframe as this +function isn't part of the `ggplot2` package). ```{r plot-rgb-image} plotRGB(RGB_stack_HARV, @@ -258,15 +255,15 @@ the image might improve clarity and contrast using `stretch="lin"` or ![](fig/dc-spatial-raster/imageStretch_dark.jpg){alt='Image Stretch'} -When the range of pixel brightness values is closer to 0, a darker image is rendered by default. We can stretch -the values to extend to the full 0-255 range of potential values to increase the visual contrast of the image. +When the range of pixel brightness values is closer to 0, a darker image is +rendered by default. We can stretch the values to extend to the full 0-255 +range of potential values to increase the visual contrast of the image. ![](fig/dc-spatial-raster/imageStretch_light.jpg){alt='Image Stretch light'} -When the range of pixel brightness values is closer to 255, a -lighter image is rendered by default. We can stretch the values to extend to -the full 0-255 range of potential values to increase the visual contrast of -the image. +When the range of pixel brightness values is closer to 255, a lighter image is +rendered by default. We can stretch the values to extend to the full 0-255 +range of potential values to increase the visual contrast of the image. ```{r plot-rbg-image-linear} plotRGB(RGB_stack_HARV, @@ -282,18 +279,17 @@ plotRGB(RGB_stack_HARV, stretch = "hist") ``` -In this case, the stretch doesn't enhance the contrast our image significantly -given the distribution of reflectance (or brightness) values is distributed well -between 0 and 255. +In this case, the stretch doesn't enhance the contrast our image significantly +given the distribution of reflectance (or brightness) values is distributed +well between 0 and 255. ::::::::::::::::::::::::::::::::::::::: challenge ## Challenge - NoData Values -Let's explore what happens with NoData values when working with -RasterStack objects and using the -`plotRGB()` function. We will use the -`HARV_Ortho_wNA.tif` GeoTIFF file in the +Let's explore what happens with NoData values when working with RasterStack +objects and using the `plotRGB()` function. We will use the +`HARV_Ortho_wNA.tif` GeoTIFF file in the `NEON-DS-Airborne-Remote-Sensing/HARVRGB_Imagery/` directory. 1. View the files attributes. Are there `NoData` values assigned for this file? @@ -303,29 +299,29 @@ RasterStack objects and using the 5. Plot the object as a true color image. 6. What happened to the black edges in the data? 7. What does this tell us about the difference in the data structure between - `HARV_Ortho_wNA.tif` and `HARV_RGB_Ortho.tif` (R object `RGB_stack`). How can + `HARV_Ortho_wNA.tif` and `HARV_RGB_Ortho.tif` (R object `RGB_stack`). How can you check? ::::::::::::::: solution ## Answers -1) First we use the `GDALinfo()` function to view the - data attributes. +1) First we use the `describe()` function to view the data attributes. ```{r} -GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_Ortho_wNA.tif") +describe("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_Ortho_wNA.tif") ``` -2) From the output above, we see that there are `NoData` values - and they are assigned the value of -9999. +2) From the output above, we see that there are `NoData` values and they are +assigned the value of -9999. 3) The data has three bands. -4) To read in the file, we will use the `stack()` function: +4) To read in the file, we will use the `rast()` function: ```{r} -HARV_NA <- stack("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_Ortho_wNA.tif") +HARV_NA <- + rast("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_Ortho_wNA.tif") ``` 5) We can plot the data with the `plotRGB()` function: @@ -336,12 +332,14 @@ plotRGB(HARV_NA, ``` 6) The black edges are not plotted. -7) Both data sets have `NoData` values, however, in the RGB\_stack the NoData value is not - defined in the tiff tags, thus R renders them as black as the reflectance - values are 0. The black edges in the other file are defined as -9999 and R renders them as NA. + +7) Both data sets have `NoData` values, however, in the RGB\_stack the NoData +value is not defined in the tiff tags, thus R renders them as black as the +reflectance values are 0. The black edges in the other file are defined as +-9999 and R renders them as NA. ```{r} -GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif") +describe("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif") ``` ::::::::::::::::::::::::: @@ -352,59 +350,50 @@ GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.t ## Data Tip -We can create a RasterStack from -several, individual single-band GeoTIFFs too. We will do this in -a later episode, +We can create a raster object from several, individual single-band GeoTIFFs +too. We will do this in a later episode, [Raster Time Series Data in R](12-time-series-raster/). :::::::::::::::::::::::::::::::::::::::::::::::::: -## RasterStack vs RasterBrick in R +## SpatRaster in R + +The R SpatRaster object type can handle rasters with multiple bands. +The SpatRaster only holds parameters that describe the properties of raster +data that is located somewhere on our computer. -The R RasterStack and RasterBrick object types can both store multiple bands. -However, how they store each band is different. The bands in a RasterStack are -stored as links to raster data that is located somewhere on our computer. A -RasterBrick contains all of the objects stored within the actual R object. -In most cases, we can work with a RasterBrick in the same way we might work -with a RasterStack. However a RasterBrick is often more efficient and faster -to process - which is important when working with larger files. +A SpatRasterDataset object can hold references to sub-datasets, that is, +SpatRaster objects. In most cases, we can work with a SpatRaster in the same +way we might work with a SpatRasterDataset. ::::::::::::::::::::::::::::::::::::::::: callout ## More Resources -You can read the help for the `brick()` function by typing `?brick`. +You can read the help for the `rast()` and `sds()` functions by typing `?rast` +or `?sds`. :::::::::::::::::::::::::::::::::::::::::::::::::: -We can turn a RasterStack into a RasterBrick in R by using -`brick(StackName)`. Let's use the `object.size()` function to compare RasterStack and RasterBrick objects. First we will check -the size of our RasterStack object: -```{r raster-bricks} -object.size(RGB_stack_HARV) -``` - -Now we will create a RasterBrick object from our RasterStack data and view its size: +We can build a SpatRasterDataset using a SpatRaster or a list of SpatRaster: ```{r} -RGB_brick_HARV <- brick(RGB_stack_HARV) - -object.size(RGB_brick_HARV) +RGB_sds_HARV <- sds(RGB_stack_HARV) +RGB_sds_HARV <- sds(list(RGB_stack_HARV, RGB_stack_HARV)) ``` -Notice that in the RasterBrick, all of the bands are stored within the actual -object. Thus, the RasterBrick object size is much larger than the -RasterStack object. +We can retrieve the SpatRaster objects from a SpatRasterDataset using +subsetting: -You use the `plotRGB()` function to plot a RasterBrick too: - -```{r plot-brick} -plotRGB(RGB_brick_HARV) +```{r} +RGB_sds_HARV[[1]] +RGB_sds_HARV[[2]] ``` + ::::::::::::::::::::::::::::::::::::::: challenge ## Challenge: What Functions Can Be Used on an R Object of a particular class? @@ -414,7 +403,7 @@ We can view various functions (or methods) available to use on an R object with 1. What methods can be used on the `RGB_stack_HARV` object? 2. What methods can be used on a single band within `RGB_stack_HARV`? -3. Why do you think there is a difference? +3. Why do you think there isn't a difference? ::::::::::::::: solution @@ -430,11 +419,10 @@ methods(class=class(RGB_stack_HARV)) 2) And compare that with the methods available for a single band: ```{r} -methods(class=class(RGB_stack_HARV[1])) +methods(class=class(RGB_stack_HARV[[1]])) ``` -3) There are far more things one could or want to ask of a full stack than of - a single band. +3) A SpatRaster is the same no matter its number of bands. @@ -447,8 +435,8 @@ methods(class=class(RGB_stack_HARV[1])) :::::::::::::::::::::::::::::::::::::::: keypoints - A single raster file can contain multiple bands or layers. -- Use the `stack()` function to load all bands in a multi-layer raster file into R. -- Individual bands within a stack can be accessed, analyzed, and visualized using the same functions as single bands. +- Use the `rast()` function to load all bands in a multi-layer raster file into R. +- Individual bands within a SpatRaster can be accessed, analyzed, and visualized using the same functions no matter how many bands it holds. ::::::::::::::::::::::::::::::::::::::::::::::::::