Scores for FLUXNET reference data in NetCDF format
This function compares model output against
measurements from FLUXNET or from the Baseline Surface Radiation Networking in
NetCDF format. The performance of a model is expressed through scores that range
from zero to one, where increasing values imply better performance. These scores
are computed in five steps:
(i) computation of a statistical metric,
(ii) nondimensionalization,
(iii) conversion to unit interval,
(iv) spatial integration, and
(v) averaging scores computed from different statistical metrics.
The latter includes the bias, root-mean-square error, phase shift,
inter-annual variability, and spatial distribution. The corresponding equations
are documented in amber-package.
scores.fluxnet.nc(long.name, nc.mod, nc.ref, mod.id, ref.id, unit.conv.mod,
unit.conv.ref, variable.unit, score.weights, rotate.me = TRUE,
irregular = FALSE,
my.projection = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0",
shp.filename = system.file("extdata/ne_110m_land/ne_110m_land.shp",
package = "amber"), my.xlim = c(-180, 180), my.ylim = c(-60, 85),
plot.width = 8, plot.height = 3.8, numCores = 2,
outputDir = FALSE, phaseMinMax = "phaseMax")long.name |
A string that gives the full name of the variable, e.g. 'Gross primary productivity' |
nc.mod |
A string that gives the path and name of the netcdf file that contains the model output, e.g. '/home/model_gpp.nc' |
nc.ref |
A string that gives the path and name of the netcdf file that contains the reference data output, e.g. '/home/reference_gpp.nc' |
mod.id |
A string that identifies the source of the reference data set, e.g. 'CanESM2' |
ref.id |
A string that identifies the source of the reference data set, e.g. 'MODIS' |
unit.conv.mod |
A number that is used as a factor to convert the unit of the model data, e.g. 86400 |
unit.conv.ref |
A number that is used as a factor to convert the unit of the reference data, e.g. 86400 |
variable.unit |
A string that gives the final units using LaTeX notation, e.g. 'gC m$^-2$ day$^-1$' |
score.weights |
R object that gives the weights of each score (S_{bias}, S_{rmse}, S_{phase}, S_{iav}, S_{dist}) that are used for computing the overall score, e.g. c(1,2,1,1,1) |
rotate.me |
logical: TRUE if you want longitudes to range from -180 to 180 degrees and FALSE if you want longitudes to range from 0 to 360 degrees |
irregular |
logical: TRUE if data is on an irregular grid and FALSE if data is on a regular grid |
my.projection |
A string that gives the projection of the irregular grid |
shp.filename |
A string that gives the coastline shapefile |
my.xlim |
An R object that gives the longitude range that you wish to plot, e.g. c(-180, 180) |
my.ylim |
An R object that gives the longitude range that you wish to plot, e.g. c(-90, 90) |
plot.width |
Number that gives the plot width, e.g. 8 |
plot.height |
Number that gives the plot height, e.g. 4 |
numCores |
An integer that defines the number of cores, e.g. 2 |
outputDir |
A string that gives the output directory, e.g. '/home/project/study'. The output will only be written if the user specifies an output directory. |
phaseMinMax |
A string (either 'phaseMax' or 'phaseMin') that determines whether to assess the seasonal peak as a maximum or a minimum. The latter may be appropriate for variables that tend to be negative, such as net longwave radiation or net ecosystem exchange. |
(1) Figures in PDF format that show global maps of the model data at the location of FLUXNET sites (mean, mod.mean; interannual-variability, mod.iav; month of annual cycle maximum, mod.max.month), the reference data (mean, ref.mean; interannual-variability, ref.iav; month of annual cycle maximum, ref.max.month), statistical metrics (bias, bias; centralized root mean square error, crmse; time difference of the annual cycle maximum, phase), and scores (bias score, bias.score; root mean square error score, rmse.score; inter-annual variability score iav.score; annual cycle score (phase.score).
(2) Four text files: (i) score values and (ii) score inputs for each individual site, and (iii) score values and (iv) score inputs averaged across sites. when averaging over all station.
library(amber)
library(classInt)
library(doParallel)
library(foreach)
library(Hmisc)
library(latex2exp)
library(ncdf4)
library(parallel)
library(raster)
library(rgdal)
library(rgeos)
library(scico)
library(sp)
library(stats)
library(utils)
library(viridis)
library(xtable)
long.name <- 'Gross primary productivity'
nc.mod <- system.file('extdata/modelRegular', 'gpp_monthly.nc', package = 'amber')
nc.ref <- system.file('extdata/referenceRegular', 'gpp_FLUXNET.nc', package = 'amber')
mod.id <- 'CLASSIC' # define a model experiment ID
ref.id <- 'FLUXNET' # give reference dataset a name
unit.conv.mod <- 86400*1000 # optional unit conversion for model data
unit.conv.ref <- 86400*1000 # optional unit conversion for reference data
variable.unit <- 'gC m$^{-2}$ day$^{-1}$' # unit after conversion (LaTeX notation)
score.weights <- c(1,2,1,1,1) # score weights of S_bias, S_rmse, S_phase, S_iav, S_dist
# global plot
scores.fluxnet.nc(long.name, nc.mod, nc.ref, mod.id, ref.id, unit.conv.mod,
unit.conv.ref, variable.unit, score.weights)
# regional plot
scores.fluxnet.nc(long.name, nc.mod, nc.ref, mod.id, ref.id, unit.conv.mod,
unit.conv.ref, variable.unit, score.weights, my.xlim = c(-150, -60), my.ylim = c(20, 60),
plot.width = 6, plot.height = 3.8)
scores.fluxnet.nc(long.name, nc.mod, nc.ref, mod.id, ref.id, unit.conv.mod,
unit.conv.ref, variable.unit, score.weights)
# (2) Example for data on a rotated grid
nc.mod <- system.file('extdata/modelRotated', 'gpp_monthly.nc', package = 'amber')
nc.ref <- system.file('extdata/referenceRegular', 'gpp_FLUXNET.nc', package = 'amber')
mod.id <- 'CLASSIC' # define a model experiment ID
ref.id <- 'FLUXNET' # give reference dataset a name
unit.conv.mod <- 86400*1000 # optional unit conversion for model data
unit.conv.ref <- 86400*1000 # optional unit conversion for reference data
variable.unit <- 'gC m$^{-2}$ day$^{-1}$' # unit after conversion (LaTeX notation)
score.weights <- c(1,2,1,1,1) # score weights of S_bias, S_rmse, S_phase, S_iav, S_dist
rotate.me <- FALSE
irregular <- TRUE
my.projection <-'+proj=ob_tran +o_proj=longlat +o_lon_p=83. +o_lat_p=42.5 +lon_0=263.'
# shp.filename <- system.file('extdata/ne_50m_admin_0_countries/ne_50m_admin_0_countries.shp',
# package = 'amber')
shp.filename <- system.file("extdata/ne_110m_land/ne_110m_land.shp", package = "amber")
my.xlim <- c(-171, 0) # longitude range that you wish to plot
my.ylim <- c(32, 78) # latitude range that you wish to plot
plot.width <- 7 # plot width
plot.height <- 3.8 # plot height
numCores = 2
scores.fluxnet.nc(long.name, nc.mod, nc.ref, mod.id, ref.id, unit.conv.mod,
unit.conv.ref, variable.unit, score.weights, rotate.me, irregular,
my.projection,
shp.filename, my.xlim, my.ylim)
#donttestPlease choose more modern alternatives, such as Google Chrome or Mozilla Firefox.