PhiWeek_RClient_2021.Rmd

---
title: "PhiWeek_2021: Introduction to openEO platform via the R-Client"
author: 
 - Peter Zellner, Eurac Research
 - Florian Lahn, EFTAS
 - Matthias Mohr, University of Münster
date: "15/10/2021"
output:
  html_document:
    toc: true
    toc_depth: 3
    toc_float: TRUE
---

![](https://webassets.eurac.edu/31538/1629894586-open-eo-platform.png?fit=max&fm=jpg&w=1000){width="400"}

![](https://www.eurac.edu/_next/image?url=https%253A%252F%252Fwebassets.eurac.edu%252F31538%252F1634027841-logoproject.png%253Fauto%253Dformat%2526fit%253Dclip%2526fm%253Djpg%2526h%253D300%2526w%253D300&w=640&q=85){width="116"}

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

# Setup

## Before you start...

-   the openEO getting started guide <https://openeo.org/documentation/1.0/r/>
-   the openEO platform getting started guide <https://docs.openeo.cloud/getting-started/r/>
-   and the R-Client repo for logging issues and collaborating <https://github.com/Open-EO/openeo-r-client>

## Topics

-   Installation
-   Authentication and Log-In
-   Exploring the Backend (widgets, using metadata)
-   Example Process - Usability R (tidyverse, geospatial)
-   Outlook User Defined Functions R-UDF

## Install the openEO R-Client
The openEO R-Client is now available on CRAN. The latest stable version can be installed via `install.packages("openeo")`.
In case you want to use a specific version of the R-Client you can install it from the [openEO github repository](https://github.com/Open-EO/openeo-r-client). It can be installed via `remotes::install_github()`. For the dev verstion set `remotes::install_gitub(remotes::install_github(repo="Open-EO/openeo-r-client",ref="develop", dependencies=TRUE)`.

```{r install, eval=FALSE}
install.packages("openeo")
packageVersion("openeo") >= '1.0.1'
```

## Useful libraries

Load the openEO R-Client library and some other useful libraries.

```{r libs, message=FALSE, warning=FALSE}
library(openeo)
library(stars)
library(sf)
library(mapview)
library(mapedit)
library(dplyr)
library(ggplot2)
library(plotly)
```

# Login

First you can connect to the backend to explore the available collections (as seen in the Connections Pane) and the processes.

```{r connect}
host = "https://openeo.cloud"
con = openeo::connect(host)
```

In the code snippets above we did not need to log in since we just queried publicly available back-end information. However, to run non-trivial processing queries one has to authenticate so that permissions, resource usage, etc. can be managed properly.

To handle authentication, openEO leverages [OpenID Connect (OIDC)](https://openid.net/connect/). It offers some interesting features (e.g. a user can securely reuse an existing account), but is a fairly complex topic, discussed in more depth in the [general authentication documentation for openEO Platform](https://docs.openeo.cloud/authentication/#authentication).

Unfortunately, you need to request a Client ID and a Client Secret for this from the openEO Platform support due to the R client not being officially supported by openEO Platform! Once you have received the Client ID and a Client Secret, you can can continue with the instructions below.

```{r login_oidc, eval=FALSE}
creds_rclient = read.csv2("/home/pzellner@eurac.edu/pwd/openeo_platform_rclient.pwd")
client_id = creds_rclient$client_id # "request at openeo platform support" 
secret = creds_rclient$secret # "request at openeo platform support" 
openeo::login(login_type = "oidc", 
              provider = "egi", 
              config = list(client_id = client_id, secret = secret))
#describe_account()
```

# Discover the backend

Once your connected you can discover processes and collections. This gives you valuable metadata on collections and information on how to use the processes.

## Discover collections

The collections appear in the Connections pane of R-Studio. You can also list them via `list_collections()` and retrieve some first info.

```{r list_collections}
names(openeo::list_collections()) %>% as.data.frame()
```

To get more details you can use the `collection_viewer()` which opens the description of the collection in the Viewer pane

```{r collection_viewer, eval=FALSE}
collection_viewer("SENTINEL2_L1C_SENTINELHUB")
```

To use the collection metadata for further coding you can use `decribe_collection()`.

```{r describe_collection}
coll_meta = openeo::describe_collection("SENTINEL2_L1C_SENTINELHUB")
ext = coll_meta$extent$spatial # you could use this to draw a bbox of the collection
st_bbox(obj = c(xmin = ext[1], ymin = ext[2], xmax = ext[3], ymax = ext[4]), crs = 4326) %>% mapview()
```

## Discover processes

To get an overview of the available processes use `list_processes()`

```{r list_processes}
names(list_processes()) %>% as.data.frame()
```

To get detailed information about a process you can again use the `process_viewer()`

```{r process_viewer, eval=FALSE}
process_viewer("atmospheric_correction")
```

Or `describe_process()`

```{r describe_process}
describe_process("atmospheric_correction")
```

# Calculate NDVI: From Sentinel 2 data atmospherically corrected on the fly

-   Load Sentinel 2 L1C data
-   Define the area of interest, temporal range and bands
-   Atmospherically correct the values on the fly
-   Calculate the NDVI
-   Download the result

## Definitions

First define an area of interest interactively (either a point or area). You can also use multiple features and loop through them later on.

```{r aoi_process}
# this is interactive
pnts = mapedit::drawFeatures()

# this is static for areas - Frascati
# pnts = c(xmin = 12.64805, ymin = 41.80398, xmax = 12.6517, ymax = 41.80603)
# pnts = st_bbox(pnts, crs = st_crs(4326))
mapview(pnts)

bbox = st_bbox(pnts)
bbox = list(west = bbox[[1]],
            east = bbox[[3]],
            south = bbox[[2]],
            north = bbox[[4]])
```

Define the collection (data set), time range and bands.

```{r defs_process}
collection = "SENTINEL2_L1C_SENTINELHUB"
time_range = list("2018-01-01T00:00:00.000Z", 
                  "2018-12-31T00:00:00.000Z")
bands = c("B04", "B08", "CLP", "B09", "B8A", "B11",
          "sunAzimuthAngles", "sunZenithAngles", "viewAzimuthMean", "viewZenithMean")
#bands = c("B04", "B08")
```

## Build the processing chain

Start building the processing chain. First, load the available processes.

```{r load_processes, eval=FALSE}
p = processes()
# names(p)
```

Load the collection.

```{r load_data, eval=FALSE}
data = p$load_collection(id = collection, 
                         spatial_extent = bbox,
                         temporal_extent = time_range, 
                         bands = bands) %>% 
  p$atmospheric_correction(method = "smac")
  
```

Atmospherical Correction. The different processes can be chained using the pipe `%>%`.

```{r atmopspherical_correction, eval=FALSE}
boa_corr = p$atmospheric_correction(data = data, method = "smac") # or use "iCor"
```

Calculate NDVI. Reduce the "bands" dimension with the well known NDVI formula.

```{r calc_ndvi, eval=FALSE}
ndvi_calc = p$reduce_dimension(data = boa_corr, 
                               dimension = "bands", 
                               reducer = function(data, context) {
                                 red = data[1]
                                 nir = data[2]
                                 (nir-red)/(nir+red)})

```

Aggregate to temporal periods.

```{r temporal_period, eval=FALSE}
# process_viewer("aggregate_temporal")

intervals = list(c('2018-01-02', '2018-02-01'),
                 c('2018-02-01', '2018-03-01'),
                 c('2018-03-01', '2018-04-01'),
                 c('2018-04-01', '2018-05-01'),
                 c('2018-05-01', '2018-06-01'), 
                 c('2018-06-01', '2018-07-01'), 
                 c('2018-07-01', '2018-08-01'),
                 c('2018-08-01', '2018-09-01'),
                 c('2018-09-01', '2018-10-01'), 
                 c('2018-10-01', '2018-11-01'),
                 c('2018-11-01', '2018-12-30'))

labels = lapply(intervals, function(x){x[[1]]}) %>% unlist()

ndvi_mnth = p$aggregate_temporal(data = ndvi_calc,
                                 intervals = intervals,
                                 reducer = function(data, context){p$median(data)},
                                 labels = labels,
                                 dimension = "t")
```

Save the result.

```{r save_result, eval=FALSE}
result = p$save_result(data = ndvi_mnth, format="NetCDF")
```

No processing has happened so far. Only the workflow/process graph has been defined.

```{r vis_graph, eval=FALSE}
graph_info = create_user_process(result, id = "test", submit = FALSE)
parse_graph(graph_info)
```

## Compute the result

synchronous call (result is computed directly)

```{r out_name}
out_name =  "/home/pzellner@eurac.edu/git_projects/sample-notebooks/r_ndvi_mnth_int.nc"
```

```{r compute_result, eval=FALSE}
a = Sys.time()
compute_result(result,
               format = "netCDF",
               output_file = out_name, 
               con = con)
b = Sys.time()-a
b
```

## Load the result

```{r load_result}
# load the data into r
ndvi = read_ncdf(out_name)

# check which projection your data has and assign it
#system(paste0("gdalinfo ", out_name))
st_crs(ndvi) = st_crs(32633)

# look at the time dimension
stars::st_get_dimension_values(ndvi, "t")
```

plot some time slices (select and deselect in the viewer)

```{r plot_area}
brks = seq(-0, 1, 0.1)
mapview(ndvi %>% slice("t", 3), at = brks) + 
  mapview(ndvi %>% slice("t", 5), at = brks) + 
  mapview(ndvi %>% slice("t", 7), at = brks) +
  mapview(ndvi %>% slice("t", 9), at = brks)
```

get a point for plotting a pixel timeseries

```{r get_pixel}
# define a point
pixel = mapedit::drawFeatures()

# static assignment
# pixel = data.frame(x = 12.64956, y = 41.80526)
# pixel = st_as_sf(pixel, coords = c("x", "y"), crs = st_crs(4326))

mapview(pixel) + mapview(st_bbox(ndvi))
```

subset to a point and plot a pixel time series

```{r plot_}
# subset to point
pixel = st_transform(x = pixel, crs = st_crs(ndvi))
ndvi_ts = ndvi[pixel]

# generate a data frame from the values and dates
ndvi_ts_df = data.frame(value = ndvi_ts %>% pull() %>% c() %>% as.vector(), 
                        dates = as.Date(st_get_dimension_values(ndvi_ts, "t")))

# plot the timeseries
plot_ts = ggplot(ndvi_ts_df, aes(x = dates, y = value)) + 
  geom_line() + 
  geom_point()
plot_ts_plotly = plotly::ggplotly(plot_ts)
plot_ts_plotly

```

# UDF Outlook

User defined functions (UDF) allow to run arbitrary R-Code on openEO platform.

Main Concepts:

-   Easy data access for local prototyping or training
-   Refine R-Code to produce desired results
-   Run R-Code in openEO platform on a larger scale