Library to work load feedback files and perform a verification in R

# Download and install the FFV2 package
git clone https://gitlab.com/rfxf/ffv2.git
R CMD INSTALL FFV2

* Functions to work with feedback files in R
* Functions for individual verification tasks
* Script control via namelists to a large extent
* Verification can be run interactively or in batch mode
* Verification scripts are parallelized in many places
* Results in score files to be used in dedicated visualization apps
* Shell scripts in "/demo" folder are recommended to run verification jobs 
  on a regular basis

Domain average verification

Domain average verification summary

Domain average time-series verification

Domain average station based verification

Domain average TEMP verification

Feedback Files

* Filename should contain valid-date  (YYYYMMDDHH)
* Filename should indicate observation type (e.g. SYNOP, TEMP, etc.)
* Filename should indicate forecast type (ver, vep)
* All files of one experiment in one separate directory
* Different directories for different experiments

Other

* Additionally the verification can be refined using
* Station lists
* Spatial polygon lists
* Observation whitelists
* Datelists

Domain average verification

# Calling the R verification script with arguments 
# namelist and number of cores

Rscript YOURPATH/FFV2/demo/starter_scores_by_date.R \ 
        YOURPATH/FFV/namelists/namelist_verSYNOP.nl SYNOP DET 6

* Rscript to run a script in batch mode
* SYNOP or TEMP verification
* Ensemble (EPS) or deterministic (EPS) 
* Path and name of namelist
* Number of cores (default 1)
* Results in one intermediate score file per feedback file
* Existing score files will not be replaced

Domain average verification

# Calling the R aggregation script with arguments
# namelist and number of cores

Rscript YOURPATH/FFV2/demo/starter_aggregate.R \
        YOURPATH/FFV/namelists/namelist_verSYNOP.nl SYNOP DET 6

* Rscript to run a script in batch mode
* SYNOP or TEMP verification
* Ensemble (EPS) or deterministic (EPS) 
* Path and name of namelist
* Number of cores (default 1)
* Results in at least one final score file per verification job
* Final scores can be sent to visualization app
* Existing final score file will be replaced

Station based verification

# Calling the R station verification script with arguments 
# namelist and number of cores
Rscript YOURPATH/FFV2/demo/starter_scores_by_station.R \
        YOURPATH/FFV/namelists/namelist_verSYNOP.nl SYNOP DET CONT 6
Rscript YOURPATH/FFV2/demo/starter_scores_by_station.R \
        YOURPATH/FFV/namelists/namelist_verSYNOP.nl SYNOP DET CAT 6

* Rscript to run a script in batch mode
* Script for SYNOP or TEMP verification
* Ensemble (EPS) or deterministic (EPS)
* Continuous (CONT) or categorical (CAT) scores 
* Path and name of namelist
* Number of cores (default 1)
* Results directly in one final score file per verification job
* Final scores can be sent to visualization app
* Most memory demanding
* Existing score file will be replaced

SYNOP deterministic

varnoContinuous 'T2M,RH,RH2M,PS,TD2M,FF,DD,RAD_GL_1h,RAD_DF_1h,
                 TMIN_12h,TMAX_12h,JJ_12h,N,N_L,N_M,N_H'
pecthresholds   list('PS'=list('lower'=c(-500),'upper'=c(500)),
                     'FF'=list('lower'=c(-3),'upper'=c(3)))
catthresholds   list('GUST_1h'=c(5,10,15),'RR_1h'=c(0.1,1,2))
subdomains      'ALL,GER'
filePattern     'verSYNOP.'
sigTest         'T'
veriSetting     'COSMO-DE'
expType         'standard'
experiment      'T'
expIds          'Exp1,Exp2'
expDescr        'Exp1,Exp2'
fileDescr       'Exp_Exp1_vs_Exp2_username'
fdbkDirs        '/Directory/Of/Exp1/,/Directory/Of/Exp2/'
outDir          '/Output/Directory/'

(Minimum namelist required (pecthresholds & catthresholds can be omitted))

Namelist key expIds

* ff usually don't have meaningful names for the forecast model
* Experiment names have to be set in the namelist via the key `expIds`
* Value has to be a comma separated character string
* Identical names will cause scripts to crash
* Longish names will make plots look awkward
* Names will be used for the results and plots

Example for `expIds`  'CDE,ILAM'

Namelist key experiment

Verification distinguishes between experiment mode an standard mode
using the namelist key `experiment` with values 'TRUE' or 'FALSE'

Experiment mode

* Evaluates only ff available for all experiments
* Verification period is that of shortest experiment

Standard mode

* Evaluates all ff
* Intermediate scores even for dates with missing experiments
* Allows to aggregate scores using a subset of experiments

Namelist keys expDescr & fileDescr

expDescr

Used to name the intermediate score files. 
This name is than looked after in the aggregation script

fileDescr

Used to name the final score file.
Note: Too long names might cause system errors

Namelist keys outdir & shinyServer & shinyPath

Set `outdir` as the working directory. Namelists, scores and logfiles can be found there


Set `shinyPath` to copy (scp) the final score files in the visualization app directory

Set `shinySever` for the server name of the shiny apps


Both keys (shinyPath & shinyServer) have to be set, even if the app is 
running on the same server
Passwordless login to the server is required
If both keys are not set, the final files remain in the working dir. (`outdir`)

Example:

shinyServer 'fdbk@remora.hnms.gr'
shinyPath   '/home/fdbk/ShinyApps/'

Namelist key expType

Used to perform hindcast or standard verification
In hindcast mode (set `expType` to 'singlerun'):

* Lead-times > 23 hours will be converted to time of day (0-23)
* Usually very fast, as each ff contains one forecast only
* All other features (sign. test, conditions, etc.) available
* Time-series and station based verification are produced as well

Namelist key subdomains

Used to further devide the verification domain in sub-domains.

Valid values are


* ALL - Entire domain
* LAT - NH,SH,TR
* SYNOP - predefined set of subdomains for SYNOP verif.
* WMO  - predifined set of stations form TEMP veri.
* CDE,CEU - CSOMO-DE, COSMO-EU domain
* GER - Germany
* USER - see next slides

Namelist key domainTable I

Points to a text file that holds polygon set(s) of sub-domain(s).
Example syntax:

name    lon     lat
NORD    8       50.001
NORD    15      50.001
NORD    15      55
NORD    8       55
SUED    8       45
SUED    15      45
SUED    15      50
SUED    8       50

In addition the key `subdomains` needs to contain 'USER'
Polygons can be complex

Namelist key domainTable II

Points to a text file that holds station list(s).
Example syntax:

    name    id
    GER     Q887
    GER     10837
    GER     10184
    CH      06670
    CH      06612
    CH      06610

Key `subdomains` needs to contain 'USER'

Namelist key whitelists

If the key `whitelists` is set, only the stations for the variables in the list are used.
All other variables remain unaffected.
Example for allowing only 3 stations for N and N_H:

"statid","varno"
"07024","N"
"07412","N"
"07134","N"
"07791","N_H"
"07770","N_H"
"07666","N_H"

Namelist key condition

Arbitrary number of conditions can be set
This implementation allow for conditions on variables only!
The condition number will be added to the experiment names (expIDs)

condition1  "list(N='obs==0',N='abs(veri_data-obs)<1')"
condition2  "list(N='obs==0',N='abs(veri_data-obs)>=1')"
condition3  "list(N='obs%between%c(1,4)',N='abs(veri_data-obs)<1')"
condition4  "list(N='obs%between%c(1,4)',N='abs(veri_data-obs)>=1')"
condition5  "list(N='obs%between%c(5,7)',N='abs(veri_data-obs)<1')"
condition6  "list(N='obs%between%c(5,7)',N='abs(veri_data-obs)>=1')"
condition7  "list(N='obs==8',N='abs(veri_data-obs)<1')"
condition8  "list(N='obs==8',N='abs(veri_data-obs)>=1')"

Namelist key sigTest

Significance test is performed if `sigTest` is set 'TRUE'
Implemented for SYNOP & TEMP (det.) domain average scores
Not implemented for station based scores
Test on sign. in diff. between median of score distributions (t-test, 95% c.l.)
Only for individual runs to avoid correlated scores
Costs extra run-time and memory
Returns as character string, for each exp. combination, 
    if test was positive (2), negative (1) or not possible (0)
Result is encoded in visualization app and used for coloring

Namelist keys veri_run_type & veri_ens_member

Type of forecasts & analyses can be selected using the keys

`veri_run_tpye`
based of table 20, ff documentation, forecast and different types of analysis
can be selected (default '0,2' for ICON or '0,4' for COSMO)

`veri_ens_member`
based on table 22, ff documentation, this key can be used to select between the
determinisic run ('-1'), ensemble members (e.g. '1' - '20') or the ensemble mean ('-7') 

Namelist keys state & mimicVersus

The default data quality settings (next slide) can be overridden with

`state`
Defines the status of the observations that will be used for verification
according to table 3 of the ff documentation

`mimicVersus`
ignores most of the pre-defined quality settings and uses the VERSUS
settings:
    
|fcst-obs| >= 50 m/s for surface windspeed
|fcst-obs| >= 25 hPa for surface pressure
|fcst-obs| >= 30 K for 2m dew-point or 2m temperature

The verification uses all observations that pass the DA quality control
as they are, with a few, hard-coded exceptions:

* Precip. observations >100mm/h are not used
* Precip. forecasts are rounded to the first decimal place
* Double reports from some USA stations are removed
* Wind observations >100m station height are used, although rejected by DA (COSMO)
* Wind observations with more than 100m difference between model and station height are not used
* Wind observation are discarded is difference to forecast is >40m/s
* Wind (and humidity sensitive) observations are not used if sso_stdh>70 (and >200)
* Temperature, humidity and pressure observations with more than 150m difference between model and station height are not used
* 2m Dew-point observation are discarded is difference to forecast is >=30K
* Wind direction is discarded if obs or fcst <3m/s
* Global radiation obs. >2000W are not used
* RH observations above 300hPa are used only from Vaisala sondes
* Ceiling observation are discarded is difference to forecast is >100m
* Observed gusts >75m/s are not used
* Diffuse radiation is set to 0 in case global radiation is 0

Problems?

* All verification scripts print some output to 
  indicate errors or the current status.
* Namelist syntax OK?
* Passwordless login to shiny server working?
* Missing data in feedback files?
* Too little memory requested?

Recommendations

* When looking at scores, also look at size of data 'LEN'
* Look at time-series of scores to see if dates are missing
* Look at station based verification to see if stations are missing
* If entire dates are missing, was it because data alignment
  was not possile (see *veri* script output)?
* Read the documentation for more namelist options and
  examples, e.g. on ensemble verification

R data.table

* Package on CRAN https://cran.r-project.org/web/packages/data.table/index.html
* Good start on https://github.com/Rdatatable/data.table/wiki
* Faster version of R data.frame
* Very efficient calculations with tables
* Frequently used in ff verification scripts

# SYNTAX
DT[i,j,k]

DT - data table (columns of equal length, mixed variable type allowed)
i - filter on rows
j - function using columns
k - group by categories in columns 

Using R data.table in verification

require(data.table)
DT = data.table(varno=c("T2M","TD2M"),obs=rnorm(10000),veri_data=rnorm(10000),statid=1:33)
print(DT)

##         varno         obs   veri_data statid
##        <char>       <num>       <num>  <int>
##     1:    T2M -0.46479241 -1.67122751      1
##     2:   TD2M -1.44929950  2.13268556      2
##     3:    T2M -0.79650888 -0.51635272      3
##     4:   TD2M -0.77189308  0.51123171      4
##     5:    T2M  0.38769547 -0.53496822      5
##    ---                                      
##  9996:   TD2M  1.22625182 -0.01656439     30
##  9997:    T2M  0.55240318  0.20173571     31
##  9998:   TD2M -1.41817025 -0.46954168     32
##  9999:    T2M -0.38948671  0.31524668     33
## 10000:   TD2M  0.06973799  0.20941130      1

Using R data.table in verification

Calculate BIAS
if |obs-fcst|<10°C
group by variable

require(data.table)
DT = data.table(varno=c("T2M","TD2M"),obs=rnorm(10000),veri_data=rnorm(10000),statid=1:33)
SCORES = DT[abs(veri_data-obs)<10,list(BIAS=mean(veri_data-obs)),by=c("varno")]
print(SCORES)

##     varno       BIAS
##    <char>      <num>
## 1:    T2M 0.01513599
## 2:   TD2M 0.01398364

Actual verification problems are much larger!
Example COSMO:
27 forecast steps, 1 runs/day, 30 days, 12 variables, 500 stations, 2 experiments
Data table with ~ 20 million rows