Sentinel 2 Notebook
Ivan Lizarazo
22.10.2020
1. Introduction
This is an R Markdown notebook which illustrates how to read, subset, and process a Sentinel-2 multispectral image in R.
This notebook aims to help Geomatica Basica students at Universidad Nacional to get to know R geospatial functionalities.
Previously, I downloaded a Sentinel 2 image over Bogota from the GloVis website.
After uncompressing the file, the content of the downloaded file is:
2. Preliminaries
First, clean the memory:
rm(list=ls())Then, load the required libraries:
library(rgdal)
library(gdalUtils)library(raster)
library(sf)## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0##
## Attaching package: 'sf'## The following object is masked from 'package:gdalUtils':
##
## gdal_rasterizelibrary(sp)
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──## ✓ ggplot2 3.3.2 ✓ purrr 0.3.4
## ✓ tibble 3.0.4 ✓ dplyr 1.0.2
## ✓ tidyr 1.1.0 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x tidyr::extract() masks raster::extract()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x dplyr::select() masks raster::select()#library(getSpatialData)folder <- "./datos/S2B_MSIL1C_20200104T152639_N0208_R025_T18NWL_20200104T184611.SAFE/GRANULE/L1C_T18NWL_A014780_20200104T153054/IMG_DATA/"(bandas <- list.files(folder))## [1] "T18NWL_20200104T152639_B01.jp2" "T18NWL_20200104T152639_B02.jp2"
## [3] "T18NWL_20200104T152639_B03.jp2" "T18NWL_20200104T152639_B04.jp2"
## [5] "T18NWL_20200104T152639_B05.jp2" "T18NWL_20200104T152639_B06.jp2"
## [7] "T18NWL_20200104T152639_B07.jp2" "T18NWL_20200104T152639_B08.jp2"
## [9] "T18NWL_20200104T152639_B09.jp2" "T18NWL_20200104T152639_B10.jp2"
## [11] "T18NWL_20200104T152639_B11.jp2" "T18NWL_20200104T152639_B12.jp2"
## [13] "T18NWL_20200104T152639_B8A.jp2" "T18NWL_20200104T152639_TCI.jp2"bandas[1]## [1] "T18NWL_20200104T152639_B01.jp2"(ruta_bandas <- paste0(folder, bandas))## [1] "./datos/S2B_MSIL1C_20200104T152639_N0208_R025_T18NWL_20200104T184611.SAFE/GRANULE/L1C_T18NWL_A014780_20200104T153054/IMG_DATA/T18NWL_20200104T152639_B01.jp2"
## [2] "./datos/S2B_MSIL1C_20200104T152639_N0208_R025_T18NWL_20200104T184611.SAFE/GRANULE/L1C_T18NWL_A014780_20200104T153054/IMG_DATA/T18NWL_20200104T152639_B02.jp2"
## [3] "./datos/S2B_MSIL1C_20200104T152639_N0208_R025_T18NWL_20200104T184611.SAFE/GRANULE/L1C_T18NWL_A014780_20200104T153054/IMG_DATA/T18NWL_20200104T152639_B03.jp2"
## [4] "./datos/S2B_MSIL1C_20200104T152639_N0208_R025_T18NWL_20200104T184611.SAFE/GRANULE/L1C_T18NWL_A014780_20200104T153054/IMG_DATA/T18NWL_20200104T152639_B04.jp2"
## [5] "./datos/S2B_MSIL1C_20200104T152639_N0208_R025_T18NWL_20200104T184611.SAFE/GRANULE/L1C_T18NWL_A014780_20200104T153054/IMG_DATA/T18NWL_20200104T152639_B05.jp2"
## [6] "./datos/S2B_MSIL1C_20200104T152639_N0208_R025_T18NWL_20200104T184611.SAFE/GRANULE/L1C_T18NWL_A014780_20200104T153054/IMG_DATA/T18NWL_20200104T152639_B06.jp2"
## [7] "./datos/S2B_MSIL1C_20200104T152639_N0208_R025_T18NWL_20200104T184611.SAFE/GRANULE/L1C_T18NWL_A014780_20200104T153054/IMG_DATA/T18NWL_20200104T152639_B07.jp2"
## [8] "./datos/S2B_MSIL1C_20200104T152639_N0208_R025_T18NWL_20200104T184611.SAFE/GRANULE/L1C_T18NWL_A014780_20200104T153054/IMG_DATA/T18NWL_20200104T152639_B08.jp2"
## [9] "./datos/S2B_MSIL1C_20200104T152639_N0208_R025_T18NWL_20200104T184611.SAFE/GRANULE/L1C_T18NWL_A014780_20200104T153054/IMG_DATA/T18NWL_20200104T152639_B09.jp2"
## [10] "./datos/S2B_MSIL1C_20200104T152639_N0208_R025_T18NWL_20200104T184611.SAFE/GRANULE/L1C_T18NWL_A014780_20200104T153054/IMG_DATA/T18NWL_20200104T152639_B10.jp2"
## [11] "./datos/S2B_MSIL1C_20200104T152639_N0208_R025_T18NWL_20200104T184611.SAFE/GRANULE/L1C_T18NWL_A014780_20200104T153054/IMG_DATA/T18NWL_20200104T152639_B11.jp2"
## [12] "./datos/S2B_MSIL1C_20200104T152639_N0208_R025_T18NWL_20200104T184611.SAFE/GRANULE/L1C_T18NWL_A014780_20200104T153054/IMG_DATA/T18NWL_20200104T152639_B12.jp2"
## [13] "./datos/S2B_MSIL1C_20200104T152639_N0208_R025_T18NWL_20200104T184611.SAFE/GRANULE/L1C_T18NWL_A014780_20200104T153054/IMG_DATA/T18NWL_20200104T152639_B8A.jp2"
## [14] "./datos/S2B_MSIL1C_20200104T152639_N0208_R025_T18NWL_20200104T184611.SAFE/GRANULE/L1C_T18NWL_A014780_20200104T153054/IMG_DATA/T18NWL_20200104T152639_TCI.jp2"3. Reading individual bands
The common way of reading raster files from R is using the raster library functionalities: - raster (for individual bands) - stack (for merged bands)
getGDALVersionInfo()## [1] "GDAL 2.4.2, released 2019/06/28"Older GDAL versions (i.e. GDAL <= 2.4.4) are not capable of reading JPEG-2000 files. Thus, we need to translate such files into the very-well-known GeoTiff format:
gdal_translate(ruta_bandas[2], "./datos/s2_bogota/B02.tif")## NULLgdal_translate(ruta_bandas[3], "./datos/s2_bogota/B03.tif")gdal_translate(ruta_bandas[4], "./datos/s2_bogota/B04.tif")gdal_translate(ruta_bandas[8], "./datos/s2_bogota/B08.tif")After converting the jp2 files into tif files, the result is this: 
ruta_tif <- "./datos/s2_bogota/"(lista_tif <- paste0(ruta_tif, list.files(ruta_tif)))## [1] "./datos/s2_bogota/B02.tif" "./datos/s2_bogota/B03.tif"
## [3] "./datos/s2_bogota/B04.tif" "./datos/s2_bogota/B08.tif"Let’s create a raster stack using the raster library:
s2_stack <- stack(lista_tif)s2_stack## class : RasterStack
## dimensions : 10980, 10980, 120560400, 4 (nrow, ncol, ncell, nlayers)
## resolution : 10, 10 (x, y)
## extent : 499980, 609780, 490200, 6e+05 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## names : B02, B03, B04, B08
## min values : 0, 0, 0, 0
## max values : 65535, 65535, 65535, 655354. Subsetting raster stack
We need a file for our area of interest. For this, we can use the geojson.io application to create a polygon. The corresponding data was stored in a local directory.
mosca <- st_read("./datos/mosquera.geojson")## Reading layer `mosquera' from data source `/Users/ivanlizarazo/Documents/ivan/GB/notebooks/datos/mosquera.geojson' using driver `GeoJSON'
## Simple feature collection with 1 feature and 0 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: -74.24011 ymin: 4.636997 xmax: -74.17248 ymax: 4.703722
## CRS: 4326(crs_s2 <- crs(s2_stack))## CRS arguments:
## +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0Let’s reproject the mosca object into the same crs of the raster stack:
mosca2 <- st_transform(mosca, crs_s2) Now, convert the simple feature into a spatial polygon data frame:
aoi <- as_Spatial(mosca2)Now, subset the sentinel 2 stack to cover the area of interest:
s2_crop <- crop(s2_stack, aoi, filename = "./datos/s2_mosca.tif", overwrite=TRUE)s2_crop## class : RasterBrick
## dimensions : 737, 751, 553487, 4 (nrow, ncol, ncell, nlayers)
## resolution : 10, 10 (x, y)
## extent : 584280, 591790, 512590, 519960 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : /Users/ivanlizarazo/Documents/ivan/GB/notebooks/datos/s2_mosca.tif
## names : s2_mosca.1, s2_mosca.2, s2_mosca.3, s2_mosca.4
## min values : 640, 472, 322, 450
## max values : 5428, 5884, 8544, 10184s2_resc <- s2_crop * 0.0001s2_resc## class : RasterBrick
## dimensions : 737, 751, 553487, 4 (nrow, ncol, ncell, nlayers)
## resolution : 10, 10 (x, y)
## extent : 584280, 591790, 512590, 519960 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : /private/var/folders/l_/1nv5l0dx78g_d3k3_t85lrjr0000gn/T/Rtmpq2Pl6n/raster/r_tmp_2020-10-22_101515_10366_83433.grd
## names : s2_mosca.1, s2_mosca.2, s2_mosca.3, s2_mosca.4
## min values : 0.0640, 0.0472, 0.0322, 0.0450
## max values : 0.5428, 0.5884, 0.8544, 1.0184Let’s calculate the normalized difference vegetation index (NDVI):
s2_ndvi <- (s2_resc$s2_mosca.4 - s2_resc$s2_mosca.3)/(s2_resc$s2_mosca.4 + s2_resc$s2_mosca.3)s2_ndvi ## class : RasterLayer
## dimensions : 737, 751, 553487 (nrow, ncol, ncell)
## resolution : 10, 10 (x, y)
## extent : 584280, 591790, 512590, 519960 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : layer
## values : -0.2767748, 0.8522812 (min, max)Let’s normalize the NDVI values into the 0-255 range:
normalize <- function(x) {
min <- raster::minValue(x)
max <- raster::maxValue(x)
return(255* (x - min) / (max - min))
}s2_ndvi2 <- normalize(s2_ndvi)s2_ndvi2## class : RasterLayer
## dimensions : 737, 751, 553487 (nrow, ncol, ncell)
## resolution : 10, 10 (x, y)
## extent : 584280, 591790, 512590, 519960 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : layer
## values : 0, 255 (min, max)rf <- writeRaster(s2_ndvi2, filename="./datos/ndvi_mosca.tif", format="GTiff", overwrite=TRUE)library(leaflet)
library(RColorBrewer)
pal <- colorNumeric(c("red", "orange", "yellow", "green", "darkgreen"), values(s2_ndvi2),
na.color = "transparent")leaflet() %>% addTiles() %>%
addRasterImage(s2_ndvi2, colors = pal, opacity = 0.8) %>%
addLegend(pal = pal, values = values(s2_ndvi2),
title = "NDVI from Sentinel 2 - 04.01.2020")5. Reproducibility
sessionInfo()## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] RColorBrewer_1.1-2 leaflet_2.0.3.9000 forcats_0.5.0 stringr_1.4.0
## [5] dplyr_1.0.2 purrr_0.3.4 readr_1.3.1 tidyr_1.1.0
## [9] tibble_3.0.4 ggplot2_3.3.2 tidyverse_1.3.0 sf_0.9-6
## [13] raster_3.3-13 gdalUtils_2.0.3.2 rgdal_1.5-18 sp_1.4-4
##
## loaded via a namespace (and not attached):
## [1] httr_1.4.2 jsonlite_1.7.1 foreach_1.5.1 R.utils_2.10.1
## [5] modelr_0.1.6 assertthat_0.2.1 cellranger_1.1.0 yaml_2.2.1
## [9] pillar_1.4.6 backports_1.1.10 lattice_0.20-41 glue_1.4.2
## [13] digest_0.6.26 rvest_0.3.5 colorspace_1.4-1 htmltools_0.5.0
## [17] R.oo_1.24.0 pkgconfig_2.0.3 broom_0.5.6 haven_2.2.0
## [21] scales_1.1.1 generics_0.0.2 farver_2.0.3 ellipsis_0.3.1
## [25] withr_2.3.0 cli_2.1.0 magrittr_1.5 crayon_1.3.4
## [29] readxl_1.3.1 evaluate_0.14 R.methodsS3_1.8.1 fs_1.4.1
## [33] fansi_0.4.1 nlme_3.1-147 xml2_1.3.2 class_7.3-17
## [37] tools_3.6.3 hms_0.5.3 lifecycle_0.2.0 munsell_0.5.0
## [41] reprex_0.3.0 compiler_3.6.3 e1071_1.7-4 rlang_0.4.8
## [45] classInt_0.4-3 units_0.6-7 grid_3.6.3 iterators_1.0.13
## [49] rstudioapi_0.11 htmlwidgets_1.5.2 crosstalk_1.1.0.1 base64enc_0.1-3
## [53] rmarkdown_2.4.6 gtable_0.3.0 codetools_0.2-16 DBI_1.1.0
## [57] R6_2.4.1 lubridate_1.7.9 knitr_1.30 KernSmooth_2.23-17
## [61] stringi_1.5.3 Rcpp_1.0.5 vctrs_0.3.4 png_0.1-7
## [65] dbplyr_1.4.3 tidyselect_1.1.0 xfun_0.18