---
title: "Introduction to eva3dm"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to eva3dm}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
library(eva3dm)
library(riem)
library(terra)
```

The recommended workflow to evaluate your model results will include 4 steps:


1. Pre-processing of observations

2. Pre-processing of model output

3. Model Evaluation

4. Visualization

The exact steps needed for evaluation will be different for different models, different observation data-sets, and different variables. However, in this vignette is presented an example of evaluation of temperature at surface level (T2) from WRF-Chem for a simulation at 3-km of resolution for January 2016 using data from METeorological Aerodrome Report (METAR).

1 - Pre-processing of observations:
We are using the r-package `riem` to download METAR data:

```{r metar}
start_date  <- "2016-01-01"
end_date    <- "2016-02-01"
sites       <- c("SBGR","SBKP","SBMT","SBSJ","SBSP","SBST","SBTA")

METAR  <- data.frame(date = seq.POSIXt(as.POSIXct(start_date), 
                                       as.POSIXct(end_date), 
                                       by = "hour"))
for(site in sites){
  cat('Trying to download METAR from:',site,'...\n')

  DATA <- tryCatch(riem::riem_measures(station    = site,
                                       date_start = start_date,
                                       date_end   = end_date,
                                       elev       = FALSE,
                                       latlon     = FALSE),
                   error = NULL)
  
  if(is.null(DATA)){
    cat('fail to download, loading some data ...\n')
    METAR <- readRDS(paste0(system.file("extdata",package="eva3dm"),
                            "/METAR_MASP_jan_2016.Rds"))
    break
  }
  
  DATA <- data.frame(date = DATA$valid,
                     T2   = DATA$tmpf)
  names(DATA) <- c('date', site)
  METAR       <- merge(x     = METAR, 
                       y     = DATA, 
                       by    = "date", 
                       all   = T, 
                       sort  = TRUE)
}
```
Lets check if the observations are ok

```{r check, fig.width = 7, fig.height = 4}
plot(METAR$date, METAR[,2], ty = 'l',xlab = '',ylab = 'T2', main = 'METAR OBS')
head(METAR)
```
In this case the temperature units needs to be converted from Fahrenheit to Celsius and the data needs to be grouped or averaged in time, in this case we are using the `hourly()` to perform a evaluation against hourly data.
```{r observation, fig.width = 7, fig.height = 4}

METAR[,-1] <- 5/9 * (METAR[,-1]-32)
METAR      <- hourly(METAR)

plot(METAR$date, METAR[,2], ty = 'l',xlab = '',ylab = 'T2', main = 'METAR processed OBS')
head(METAR)

```

For the next step a site-list with latitude and longitude of each surface station is needed, a list of all METAR stations is available in eva3dm package.
```{r site-list}
site_list <- readRDS(paste0(system.file("extdata",package="eva3dm"),"/sites_METAR.Rds"))
head(site_list)
```

2 - Pre-processing of model output:
The model output are large files and does not allow the direct comparison with observations. The package provides a function that extracts time-series from geographic locations from a site-list. The function `extract_serie()` extracts the time-series from a list of files and saves the output in R file (in this example we are reading the R file generated by `extract_serie()`). An finally, the temperature extracted from the model is converted from Kelvin to Celsius.
```{r model}
## to extract time-series from WRF-Chem model
## wrf_files <- dir(pattern = "wrfout_d03")
## extract_serie(filelist = wrf_files, point = site_list, variable="T2", prefix="model.d03", field="3d")
model_d03 <- readRDS(paste0(system.file("extdata",package="eva3dm"),"/model.d03.T2.Rds"))
model_d03[-1] <- model_d03[-1] - 273.15
```
NOTE: to extract time-series from other models (not WRF or WRF-Chem), the arguments `latitude` and `longitude` from the function `extract_serie()` need to match the name of these variables for latitude and longitude in the model output. If the model do not include the Times variable, `use_TFLAG` must be set to `TRUE` to extract and convert the time from TFLAG variable (for example for CMAQ / CAMx models) or `use_datesec` must be set to `TRUE` to extract and convert time from the datasec variable (for example for WACCM / CAM-Chem models).

3 - Model Evaluation: 
The functions `eva()` can be used to calculate the statistical indexes from time-series for individual stations or all stations combined.
```{r evaluation}
table <- data.frame()
for(site in sites){
  table <- eva(mo = model_d03, ob = METAR, site = site, table = table)
}
table <- eva(mo = model_d03, ob = METAR, site = 'ALL', table = table)
print(table)
```

4 - Visualization:
The function `overlay()` can be used to plot the results from `eva()` on a map. However, first the data need to be georeferenced using `%at%`.
```{r visualize, fig.width = 7, fig.height = 5.5}
spatial_table <- table %at% site_list
overlay(spatial_table, z = 'MB', main = 'T2 main bias (MB)',expand = 1.6,lim = 0.1)
masp <- terra::vect(paste0(system.file("extdata",package="eva3dm"),"/masp.shp"))
BR   <- terra::vect(paste0(system.file("extdata",package="eva3dm"),"/BR.shp"))
terra::lines(BR)
terra::lines(masp, col = 'gray')
legend_range(spatial_table$MB,y = table["ALL","MB"])
```
The WRF-Chem temperature in this example showed a very low MB, except for the SBSJ station (MB=-0.7), which is in the criteria for acceptable performance (|MB| < 0.5) and the criteria for evaluation for complex terrain for the SBSJ case (|MB| < 1.0). Another criteria for temperature are IAO > 0.8 and ME < 2.0 (ME < 3, for complex terrain), in general the model presented a reasonable performance for temperature. More information and references can be found on the documentation of `stat()` function.