Import Raster Data

This Notebook will demonstrate how to import raster GIS data into R.

Import a Landsat 8 Image

Because Landsat images are multi-band we import them using the brick():

library(raster)

Loading required package: sp

yose_l8_rst <- raster::brick("./data/yose_l8_20180822_b2345.tif")
yose_l8_rst

class      : RasterBrick 
dimensions : 2555, 2014, 5145770, 4  (nrow, ncol, ncell, nlayers)
resolution : 30, 30  (x, y)
extent     : 246315, 306735, 4152885, 4229535  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs

Error in gregexpr(calltext, singleline, fixed = TRUE) : 
  regular expression is invalid UTF-8

 
source     : D:/Workshops/R-Spatial/rspatial_mod/outputs/rspatial_data/data/yose_l8_20180822_b2345.tif 
names      : yose_l8_20180822_b2345.1, yose_l8_20180822_b2345.2, yose_l8_20180822_b2345.3, yose_l8_20180822_b2345.4 
min values :                     7208,                     6100,                     5357,                     4330 
max values :                    29840,                    30960,                    32595,                    30481 

View a summary of the bands:

summary(yose_l8_rst)

summary is an estimate based on a sample of 1e+05 cells (1.94% of all cells)

        yose_l8_20180822_b2345.1 yose_l8_20180822_b2345.2 yose_l8_20180822_b2345.3 yose_l8_20180822_b2345.4
Min.                        7411                     6366                     5609                     4933
1st Qu.                     8589                     8072                     7749                    12439
Median                      9254                     8908                     8965                    13980
3rd Qu.                    10419                    10456                    10960                    15809
Max.                       26650                    27824                    29578                    27484
NA's                           0                        0                        0                        0

We can tell from the values that Landsat pixel values are 16-bit (0..65536). Remind ourselves of the band combinations for Landsat 8:

Plot it with plotRGB

Plot a Pseudo True Color:

plotRGB(yose_l8_rst, r=3, g=2, b=1, scale = 2 ^ 16)

Isn’t that awful? The problem is the range of brightness values is fairly narrow. The Landsat sensor isn’t like the camera in your phone, which artificially stretches values to make them look pleasing to the eye.

To make it look better, we can pass the stretch argument to have it make a contrast stretch:

plotRGB(yose_l8_rst, r=3, g=2, b=1, scale = 2 ^ 16, stretch = 'lin')

Plot a False Color Composite:

plotRGB(yose_l8_rst, r=4, g=3, b=2, scale = 2 ^ 16, stretch = 'lin')

You can do a lot more with satellite imagery using RStoolbox.

Plot it with tmap

To overlay the park boundary, we’ll use tmap:

library(sf)

Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1

library(tmap)

Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

epsg_utm11n_wgs84 <- 32611
yose_bnd_utm <- sf::st_read(dsn="./data", layer="yose_boundary") %>% 
  st_transform(epsg_utm11n_wgs84)

Reading layer `yose_boundary' from data source `D:\Workshops\R-Spatial\rspatial_mod\outputs\rspatial_data\data' using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 11 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -119.8864 ymin: 37.4947 xmax: -119.1964 ymax: 38.18515
Geodetic CRS:  North_American_Datum_1983

tm_shape(yose_l8_rst) +
  tm_rgb(r=4, g=3, b=2, interpolate = FALSE, max.value = 2 ^ 16) +
tm_shape(yose_bnd_utm) +
  tm_borders(col="red", lwd=2) +
tm_layout(legend.show = FALSE)

stars object downsampled to 888 by 1126 cells. See tm_shape manual (argument raster.downsample)

End

Congratulations, you’ve completed the Notebook!

To view your Notebook at HTML, save it (again), then click the ‘Preview’ button in the RStudio toolbar.