ReLTER_demo.Rmd

---
title: "Demonstrating the Use of the `ReLTER` package for Earth Observation (EO)"
author: "Micha Silver and Alessandro Oggioni"
date: "1/02/2022"
output:
  github_document:
    toc: yes
  pdf_document:
    toc: yes
  html_document: default
---

This code demonstrates the use of the new `ReLTER` package. For more details, see 
Allesando Oggioni's [github page](https://github.com/oggioniale/ReLTER)

Alessandro Oggioni, Micha Silver, Luigi Ranghetti & Paolo Tagliolato. (2021). oggioniale/ReLTER: ReLTER v1.0.0 (1.0.0). Zenodo. https://doi.org/10.5281/zenodo.5576813

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(tidy = TRUE)
```

# Install and load packages

Begin by installing packages and loading them.

```{r packages, message = FALSE, results='hide'}
# These packages are required
pkg_list  <-  c("remotes",
             "tmap", "tmaptools",
             "sf", "terra",
             "OpenStreetMap",
             "raster"
            )

# Check if already installed, install if not
installed_packages <- pkg_list %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
  install.packages(pkg_list[!installed_packages])
}
# Load Packages
lapply(pkg_list,
       function(p) {require(p,
                            character.only = TRUE,
                            quietly=TRUE)})

# Now install `ReLTER` from github and load
#remotes::install_github("oggioniale/ReLTER")
# If you want the development version, with latest functions:
remotes::install_github("oggioniale/ReLTER@dev")
library(ReLTER)

# Choose where to save outputs
Output_dir = "./Output"
if (!dir.exists(Output_dir)) {
  dir.create(Output_dir)
}
```

# Query DEIMS SDR

The `ReLTER` package interfaces with the [DEIMS SDR](https://deims.org/) database.

We can query the database in many ways: 
First retrieve the full URL to a few eLTER sites.
Sites can be selected by *country name*, *site name* or both.
Note that partial matching is also supported. So `country_name = "Austri"`
will find sites in Austria, but not Australia.

```{r deims}
eisen <- get_ilter_generalinfo(country="Austria",
                              site_name = "LTSER Platform Eisen")
eisen_deimsid <- eisen$uri

# Using abbreviated "United K" to differentiate from United States
cairngorms <- get_ilter_generalinfo(country = "United K",
                                   site_name = "Cairngorms National")
cairngorms_deimsid <- cairngorms$uri
```

### General metadata

In `ReLTER` there are functions to grab metadata for the sites. Metadata is
available for a few categories.

These categories are available:

  - 'Affiliations'
  - 'Boundaries'    (spatial layer)
  - 'Contacts'
  - 'EnvCharacts'   (environmental characteristics)
  - 'General',
  - 'Infrastructure'
  - 'Parameters'    (which parameters are collected)
  - 'RelateRes'     (related research)
  - 'ResearchTop'   (research topics)

Here are a few basic examples:

```{r general-info}
response <- get_site_info(eisen_deimsid, category = "ResearchTop")
response$researchTopics

response <- get_site_info(cairngorms_deimsid, category = "Affiliations")
response$affiliation.projects
```
# Spatial queries

Now use the DEIMS ID, acquired above, to get the boundary of a site, 
by setting `category = "Boundaries"`.

### Get boundary of site

```{r boundary}
# Acquire boundary for site
eisen_boundary <- get_site_info(eisen_deimsid, "Boundaries")

# Prepare OSM background tile and plot
osm <- read_osm(eisen_boundary, ext = 1.2)
tmap_mode("plot")

# For interactive maps use:
# tmap_mode("view")
# Then these basemaps are available:
# tm_basemap("Stamen.TerrainBackground") +
# tm_basemap("OpenStreetMap") +

tm_shape(osm) +
	tm_rgb() + 
tm_shape(eisen_boundary) +
  tm_polygons(col = "skyblue", alpha = 0.25, border.col = "blue")
```

### Save boundary as shapefile/geopackage for later

The code below saves the boundary to a file for use in other GIS software.

```{r st_write, message=FALSE}
# Edit here to choose your output directory 
boundary_file <- file.path(Output_dir, "eisen_boundary.gpkg")

# Remove country column since it is a list
# (Some sites extend across country boundaries)
eisen_boundary <- subset(eisen_boundary, select = -country)
st_write(eisen_boundary, dsn = boundary_file, append = FALSE)
```

# Dependency on quality of data in DEIMS SDR

`ReLTER` relies on the DEIMS SDR database for all site queries. Therefore, 
any errors or missing data will obviously be echoed in the `ReLTER` results.
These errors include:

  - Missing information
  - Duplicate names
  - Missing boundary shapefile
  
Here are a few examples:

```{r deims-errors, collapse = TRUE}
eisen_contact <- get_site_info(eisen_deimsid, "Contact")
names(eisen_contact)
# No contact information :-(

kiskun <- get_ilter_generalinfo(country_name = "Hungary",
                                site_name = "KISKUN LTER")
kiskun_deimsid <- kiskun$uri
length(kiskun_deimsid)
# Multiple sites with similar name :-(
# Which to choose? View the list...
kiskun$title
kiskun_deimsid <- kiskun$uri[5]
length(kiskun_deimsid)
kiskun_boundary <- get_site_info(kiskun_deimsid, "Boundaries")

# Oops, no boundary for this site!
```
# Datasets from Open DataScience Europe

The recent efforts by the [Geoharmonizer](https://opendatascience.eu/geoharmonizer-project/) program have resulted in a consolidated set of freely available raster data (gridded) datasets.
All rasters are formatted as Cloud Optimzed Geotiff (COG). These can be viewed on the
[ODS web portal](https://maps.opendatascience.eu/).

The code below demonstrates how to access various data from ODS from within *R*,
and to clip to site boundaries. The datasets currently implemented are:

  - "landcover"     (Landcover at 30 meter resolution from Landsat)
  - "clc2018"       (Corine landcover from 2018)
  - "osm_buildings" (Open Street Maps buildings)
  - "natura2000"
  - "ndvi_spring"
  - "ndvi_summer"
  - "ndvi_autumn"
  - "ndvi_winter"

### High resolution landcover

```{r ods-landcover}
# Use boundary and OSM tile from above
eisen_landcover <- get_site_ODS(eisen_deimsid, "landcover")

tm_shape(osm) +
	tm_rgb() + 
  tm_shape(eisen_landcover) +
  tm_raster(style = "pretty", palette = "RdYlBu", alpha=0.75)
```

### Compare with Corine (lower resolution) Landcover

```{r ods-corine}
eisen_corine <- get_site_ODS(eisen_deimsid, "clc2018")

tm_shape(osm) +
	tm_rgb() + 
  tm_shape(eisen_corine) +
  tm_raster(style = "pretty", palette = "Spectral", alpha=0.75)
```

### NDVI during the spring

```{r ods-ndvi}
# Takes several minutes
eisen_ndvi <- get_site_ODS(eisen_deimsid, "ndvi_spring")
tm_shape(osm) +
	tm_rgb() + 
  tm_shape(eisen_ndvi) +
  tm_raster(style = "pretty", palette = "RdYlGn", alpha=0.75)
```

### Small eLTER sites

ODS data layers are at 30 meter resolution, suitable for small sites.
This code examines the Tereno site at Harsleben.

```{r ods-tereno}
# Acquire Tereno DEIMS ID and boundary
tereno <- get_ilter_generalinfo(country_name = "Germany",
                                 site_name = "Tereno - Harsleben")
tereno_deimsid <- tereno$uri
tereno_boundary <- get_site_info(tereno_deimsid, "Boundaries")

# Prepare new OSM background and plot
osm <- read_osm(tereno_boundary, ext = 1.2)
tereno_ndvi <- get_site_ODS(tereno_deimsid, "ndvi_autumn")
tm_shape(osm) +
	tm_rgb() + 
  tm_shape(tereno_ndvi) +
  tm_raster(style = "pretty", palette = "RdYlGn", alpha=0.75)
```

Again compare Landsat based landcover (30 m.) with Corine 2018 (100 m.)

```{r ods-tereno-landcover}
tereno_landcover <- get_site_ODS(tereno_deimsid, "landcover")
tm_shape(osm) +
	tm_rgb() + 
  tm_shape(tereno_landcover) +
  tm_raster(style = "pretty", palette = "Spectral", alpha=0.75)
tereno_corine <- get_site_ODS(tereno_deimsid, "clc2018")
tm_shape(osm) +
	tm_rgb() + 
  tm_shape(tereno_corine) +
  tm_raster(style = "pretty", palette = "Spectral", alpha=0.75)
```


### Save to a Geotiff file for use in other GIS software

```{r ods-save, echo=TRUE, eval=FALSE, message=FALSE}
# Edit here to choose your output directory 
landcover_file <- file.path(Output_dir, "tereno_landcover.tif")
writeRaster(tereno_landcover, landcover_file, overwrite=TRUE)
```


### Copernicus building area in Saldur river catchment site

This code chunk shows how to acquire OpenStreetMap building footprints
from the ODS datasets. OSM building footprints are available by setting
the `dataset` option to "osm_buildings". Here is an example for the 
Saldur river basin in Italy.

```{r ods-saldur-osm_buildings}
saldur <- get_ilter_generalinfo(country="Italy",
                              site_name = "Saldur River Catchment")
saldur_deimsid <- saldur$uri
saldurRiver_osmLandUse <- get_site_ODS(
  deimsid = saldur_deimsid, dataset = "osm_buildings"
)
saldur_boundary <- get_site_info(
  deimsid = saldur_deimsid, "Boundaries"
)
# Prepare OSM background tile and plot
osm <- read_osm(saldur_boundary, ext = 1.2)
# Hillshade of Saldur river bounding box
saldur_hs <- raster::raster("ReLTER_demo_files/saldur_hillshade.tif")

tm_shape(osm) +
  tm_rgb() +
  tm_compass(type = "arrow", position = c("right", "bottom"), text.size = 1) +
  tm_scale_bar(position = c(0.6, "bottom"), text.size = .8) +
  tm_credits("Data from lcv building Copernicus", position = c("left", "top"), size = 1) +
  tm_layout(legend.position = c("left", "bottom")) +
  tm_shape(saldur_hs) +
  tm_raster(palette = "-Greys", style = "cont", legend.show = FALSE, alpha = .4) +
  tm_shape(saldur_boundary) +
  tm_polygons(col = "skyblue", alpha = 0.2, border.col = "gray") +
  tm_shape(saldurRiver_osmLandUse) +
  tm_raster(style = "cont") +
  tm_layout(legend.outside = TRUE)
```

# Datasets from MODIS

This capability of `ReLTER` relies on another package `MODIStsp`
  L. Busetto, L. Ranghetti (2016) 
  *MODIStsp: An R package for automatic preprocessing of MODIS Land Products time series*, 
  Computers & Geosciences, Volume 97, Pages 40-48, ISSN 0098-3004,
  https://doi.org/10.1016/j.cageo.2016.08.020. URL
  https://github.com/ropensci/MODIStsp/.

### What Products are available? What bands?

```{r modis-products}
get_site_MODIS(show_products = TRUE)
get_site_MODIS(show_bands = "Vegetation_Indexes_Monthly_1Km (M*D13A3)" )
get_site_MODIS(show_bands = "LST_3band_emissivity_8day_1km (M*D21A2)" )
```

### Registering on EarthData website

To acquire MODIS data you must be registered at:
https://urs.earthdata.nasa.gov/
Then use your earthdata username and password in the code below

```{r earthdata}
# Username and password saved in advance
# For example:

# creds <- list("username"="homer", password="Simpsons")
# saveRDS(creds, "earthdata_credentials.rds)

# Then...
creds <- readRDS("~/work/EU_Projects/EO_demo/earthdata_credentials.rds")
```

### Download a time series of MODIS NDVI for Eisenwurzen

```{r modis-ndvi, message=FALSE, fig.width=6}
# Set which product and which bands to get
product = "Vegetation_Indexes_Monthly_1Km (M*D13A3)"
bands = "NDVI"    # Also EVI is available

eisen_ndvi <- get_site_MODIS(eisen_deimsid,               # Which site
                          earthdata_user = creds$username, # From above
                          earthdata_passwd = creds$password,
                          product = product,
                          bands = bands,
                          from_date = "2020.01.01", 
                          to_date =  "2020.12.31",        # From-To dates
                          scale = TRUE,          # Rescale COG back to real values
                          out_folder = Output_dir,        # Where to save outputs
                          save_ts_dir = Output_dir
                          )

# Plot was saved to:
plot_file <- paste("time_series", paste(bands, collapse="_"), sep="_")
plot_path <- file.path(Output_dir, paste0(plot_file, ".png"))
knitr::include_graphics(plot_path)
```