---
title: "Import hydrographic data into MFDB"
output:
  html_document:
    toc: true
vignette: >
  %\VignetteIndexEntry{Hydrographic data}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
```{r, message=FALSE, echo=FALSE}
library(unittest)
# Redirect ok() output to stderr
options(unittest.output = stderr())
```
```{r, message=FALSE}
library(mfdb)
```

```{r, message=FALSE, echo=FALSE}
# Don't show chatty MFDB messages
logging::setLevel('WARN')
```

Hydrographic products are available from numerous sources. Here, we provide a working example where outputs from a hydrographic model are:

* downloaded from Copernicus Marine Service https://marine.copernicus.eu/
* aggregated over ICES rectangles
* imported into MFDB


# Download monthly hydrographic data
CMEMS Baltic Sea Physical Reanalysis product provides a physical reanalysis for the whole Baltic Sea area from January 1993 and up to minus 1-1.5 year compared to real time. The product is produced by using the ice-ocean model NEMO-Nordic (based on NEMO-3.6, Nucleus for European Modelling of the Ocean) together with a LSEIK data assimilation scheme,
https://data.marine.copernicus.eu/?option=com_csw&view=details&product_id=BALTICSEA_REANALYSIS_PHY_003_011
As a demonstration we retrieve one year of already processed monthly means of bottom temperature (bottomT), but hydrographic products are also available for separate depth layers and at higher spatiotemporal resolution, and can be aggregated before import according to the needs. Note that an account has to be created to download data from Copernicus Marine Service using the script below.


```{r eval=FALSE}
library(RCMEMS)

# point to the data product
cfg <- CMEMS.config(motu="http://my.cmems-du.eu/motu-web/Motu",
                         service.id = "BALTICSEA_REANALYSIS_PHY_003_011-TDS",
                         product.id = "dataset-reanalysis-nemo-monthlymeans",
                         variable = c("bottomT"))

# add user and psw
CMEMS.config.usr <- function(x){
    print("username")
    scan("", what="character", nmax=1, quiet=T)
    }
CMEMS.config.pwd <- function(x){
    print("password")
    scan("", what="character", nmax=1, quiet=T)
    }
cfg@user <- CMEMS.config.usr()
cfg@pwd  <- CMEMS.config.pwd()

# select year and download
y <- 2018 # year of interest
CMEMS.download(cfg,
                   ROI = c(17,20,56,58.5),
                   date.range = c(ISOdate(y,01,01), ISOdate(y,12,31)),
                   depth.range= c(10,500), # max depth Baltic 459 m
                   out.path=paste("data_provided/CMEMS_BAL_PHY_reanalysis_monthlymeans_",y,".nc",sep=""),
               debug=FALSE)
```

# Download ICES rectangles

```{r eval=FALSE}
temp <- tempfile()
download.file(url="https://gis.ices.dk/shapefiles/ICES_StatRec_mapto_ICES_Areas.zip", temp)
unzip(temp, exdir="data_provided/ices_rect")
unlink(temp)
```


# Aggregate hydrographic data by ICES rectangles

```{r eval=FALSE}
library(sf)
library(raster)
library(tidyverse)
library(rnaturalearth)
library(lwgeom)
library(ggplot2)

y <- 2018 # year of interest

# load ICES rect
ices.rect <- st_read("data_provided/ices_rect/StatRec_map_Areas_Full_20170124.dbf", quiet = T)

# load bottom temperature (bottomT)
rst <- brick(paste("data_provided/CMEMS_BAL_PHY_reanalysis_monthlymeans_",y,".nc",sep=""), varname="bottomT", lvar=4)

# calculate mean sob (raster) by ices rect (polygons)
ov <- raster::extract(rst,ices.rect, fun=mean, na.rm=T, df=T)
hydroVar <- data.frame(ID=1:length(ices.rect$ICESNAME),
                       rect=ices.rect$ICESNAME,
					   areaFull_km2=ices.rect$AREA_KM2,
					   year=y) %>%
    right_join(ov) %>%
    select(-ID) %>%
    gather("month","sob",4:15) %>%
    mutate(month = as.numeric(substring(month,7,8))) %>%
    subset(!is.na(sob)) # NA are on land or outside the raster area

land <- rnaturalearth::ne_countries(returnclass = "sf", continent="Europe", scale="large") %>%
  st_union()
ices.rect.sea <- ices.rect %>% filter(ICESNAME %in% hydroVar$rect) %>%
         st_difference(land)
ices.rect.sea$areaSea_km2 <- st_area(ices.rect.sea) %>% units::set_units(km^2)

hydroVar <- ices.rect.sea %>%
    rename(rect=ICESNAME) %>%
    select(rect, areaSea_km2) %>%
    right_join(hydroVar)

ggplot(hydroVar) +
    geom_sf() +
    geom_sf(data=st_geometry(ices.rect.sea), fill="lightblue") +
    xlim(xmin(extent(hydroVar)), xmax(extent(hydroVar))) +
    ylim(ymin(extent(hydroVar)), ymax(extent(hydroVar))) +
    theme_bw()	

# Write the results to temporary files so we can read them back in the next step

write.table(sf::st_drop_geometry(ices.rect.sea), file = 'ices.rect.sea.txt')
write.table(sf::st_drop_geometry(hydroVar), file = 'hydroVar.txt')
```

# Import bottom temperature into mfdb

Import ICES rectangles and import the hydrographic data as a survey index. Then we can retrieve the hydrographic data as a weighted mean using ICES rectangle as a weighting variable:

```{r}
mdb <- mfdb(tempfile(fileext = '.duckdb'))

ices.rect.sea <- read.table('ices.rect.sea.txt')
hydroVar <- read.table('hydroVar.txt')

# import ICES rectangles
mfdb_import_area(mdb, data.frame(
    name = as.character(ices.rect.sea$ICESNAME),
    size = ices.rect.sea$areaSea_km2))

# import area size as index for later use
mfdb_import_cs_taxonomy(mdb, "index_type", data.frame(name=c("ices_rect", "bottom_temp")))
mfdb_import_survey_index(mdb,
                         data_source="area_ices_rect",
                         data.frame(index_type = "ices_rect",
                                        year       = hydroVar$year,
                                        month      = hydroVar$month,
                                        areacell   = as.character(hydroVar$rect),
                                        value      = as.numeric(hydroVar$areaSea_km2)))

# import avg bottom temperature by month and ICES rect (under the 'length' field)
mfdb_import_survey_index(mdb,
                   data_source="bottom_temp_CMEMS",
                   data.frame(year       = hydroVar$year,
                              month      = hydroVar$month,
                              areacell   = as.character(hydroVar$rect),
                              value     = hydroVar$sob,
                              index_type = "bottom_temp",
                              stringsAsFactors = TRUE))

# extract mean bottom temperature by quarter and area (weighted mean using rectangle area)
dat <- mfdb_survey_index_mean(mdb, c(), list(
                                       ## area = mfdb_unaggregated(),
                                       area = mfdb_group('a1'=c('42G7','42G8','42G9','43G7','43G8','43G9'),
                                                         'a2'=c('44G7','44G8','44G9','45G7','45G8','45G9')),
                                       timestep = mfdb_timestep_quarterly,
                                       year = 2018,
                                       index_type = "bottom_temp",
                                       data_source = 'bottom_temp_CMEMS'),
                              scale_index = 'area_size')[[1]]
dat <- dat[,c('year', 'step', 'area', 'mean')]
dat
```

```{r, echo=FALSE}
# NB: Calculated with:
# raw <- hydroVar[hydroVar$year == 2018 & hydroVar$rect %in% c('44G7','44G8','44G9','45G7','45G8','45G9') & hydroVar$month %in% c(1,2,3),]
# sum(raw$sob * raw$areaSea_km2) / sum(raw$areaSea_km2)
ok(ut_cmp_equal(dat, read.table(blank.lines.skip = TRUE, header = TRUE, stringsAsFactors = FALSE, text = '
year step area     mean
2018    1   a1 5.131353
2018    1   a2 5.286023
2018    2   a1 5.772875
2018    2   a2 5.787779
2018    3   a1 7.259383
2018    3   a2 6.609905
2018    4   a1 6.144196
2018    4   a2 6.020029
', colClasses = c(NA, 'character', NA, NA)), tolerance = 1e-5), "extract mean bottom temperature by quarter and area")
```

```{r}
mfdb_disconnect(mdb)
options(unittest.output = NULL)
```
