Customize
Customizing means transforming raw images into useful information for particular needs. Customizing includes; (1) mosaicking, i.e. joining scenes of a region and date in a single file , (2) calculating and index to highlight the presence of process/material in the satellite image, and (3) mask the cloudy pixels to avoid misinterpreting the surface reflectance values. This demo builds on the showcase from the search vignette and so, the first section reviews the most important code from the previous vignette.
Review
As a first step of rsat’s workflow is specifying the
credentials for the the web services:
library(rsat)
set_credentials("rsat.package","UpnaSSG.2021")The showcase aims at assessing the effect of the Snowstorm
Filomena on the Iberian peninsula during January \(10^{th}\) and \(15^{th}\), \(2021\). Hence, the roi and
toi correspond to an sf polygon around the
peninsula (ip) and a vector of dates (toi)
covering the time-span:
ip <- st_sf(st_as_sfc(st_bbox(c(
xmin = -9.755859,
xmax = 4.746094,
ymin = 35.91557,
ymax = 44.02201
), crs = 4326)))
toi <- seq(as.Date("2021-01-10"),as.Date("2021-01-15"),1)The folders for the database and dataset can be created programmatically as follows:
db.path <- file.path(tempdir(),"database")
ds.path <- file.path(tempdir(),"datasets")
dir.create(db.path)
dir.create(ds.path)The minimum information to generate a new rtoi is a
name for the object, a polygon of the region of interest
(roi), and the paths to the database and dataset:
filomena <- new_rtoi(name = "filomena",
region = ip,
db_path = db.path,
rtoi_path = ds.path)To limit the amount of data and processing times, the assessment is conducted over MODIS imagery. A total number of \(24\) images are found for the region over the \(6\)-day period:
rsat_search(region = filomena, product = c("mod09ga"), dates = toi)The way to download the search results is as follows:
rsat_download(filomena)Mosaic
Mosaicking involves binding together
several images of a region from the same date. The function
mosaic() finds automatically the relevant images in the
database and joins them together in a single file. Additionally, by
default, the function crops around the roi of the
rtoi to remove unnecessary information and save space on
your hard disk drive:
rsat_mosaic(filomena)The cropping option can be disabled with the argument
warp = NULL. Here, cropping is appropriate since images
extend far beyond our region of interest.
The results are saved under rtoi_path inside
Modis/mod09ga/mosaic
(i.e. constellation_name/data_product_name/mosaic). This is the
first time in the workflow that the rtoi_path is being
used. The reason is that mosaicking is the first transformation applied
to the raw images to better fit the particular needs of the analysis.
The outcomes from the mosaic are compressed (zip) to minimize
their size:
list.files(file.path(ds.path, "filomena", "Modis/mod09ga/mosaic"), full.name = TRUE)At this point of the workflow, RGB representations of the
satellite images can be displayed on the fly, without loading the entire
image in R. By default, the plot() method for
an rtoi displays the mosaicked images when the object is
not followed by the words "preview" or "date"
(see other types of plots for an rtoi in the search
vignette):
The map shows a sample of pixels from the original image. This is a
strategy to save space in RAM memory. The downside is the
reduction in the level of detail of the image. By default, the number of
pixels on the horizontal and vertical axis is \(250\). The arguments xsize and
ysize change the size of the sample to vary the crispness
of the image. The code below doubles the number of pixels on each axis
of the image;
Clouds and snow are visually similar. False color images are
frequently used in the field in order to ease the detection of features.
False color images switch the natural color in the RGB image to
other bands. These kind of representations can be built using the
band_name argument followed by the name of the bands
replacing the red, green, and blue colors.
For instance, the rendering of swir1, nir, and
blue highlights the snow in blue color:
Index calculation
Definition
A remote sensing index is an indicator that reveals the presence of a material in a satellite image. Indexes are the result of simple math applied to the bands of an image. The computation involves the bands with a distinctively high or low reflectance for the feature at hand. Over the years, researchers have developed a wide variety of indexes for different materials or processes which can be consulted here.
For instance, the Normalized Difference Snow Index (NDSI) (see e.g., (Salomonson and Appel 2004)) highlights the snow using the green and shortwave-infrared bands (around \(1.5 \mu m\)). The subtraction of this two bands gives a large number for those pixels depicting snow. The denominator ensures that values oscillate between \(-1\) and \(1\).
\[ NDSI = \frac{Green - SWIR1}{Green + SWIR1}\]
Calculation
In R we can create a function replicating the
calculation of the NDSI index:
NDSI = function(green, swir1){
ndsi <- (green - swir1)/(green + swir1)
return(ndsi)
}rsat demands that these formulas use the band names,
e.g. red, green, blue, etc. rather than band number. Band names
and numbers differ among mission/satellites. For instance, the
green corresponds to the band number \(4\) in MODIS and Landsat-7, number \(3\) in Landsat-8 and Sentinel-2, and number
\(6\) in Sentinel-3 (see here).
Using their names enables the use of a unique custom function across
satellites/missions. rsat functions are responsible for
linking the name to the corresponding band number according to the
mission. Some widespread variables are built-in the package and the list
of variables can be printed using;
To use the NDSI() function over the series of satellite
images of the Iberian peninsula, use the function derive()
as follows;
rsat_derive(filomena, product = "mod09ga", variable = "ndsi", fun = NDSI)Again, you can plot the results without loading the scenes in
R:
The NDSI index improves the separability between clouds and snow. However, there might be some difficulties distinguishing between them in certain parts of the image. As a solution, the next step removes cloud-covered pixels.
Cloud removal
Some data providers apply algorithms over their data-sets to detect the presence of clouds (Level 1/2 products). The analysis is part of the quality assessment done during pre-processing and the results are included in the Quality Assurance (QA) band of the image. In addition to cloud coverage, the band provides information about over-saturated or filled pixels. The information is packed in this band using the bit format.
The function cloud_mask() interprets the QA
band to obtain images showing the presence/absence of clouds. Its
application is straightforward;
rsat_cloudMask(filomena)To apply the cloud mask, we need to import the NDSI images
into R using the get_raster() function. The
values of the index must be truncated between \(-1\) and \(1\) to avoid values outside the feasible
range (sun reflections on mirror-like surfaces, such as water, can lead
to misleading results):
ndsi.img <- rsat_get_raster(filomena, "mod09ga", "ndsi")
ndsi.img <- clamp(ndsi.img, -1, 1)For every image in the rtoi, the
cloud_mask() function generates a new image, called
mask, in which \(1\)s and \(NA\)s indicate clear and covered pixels.
The function identifies the mission/program and applies the appropriate
interpretation of bits to create the cloud mask. To import the result
run;
clds.msk <- rsat_get_raster(filomena, "mod09ga", "CloudMask")In MODIS, cloud-masks have a different resolution than the
multispectral image. To adjust the resolution, resample the
cloud mask to match the resolution of the NDSI images
(resample()) using the nearest neighbor method
("ngb"):
clds.msk <- resample(clds.msk, ndsi.img, method = "ngb")To apply the cloud mask, we just multiply both series of pixels. Dot multiplications are performed pixel-wise. NDSI values multiplied by \(1\) remain unaltered but those multiplied by \(NA\) become missing:
As an attempt to obtain a composite image, we extract maximum value of the NDSI for each pixel in the time series. Maximum value compositions are frequent in this field (Holben 1986). Compositing is as a way to summarize the information in a time-lapse and ignore the presence of clouds:
snow.spain <- calc(ndsi.filt, max, na.rm = TRUE)Represent the results:
The results are not completely satisfactory. The image shows unfilled gaps and this is because for those pixels there is no valid information along the time-series. The following vignette explains how to process images to fill gaps and smooth outliers.