Profiles

rsat communicates with two public data archives through their Application Programming Interfaces (APIs). The registration in the following online portals is required to get a full access to satellite images with rsat;

• USGS USGS is the sole science agency for the Department of the Interior of United States. Provide access to Modis Images. More information about USGS can be found Here.
• EarthData: A repository of NASA’s earth observation data-sets. More information about EarthData can be found here.
• Dataspace, a web service giving access to Copernicus’ scientific data hub. Please go here to find more details about the data hub.

For convenience, try to use the same username and password for all of them. To satisfy the criteria in both web services make sure that the username is \(4\) characters long and includes a period, number or underscore. The password must be \(12\) character long and should include characters with at least one capital letter, and numbers.

Credentials management

Two functions deal with the usernames and passwords for the various APIs considered in rsat; set_credentials() and print_credentials(). The former defines which username (first argument) and password (second argument) should be used for which portal (third argument):

library(rsat)
set_credentials("rsat.package","UpnaSSG.2021", "scihub")
set_credentials("rsat.package","UpnaSSG.2021", "earthdata")

The latter prints the usernames and passwords saved during the current R session:

print_credentials()

If usernames and passwords are the same for all the portals, they can be set with a single instruction:

set_credentials("rsat.package","UpnaSSG.2021")

Simple search

data("ex.navarre")
toi <- as.Date("2021-01-01") + 0:15
rcd <- rsat_search(region = ex.navarre,
                   product = c("mod09ga"),
                   dates = toi)

rcd <- rsat_search(region = ex.navarre,
                   product = c("mod09ga", "SY_2_SYN___"),
                   dates = toi)
class(rcd)

unique(sat_name(rcd))
names(rcd)[1]
unique(dates(rcd))

Filtering search results

As mentioned earlier, a records object behaves like a vector. It works with common R methods such as c(), [], length(), subset(), or unique(). These methods allow to filter and tinker with the search results. Selecting and filtering satellite images is specially powerful when combined with visualization (see below). Discarding useless images at this stage of the process can save memory space and processing time.

For instance, the code below counts the results found from each mission:

mod.rcd <- subset(rcd, "sat", "Modis")
sn3.rcd <- subset(rcd, "sat", "Sentinel-3")
length(mod.rcd)
length(sn3.rcd)

The first, second, and third argument in subset() correspond to, (1) the records object, (2) the name of the slot used for subsetting, and (3) the value of the slot we are interested in. The next lines of code show a more advance filtering example where a new records is built from images captured by both missions on the same dates:

mod.mtch <- mod.rcd[dates(mod.rcd) %in% dates(sn3.rcd)]
sn3.mtch <- sn3.rcd[dates(sn3.rcd) %in% dates(mod.rcd)]
rcd <- c(mod.mtch, sn3.mtch)

A records can be coerced to a data.frame with as.data.frame() and transformed back with as.records(). The rows and columns of the data.frame correspond to the satellite images and their slots respectively:

rcd.df <- as.data.frame(rcd)
dim(rcd.df)
names(rcd.df)
# rcd <- as.records(rcd.df)
# class(rcd)

Previewing search results

The rsat package provides a plot() method for satellite records. It displays a preview of the satellite images, i.e a low-resolution versions of the satellite images. Previewing can help to decide whether the actual image, much heavier than the preview, is worth downloading.

Cloudy images are frequently useless and they can be removed from the records object using vector-like methods (see previous section). For instance, the code below displays the first \(10\) records and, based on visual evaluation of cloud coverage, it removes the \(9^{th}\) image from the records:

prview <- plot(rcd[1:12])
prview
rcd <- rcd[-9]

Sometimes, it might be difficult to spot the location of the region of interest on the previews. The roi polygon can be passed to the plot() function using the region argument. Additionally, the plot() method in rsat is a wrapper function of tmap and it accepts other inputs from tm_raster() and tm_polygon(). These arguments should be preceded by the tm.raster and tm.polygon identifiers.

The instruction below leverages the preview to its fullest. It shows the RGB satellite image preview with the roi superposed (), with a red border color (tm.polygon.region.border.col = "red") and no fill (tm.polygon.region.alpha = 0). The compass (compass.rm = TRUE) and scale bar (scale.bar.rm = TRUE) are removed:

plot(rcd[1:6],
     region = ex.navarre,
     tm.polygon.region.border.col = "red",
     tm.polygon.region.alpha = 0,
     compass.rm = T,
     scale.bar.rm = T)

The rtoi

Sometimes, regions are located at the intersection of images or they are too large to fit in a single scene. In these situations, working with separate files might be cumbersome. As a solution, rsat provides the rtoi class object. An rtoi is a project or study-case which is associated with a region and time of interest (hence its name). An rtoi consists of the following elements:

1. name: a project identifier.
2. region: a georeferenced polygon (sf class object).
3. db_path: the path to a database, i.e. a folder for raw satellite images which have generic purposes.
4. rtoi_path: the path to a dataset, i.e. a folder for customized/processed.
5. records: a vector of satellite records relevant for the study.

As a showcase, we’ll assess the effects of the Storm Filomena over the Iberian peninsula in terms of snow coverage. The storm was an extreme meteorological event (largest since 1971, according to AEMET). The storm swept the peninsula from January \(6^{th}\) and \(11^{th}\), \(2021\). The code below generates a bounding box around the peninsula (ip) and limits the study period (toi) to the immediate dates after the storm:

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)

A new rtoi requires at least the elements from \(1\) to \(4\). The following lines generate programmatically the database and dataset folders:

db.path <- file.path(tempdir(),"database")
ds.path <- file.path(tempdir(),"datasets")
dir.create(db.path)
dir.create(ds.path)

The function new_rtoi() builds a new rtoi from the elements mentioned above:

filomena <- new_rtoi(name = "filomena",
                     region = ip,
                     db_path = db.path,
                     rtoi_path = ds.path)

The new_rtoi() function writes a set of files representing a copy of the R object (filomena.rtoi) and a shapefile of the roi (region):

rtoi.files <- list.files(file.path(ds.path, "filomena"), full.name = TRUE)
rtoi.files

This file is updated whenever the rtoi is modified. Thus, if the R session accidentally breaks or closes, the last version of the rtoi is saved in you hard drive and can be loaded with read_rtoi():

filomena <- read_rtoi(file.path(ds.path, "filomena"))

Printing the rtoi provides a summary about the region and time of interest, a summary of the relevant records for the analysis, and the paths to the database and the dataset:

print(filomena)

Search with an rtoi

You can use an rtoi to search satellite images with rsat_search():

toi <- as.Date("2021-01-10") + 0:5
rsat_search(region = filomena,
            product = c("mod09ga", "SY_2_SYN___"),
            dates = toi)

An rtoi is an S6 class object, so it behaves like any object in other programming languages. That is, if the rtoi is passed as an argument and modifies the rtoi, the object in the main environment is also updated. The function rsat_search() places the resulting records into the rtoi:

print(filomena)

Since the rtoi has been modified, the rtoi.file has also been updated:

file.info(rtoi.files[1])$ctime # creation time
file.info(rtoi.files[1])$mtime # modification time

To extract the records from the rtoi use records():

rcds <- records(filomena)
class(rcds)

If you want to apply further filters on the results, as shown in the previous section, you need to extract the records from the rtoi first and then apply c(), [], subset(), and length() at your convenience.

Displaying results

There is a plot() method too for rtois but it is more sophisticated. There is a "preview" mode that shows the low-resolution RGB version of the satellite images. The plot binds together the satellite images available for a given date and roi, and the function allows other arguments, such as product or dates, to better select the information being shown:

plot(filomena,
     "preview",
     product = "mod09ga",
     dates = "2021-01-11")

rtois are designed save information from several missions. The mode = "dates" prints a calendar indicating the availability of satellite images from one or multiple missions during the time of interest:

plot(filomena, "dates")

The following vignette explains how to download the satellite images and the role of the database when managing multiple rtois. To reduce the processing times and memory space, the showcase is restricted to MODIS images from here onward:

rcd <- records(filomena)
records(filomena) <- subset(rcd, "product", "mod09ga")