1. Introduction

This notebook illustrates how to obtain time series of remote sensing data using a simple API. This notebook is written to help Geomatica Basica students at Universidad Nacional de Colombia to become familiar with geospatial functionalities provided by the R software environment. This notebook’s code follows closely the Web Time Series Service (WTSS) R package vignette.

A WTSS server takes as input an Earth observation data cube, that has a spatial and a temporal dimension and can be multidimensional in terms of its attributes. The WTSS API has three commands, which are are (a) list_coverages, that returns a list of coverages available in the server; (b) describe_coverage, that that returns the metadata for a given coverage; (c) time_series, that returns a time series for a spatio-temporal location.

WTSS-R considers that “coverages” are equivalent to “data cubes”. Data cubes rely on the fact that Earth observation satellites revisit the same place at regular intervals. Thus, remote sensing data can be calibrated so that observations of the same place in different times are comparable. These calibrated observations can be organised in regular intervals, so that each measure from sensor is mapped into a three dimensional multivariate array in space-time.

This notebook includes the following sections:

2. Connecting to a WTSS server

Let’s start by loading the required libraries:

##devtools::install_github("e-sensing/wtss")
library(rgdal)
## Loading required package: sp
## rgdal: version: 1.4-8, (SVN revision 845)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/gdal
##  GDAL binary built with GEOS: FALSE 
##  Loaded PROJ.4 runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/proj
##  Linking to sp version: 1.3-2
library(raster)
library(sf)
## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.1     ✓ dplyr   1.0.0
## ✓ 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(tmap)
library(gstat)
library(sp)
library(wtss)
## wtss - R interface to Web Time Series Service.
## Loaded wtss v2.1.0.
##         See ?wtss for help, citation("wtss") for use in publication.
##         See demo(package = "wtss") for examples.
#install.packages("bfast")
library(bfast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
#remotes::install_github("christophsax/tsbox")
library(tsbox)
#install.packages("seasonal")
library(seasonal)
## 
## Attaching package: 'seasonal'
## The following object is masked from 'package:tibble':
## 
##     view
library(tsibble)
## 
## Attaching package: 'tsibble'
## The following object is masked from 'package:raster':
## 
##     stretch

The first step towards using the service is connecting to a server that supports the WTSS protocol. Currenlty, Brazil’s National Insitute for Space Research (INPE) runs such a service. In the package, the connection is enabled by using the URL of the service. The package informs if the connection has been correctly made.

# Connect to the WTSS server at INPE Brazil
wtss_inpe <-  wtss::WTSS("http://www.esensing.dpi.inpe.br/wtss")
## Connected to WTSS server at http://www.esensing.dpi.inpe.br/wtss

3. Listing coverages available at the WTSS server

This operation allows clients to retrieve the capabilities provided by any server that implements WTSS. It returns a list of coverage names available in a server instance.

# Connect to the WTSS server at INPE Brazil
wtss::list_coverages(wtss_inpe)
## Coverages: MOD13Q1 MOD13Q1_M

The Terra Moderate Resolution Imaging Spectroradiometer (MODIS) Vegetation Indices (MOD13Q1) Version 6 data are generated every 16 days at 250 meter (m) spatial resolution as a Level 3 product.

The MOD13Q1 product provides two primary vegetation layers. The first is the Normalized Difference Vegetation Index (NDVI) which is referred to as the continuity index to the existing National Oceanic and Atmospheric Administration-Advanced Very High Resolution Radiometer (NOAA-AVHRR) derived NDVI. The second vegetation layer is the Enhanced Vegetation Index (EVI), which has improved sensitivity over high biomass regions. The algorithm chooses the best available pixel value from all the acquisitions from the 16 day period. The criteria used is low clouds, low view angle, and the highest NDVI/EVI value.

More information on these products can be found here.

4. Describing a coverage from the WTSS server

This operation returns the metadata for a given coverage identified by its name. It includes its range in the spatial and temporal dimensions.

# Connect to the WTSS server at INPE Brazil
metadata <- wtss::describe_coverage(wtss_inpe, name = "MOD13Q1")
## Coverage description saved in WTSS object
str(metadata)
## List of 4
##  $ url        : chr "http://www.esensing.dpi.inpe.br/wtss"
##  $ coverages  : chr [1:2] "MOD13Q1" "MOD13Q1_M"
##  $ description: list()
##  $ valid      : logi TRUE
##  - attr(*, "class")= chr [1:2] "wtss" "list"

5. Obtaining a time series

This operation requests the time series of values of a coverage attribute at a given location. Its parameters are: (a) wtss.obj: either a WTSS object (created by the operation wtss::WTSS as shown above) or a valid WTSS server URL; (b) name: Cube (coverage) name; (c) attributes: vector of band names (optional). If omitted, all bands are retrieved; (d) longitude: longitude in WGS84 coordinate system; (e)latitude: Latitude in WGS84 coordinate system; (f)start_date (optional): Start date in the format yyyy-mm-dd or yyyy-mm depending on the coverage. If omitted, the first date on the timeline is used; (g) end_date(optional): End date in the format yyyy-mm-dd or yyyy-mm depending on the coverage. If omitted, the last date of the timeline is used.

# Request a time series from the "MOD13Q1" coverage
(evi_ts <- wtss::time_series(wtss_inpe, "MOD13Q1", attributes = c("evi"), latitude = 1.59479, longitude = -74.3439))
## # A tibble: 1 x 7
##   longitude latitude start_date end_date   label   cube    time_series       
##       <dbl>    <dbl> <date>     <date>     <chr>   <chr>   <list>            
## 1     -74.3     1.59 2000-02-18 2019-09-30 NoClass MOD13Q1 <tibble [423 × 2]>

The result of the operation is a tibble which contains data and metadata. The first six columns contain the metadata: satellite, sensor, spatial and temporal information, and the coverage from where the data has been extracted. The spatial location is given in longitude and latitude coordinates for the “WGS84” ellipsoid. The time_series column contains the time series data for each spatiotemporal location. This data is also organized as a tibble, with a column with the dates and the other columns with the values for each spectral band.

# Showing the contents of a time series
evi_ts$time_series[[1]]
## # A tibble: 423 x 2
##    Index        evi
##    <date>     <dbl>
##  1 2000-02-18 0.231
##  2 2000-03-05 0.562
##  3 2000-03-21 0.549
##  4 2000-04-06 0.320
##  5 2000-04-22 0.495
##  6 2000-05-08 0.395
##  7 2000-05-24 0.185
##  8 2000-06-09 0.587
##  9 2000-06-25 0.483
## 10 2000-07-11 0.195
## # … with 413 more rows

For convenience, the WTSS package provides a convenience function for plotting the time series.

# Plotting the contents of a time series
plot(evi_ts)

6. Detecting breakpoints

Since many time series analysis functions in R require data to be made available in the “zoo” and “ts”" formats, the wtss package provides two convenience functions: \(wtss_to_zoo\) and \(wtss_to_ts\). The example below shows the detection of breakpints in a series converted to the “ts”" format using the \(bfast\) package written by Verbesselt.

# convert to ts
ts1 <- wtss::wtss_to_ts(evi_ts, band = "evi")
head(time(ts1))
## Time Series:
## Start = c(2000, 7) 
## End = c(2000, 12) 
## Frequency = 52 
## [1] 2000.115 2000.135 2000.154 2000.173 2000.192 2000.212

Before proceeding, it is important to check the range of EVI values.

summary(ts1)
##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -2503.5891     0.4043     0.4909   -40.6767     0.5606     0.7717

Notice that there are several values out of the expected range of the vegetation index. A quick fix is to select a period of time in which values are valid:

ts2 <- window(ts1, start = 2001, end = 2018)
summary(ts2)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1074  0.4310  0.5048  0.4917  0.5718  0.7717
plot(ts2, ylab="EVI")

Let’s convert the time series object to a tibble object:

newts2 <- ts_tsibble(ts2)
head(newts2)
## # A tsibble: 6 x 2 [34m 12s 488ms 688.00001218915µs] <?>
##   time                value
##   <dttm>              <dbl>
## 1 2001-01-01 00:00:00 0.536
## 2 2001-01-08 00:34:12 0.543
## 3 2001-01-15 01:08:24 0.526
## 4 2001-01-22 01:42:37 0.487
## 5 2001-01-29 02:16:49 0.455
## 6 2001-02-05 02:51:02 0.458
tail(newts2)
## # A tsibble: 6 x 2 [34m 12s 488ms 688.00001218915µs] <?>
##   time                value
##   <dttm>              <dbl>
## 1 2017-11-26 21:08:57 0.418
## 2 2017-12-03 21:43:10 0.320
## 3 2017-12-10 22:17:22 0.349
## 4 2017-12-17 22:51:35 0.442
## 5 2017-12-24 23:25:47 0.501
## 6 2018-01-01 00:00:00 0.503
class(newts2)
## [1] "tbl_ts"     "tbl_df"     "tbl"        "data.frame"
ts_plot(`Raw series` = newts2,
`Adjusted series` = ts_trend(newts2),
title = "EVI Time Series from MODIS at Lat = 1.59, Long = -74.34",
subtitle = "Temporal trend"
)

Now, let’s try bfast break detection functionality:

fit <- bfast01(ts2)
## 
## Attaching package: 'stlplus'
## The following object is masked from 'package:seasonal':
## 
##     trend
plot(fit, main="Break detected in EVI Time Series at Lat = 1.59, Long = -74.34")

This suggest an abrupt change in vegetation aprox. on 2014. Is this breakpoint a true signal or a false one?

7. Fact-checking

Well, let’s fact-check. In Leonardo Hurtado’s Master Thesis the following figure is included. It shows that, in this location, forest was cleared at the end of 2015.

library(knitr)
include_graphics("../pngs/reference_degrad1.png")

Is it a lucky shot? Well, let’s try the function in another location:

# Request a time series from the "MOD13Q1" coverage
(evi_ts <- wtss::time_series(wtss_inpe, "MOD13Q1", attributes = c("evi"), latitude = 1.021, longitude = -74.7))
## # A tibble: 1 x 7
##   longitude latitude start_date end_date   label   cube    time_series       
##       <dbl>    <dbl> <date>     <date>     <chr>   <chr>   <list>            
## 1     -74.7     1.02 2000-02-18 2019-09-30 NoClass MOD13Q1 <tibble [423 × 2]>
# Plotting the contents of a time series
plot(evi_ts)

# convert to ts
ts1 <- wtss::wtss_to_ts(evi_ts, band = "evi")
ts2 <- window(ts1,  end = 2018)
fit <- bfast01(ts2)
plot(fit)

What did Leonardo say?

library(knitr)
include_graphics("../pngs/reference_degrad2.png")

That’s all. This was a good web service experience.

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] knitr_1.28       stlplus_0.5.1    tsibble_0.9.1    seasonal_1.7.1  
##  [5] tsbox_0.2.1.9000 bfast_1.6.0.9000 wtss_2.1.0       gstat_2.0-6     
##  [9] tmap_3.0         forcats_0.5.0    stringr_1.4.0    dplyr_1.0.0     
## [13] purrr_0.3.4      readr_1.3.1      tidyr_1.1.0      tibble_3.0.1    
## [17] ggplot2_3.3.2    tidyverse_1.3.0  sf_0.9-3         raster_3.1-5    
## [21] rgdal_1.4-8      sp_1.4-2        
## 
## loaded via a namespace (and not attached):
##   [1] x13binary_1.1.39-2 leafem_0.1.1       colorspace_1.4-1  
##   [4] ellipsis_0.3.1     class_7.3-17       leaflet_2.0.3     
##   [7] base64enc_0.1-3    fs_1.4.1           dichromat_2.0-0   
##  [10] yaImpute_1.0-32    rstudioapi_0.11    farver_2.0.3      
##  [13] fansi_0.4.1        lubridate_1.7.9    xml2_1.3.2        
##  [16] codetools_0.2-16   jsonlite_1.6.1     tmaptools_3.0     
##  [19] broom_0.5.6        anytime_0.3.7      dbplyr_1.4.3      
##  [22] png_0.1-7          compiler_3.6.3     httr_1.4.1        
##  [25] backports_1.1.8    assertthat_0.2.1   strucchange_1.5-1 
##  [28] cli_2.0.2          htmltools_0.4.0    tools_3.6.3       
##  [31] gtable_0.3.0       glue_1.4.1         reshape2_1.4.4    
##  [34] Rcpp_1.0.4.6       cellranger_1.1.0   fracdiff_1.5-1    
##  [37] vctrs_0.3.1        urca_1.3-0         nlme_3.1-147      
##  [40] leafsync_0.1.0     crosstalk_1.1.0.1  lmtest_0.9-37     
##  [43] timeDate_3043.102  lwgeom_0.2-3       xfun_0.14         
##  [46] rvest_0.3.5        lifecycle_0.2.0    XML_3.99-0.3      
##  [49] zoo_1.8-8          scales_1.1.1       hms_0.5.3         
##  [52] parallel_3.6.3     sandwich_2.5-1     RColorBrewer_1.1-2
##  [55] yaml_2.2.1         quantmod_0.4.17    curl_4.3          
##  [58] geosphere_1.5-10   stringi_1.4.6      highr_0.8         
##  [61] tseries_0.10-47    e1071_1.7-3        TTR_0.23-6        
##  [64] bitops_1.0-6       intervals_0.15.2   rlang_0.4.6       
##  [67] pkgconfig_2.0.3    evaluate_0.14      lattice_0.20-41   
##  [70] labeling_0.3       htmlwidgets_1.5.1  tidyselect_1.1.0  
##  [73] plyr_1.8.6         magrittr_1.5       R6_2.4.1          
##  [76] generics_0.0.2     DBI_1.1.0          pillar_1.4.4      
##  [79] haven_2.2.0        withr_2.2.0        units_0.6-6       
##  [82] stars_0.4-2        xts_0.12-0         RCurl_1.98-1.2    
##  [85] abind_1.4-5        nnet_7.3-14        spacetime_1.2-3   
##  [88] modelr_0.1.6       crayon_1.3.4       utf8_1.1.4        
##  [91] KernSmooth_2.23-17 rmarkdown_2.1      grid_3.6.3        
##  [94] readxl_1.3.1       data.table_1.12.8  FNN_1.1.3         
##  [97] reprex_0.3.0       forecast_8.12      digest_0.6.25     
## [100] classInt_0.4-3     munsell_0.5.0      viridisLite_0.3.0 
## [103] quadprog_1.5-8