Type:	Package
Title:	Evaluation of 3D Meteorological and Air Quality Models
Version:	1.11
Date:	2025-06-02
Description:	Provides tools for post-process, evaluate and visualize results from 3d Meteorological and Air Quality models against point observations (i.e. surface stations) and grid (i.e. satellite) observations.
Maintainer:	Daniel Schuch <underschuch@gmail.com>
License:	MIT + file LICENSE
Imports:	terra, ncdf4, utils
Suggests:	knitr, riem, rmarkdown, testthat (≥ 3.0.0)
Encoding:	UTF-8
BugReports:	https://github.com/Schuch666/eva3dm/issues/
RoxygenNote:	7.3.2
Depends:	R (≥ 3.5.0)
URL:	https://schuch666.github.io/eva3dm/
Config/testthat/edition:	3
VignetteBuilder:	knitr
NeedsCompilation:	no
Packaged:	2025-06-02 22:34:27 UTC; schuch
Author:	Daniel Schuch [aut, cre]
Repository:	CRAN
Date/Publication:	2025-06-02 23:00:01 UTC

eva3dm: Evaluation of 3D Meteorological and Air Quality Models

Description

Provides tools for post-process, evaluate and visualize results from 3d Meteorological and Air Quality models against point observations (i.e. surface stations) and grid (i.e. satellite) observations.

Author(s)

Maintainer: Daniel Schuch underschuch@gmail.com (ORCID)

See Also

Useful links:

• https://schuch666.github.io/eva3dm/

• Report bugs at https://github.com/Schuch666/eva3dm/issues/

Returns the common columns

Description

results of 'd01 in d02' style syntax

Usage

x %IN% y

Arguments

Value

data.frame with common columns or a cropped SpatRaster

Note

A message is always displayed to keep easy to track and debug issues (with the results and the evaluation process).

Can be used to crop rast objects, such as arguments of sat() function

See Also

See select for selection based on time.

Examples

times <- seq(as.POSIXct('2024-01-01',tz = 'UTC'),
             as.POSIXct('2024-01-02',tz = 'UTC'),
             by = 'hour')
randon_stuff <- rnorm(25,10)

observation <- data.frame(date   = times,
                          site_1 = randon_stuff,
                          site_3 = randon_stuff,
                          site_4 = randon_stuff,
                          site_5 = randon_stuff,
                          site_6 = randon_stuff,
                          site_7 = randon_stuff)

model_d01 <- data.frame(date  = times,
                        site_1=randon_stuff+1,
                        site_2=randon_stuff+2,
                        site_3=randon_stuff+3,
                        site_4=randon_stuff+4)

model_d02 <- data.frame(date  = times,
                        site_1=randon_stuff-1,
                        site_3=randon_stuff-3)

# multiline
model_d01_in_d02 <- model_d01 %IN% model_d02
eva(mo = model_d01_in_d02, ob = observation, rname = 'd01 in d02')

# or single line
eva(mo = model_d01 %IN% model_d02, ob = observation, rname = 'd01 in d02')
# or
eva(mo = model_d01, ob = observation %IN% model_d02, rname = 'd01 in d02')

Combine stats and site list to overlay plot

Description

combines the stats (from individual station evaluation) and site list in a SpatVector using row.names

Usage

stat %at% site

Arguments

Value

SpatVector (terra package)

Examples

sites <- data.frame(lat = c(-22.72500,-23.64300,-20.34350),
                    lon = c(-47.34800,-46.49200,-40.31800),
                    row.names = c('Americana','SAndre','VVIbes'),
                    stringsAsFactors = F)
model<- readRDS(paste0(system.file("extdata",package="eva3dm"),"/model.Rds"))
obs  <- readRDS(paste0(system.file("extdata",package="eva3dm"),"/obs.Rds"))

# evaluation by station
stats <- eva(mo = model, ob = obs, site = "Americana")
stats <- eva(mo = model, ob = obs, site = "SAndre",table = stats)
stats <- eva(mo = model, ob = obs, site = "VVIbes",table = stats)
# evaluation using all stations
stats <- eva(mo = model, ob = obs, site = "ALL", table = stats)

print(stats)

geo_stats <- stats %at% sites

print(geo_stats)

Read and write attributes on a NetCDF file

Description

Read and write metadata information of a NetCDF files

Usage

atr(file = NA, var = "?", att = NA, action = "get", value = NA, verbose = TRUE)

Arguments

Value

string with the NetCDF attribute value

Examples

nc <- paste0(system.file("extdata",package="eva3dm"),'/wrfinput_d01')
atr(nc,0)
atr(nc,'Times')
atr(nc,'XLAT')
atr(nc,'XLONG')

atr(nc,'XLONG','MemoryOrder')
atr(nc,'XLONG','description')
atr(nc,'XLONG','units')
atr(nc,'XLONG','stagger')
atr(nc,'XLONG','FieldType')

Calculate categorical statistics in related to a threshold

Description

Calculate traditional statistics related to a threshold

Usage

cate(
  model,
  observation,
  threshold,
  cutoff = NA,
  nobs = 8,
  rname,
  to.plot = FALSE,
  col = "#4444bb",
  pch = 19,
  lty = 3,
  lcol = "#333333",
  lim,
  verbose = TRUE,
  ...
)

Arguments

Value

a data.frame including: Accuracy (A); Critical Success Index (CSI); Probability of Detection (POD); Bias(B); False Alarm Ratio (FAR); Heidke Skill Score (HSS); Pearce skill Score (PSS) in

References

Yu, S., Mathur, R., Schere, K., Kang, D., Pleim, J., Young, J., ... & Rao, S. T. (2008). Evaluation of real‐time PM2. 5 forecasts and process analysis for PM2. 5 formation over the eastern United States using the Eta‐CMAQ forecast model during the 2004 ICARTT study. Journal of Geophysical Research: Atmospheres, 113(D6).

Examples

data <- 0.02 * 1:100
set.seed(666)
model  <- abs(rnorm(100,0.01))

oldpar <- par(pty="s")
cate(model = model, observation = data, threshold = 1,
     to.plot = TRUE, rname = 'example')
par(oldpar)

Calculate daily mean, min or max

Description

function to calculate daily mean, min or max of a data.frame

Usage

daily(
  data,
  time = "date",
  var,
  stat = mean,
  min_offset = 0,
  hour_offset = 0,
  numerical = TRUE,
  verbose = TRUE
)

Arguments

Value

data.frame with time and the daily mean, min or max

Examples

# in case there is connection issue
load_data <- function(cond) {
  message(paste("conection issue, loading pre-downloaded data"))
  DATA <- readRDS(paste0(system.file("extdata",package="eva3dm"),
                         "/riem_OAKB_jan_2012.Rds"))
  return(DATA)
}

sites <- c("OAKB")
for(site in sites){
  cat('Trying to download METAR from:',site,'...\n')
  DATA <- tryCatch(riem::riem_measures(station    = sites,
                                       date_start = "2012-01-01",
                                       date_end   = "2012-02-01"),
                   error = load_data)
}
data_daily_mean <- daily(DATA,time = 'valid')
data_daily_min  <- daily(DATA[1:7],time = 'valid',stat = min)
data_daily_max  <- daily(DATA[1:7],time = 'valid',stat = max)

Model statistical evaluation

Description

Statistical (or categorical) evaluation from 2 data.frames. The input data.frames (model and observation) must contain a "date" column (containing POSIXlt). The function perform some simple case tests and perform the time pairing of observations and model data and can calculate the statistical evaluation or categorical evaluation.

Usage

eva(
  mo,
  ob,
  rname = site,
  table = NULL,
  site = "ALL",
  wd = FALSE,
  fair = NULL,
  cutoff = NA,
  cutoff_NME = NA,
  no_tz = FALSE,
  nobs = 8,
  eval_function = stat,
  select_time,
  time = "date",
  remove_ch = FALSE,
  verbose = TRUE,
  ...
)

Arguments

Value

data.frame with statistical values from stat or cate functions.

Note

fair can be a data.frame or a character string to be used for the analysis, alternatively the function

for wind direction a rotation of 360 (or -360) is applied to minimize the wind direction difference.

If site == 'ALL' (default) all the columns from observations are combined in one column (same for observation) and all the columns are evaluated together.

Special thanks to Kiarash and Libo to help to test the wind direction option.

See Also

stat for additional information about the statistical evaluation and cate for categorical evaluation.

Examples

model <- readRDS(paste0(system.file("extdata",package="eva3dm"),
                        "/model.Rds"))
obs   <- readRDS(paste0(system.file("extdata",package="eva3dm"),
                        "/obs.Rds"))

# if there is no observed data
# the function return an empty row
table <- eva(mo = model, ob = obs, site = "VVIbes")
print(table)

# if the site are not in the input data frame a message is displayed
# and the function return an empty row
table <- eva(mo = model, ob = obs, site = "Ibirapuera")
print(table)

# calculating statistical with a few observed values
table <- eva(mo = model, ob = obs, site = "Americana")
print(table)

# calculating categorical (using 2 for threshold) with a few observed values
table <- eva(mo = model, ob = obs, site = "Americana",
             eval_function = cate, threshold = 2)
print(table)

# calculating categorical (using 2 for threshold) with a few observed values
table <- eva(mo = model, ob = obs, site = "Americana",
             eval_function = cate, threshold = 10)
print(table)

# customizing the evaluation function: inclusion of p.value from stats::cor.test()
stat_p <- function(x, y, ...){
  table         <- eva3dm::stat(x, y, ...)
  cor.result    <- stats::cor.test(x, y, ... )
  table$p.value <- cor.result$p.value
  table         <- table[,c(1:4,12,5:11)]
  return(table)
}

table <- eva(mo = model, ob = obs, site = "Americana",eval_function = stat_p)
print(table)

Create a NetCDF file with the surface maximum of O3

Description

Read the values from o3 and T2, convert o3 to ug m-3 and calculate the maximum of 8-hour moving avarage from a list of files.

Usage

extract_max_8h(
  filelist,
  variable = "o3",
  field = "4d",
  prefix = "max_8h",
  units = "ug m-3",
  meta = TRUE,
  filename,
  verbose = TRUE
)

Arguments

Value

No return value

Note

The field argument '4d' / '2dz' is used to read a 4d/3d variable droping the 3rd dimention (z).

Examples

dir.create(file.path(tempdir(), "MDA8"))
folder <- system.file("extdata",package="eva3dm")
wrf_file <- paste0(folder,"/test_small_o3.nc")
extract_max_8h(filelist = wrf_file,
               prefix = paste0(file.path(tempdir(),"MDA8"),'/mean'),
               field = '3d')

Create a NetCDF file with the surface mean

Description

Read and calculate the mean value of a variable from a list of wrf output files.

Usage

extract_mean(
  filelist,
  variable = "o3",
  field = "4d",
  prefix = "mean",
  units = "ppmv",
  meta = TRUE,
  filename,
  verbose = TRUE
)

Arguments

Value

No return value

Note

The field argument '4d' / '2dz' is used to read a 4d/3d variable droping the 3rd dimention (z).

Examples

dir.create(file.path(tempdir(), "MEAN"))
folder <- system.file("extdata",package="eva3dm")
wrf_file <- paste0(folder,"/wrf.day1.o3.nc")
extract_mean(filelist = wrf_file,prefix = paste0(file.path(tempdir(),"MEAN"),'/mean'))

Extract time series of wrf file list of lat/lon

Description

Read and extract data from a list of wrf output files and a table of lat/lon points based on the distance of the points and the center of model grid points, points outside the domain (and points on domain boundary) are not extracted.

Usage

extract_serie(
  filelist,
  point,
  variable = "o3",
  field = "4d",
  level = 1,
  prefix = "serie",
  new = "check",
  return.nearest = FALSE,
  fast = FALSE,
  use_ij = FALSE,
  latitude = "XLAT",
  longitude = "XLONG",
  use_TFLAG = FALSE,
  use_datesec = FALSE,
  id = "id",
  remove_ch = FALSE,
  verbose = TRUE
)

Arguments

Value

No return value

Note

The field argument '4d' or '2dz' is used to read a 4d/3d variable droping the 3rd dimention (z).

new = TRUE create a new file, new = FALSE append the data in a old file, and new = 'check' check if the file exist and append if the file exist and create if the file doesnt exist

FOR CAMx time-series, use the options: use_TFLAG=T, latitude='latitude', longitude='longitude', new=T

FOR WACCM time-series, use the options: use_datesec=T, latitude='lat', longitude='lon', new=T

The site-list of two global data-sets (METAR and AERONET) are provided on examples and site-list for stations on Brazil (INMET and Air Quality stations).

Examples

cat('Example 1: METAR site list\n')
sites <- readRDS(paste0(system.file("extdata",package="eva3dm"),"/sites_METAR.Rds"))

cat('Example 2: Integrated Surface Dataset (ISD) site list\n')
sites <- readRDS(paste0(system.file("extdata",package="eva3dm"),"/sites_ISD.Rds"))

cat('Example 4: AERONET site list\n')
sites <- readRDS(paste0(system.file("extdata",package="eva3dm"),"/sites_AERONET.Rds"))

cat('Example 5: list of INMET stations on Brazil\n')
sites <- readRDS(paste0(system.file("extdata",package="eva3dm"),"/sites_INMET.Rds"))

cat('Example 6: list of Air Quality stations on Brazil\n')
sites <- readRDS(paste0(system.file("extdata",package="eva3dm"),"/sites_AQ_BR.Rds"))

files    <- dir(path = system.file("extdata",package="eva3dm"),
                pattern = 'wrf.day',
                full.names = TRUE)
dir.create(file.path(tempdir(),"SERIE"))
folder <- file.path(tempdir(),"SERIE")

# extract data for 3 locations
extract_serie(filelist = files, point = sites[1:3,],prefix = paste0(folder,'/serie'))

Get the distance in kilometers between two points

Description

Get the distance in kilometers between two points

Usage

get_distances(lat1, long1, lat2, long2, R = 6371)

Arguments

Value

A numeric vector with the distance in kilometers.

#' source: https://github.com/gustavobio/brclimate/blob/master/R/get_distances.R

Calculate hourly mean, min or max

Description

function to calculate Ozone Maximum Daily 8-hr Average or 8-hr moving Average for a data.frame

Usage

hourly(
  data,
  time = "date",
  var,
  stat = mean,
  min_offset = 30,
  numerical = TRUE,
  verbose = TRUE
)

Arguments

Value

data.frame including only numerical columns

data.frame with time and the hourly mean, min or max

Examples

# in case there is connection issue
load_data <- function(cond) {
  message(paste("conection issue, loading pre-downloaded data"))
  DATA <- readRDS(paste0(system.file("extdata",package="eva3dm"),
                         "/riem_OAHR_jan_2012.Rds"))
  return(DATA)
}

sites <- c("OAHR")
for(site in sites){
  cat('Trying to download METAR from:',site,'...\n')
  DATA <- tryCatch(riem::riem_measures(station    = sites,
                                       date_start = "2012-01-01",
                                       date_end   = "2012-02-01"),
                   error = load_data)
}
data_hourly_mean <- hourly(DATA,time = 'valid')
data_hourly_min  <- hourly(DATA[1:7],time = 'valid',stat = min)
data_hourly_max  <- hourly(DATA[1:7],time = 'valid',stat = max)

Interpolation (project and resample)

Description

function to project and interpolate rast

Usage

interp(x, y, method = "bilinear", mask, verbose = FALSE)

Arguments

Value

SpatRaster (terra package)

Examples

model_o3 <- terra::rast(paste0(system.file("extdata",package="eva3dm"),
                              "/camx_no2.Rds"))
omi_o3   <- terra::rast(paste0(system.file("extdata",package="eva3dm"),
                              "/omi_no2.Rds"))

# interpolate omi O3 column to model grid
omi_o3c_interp_model <- interp(omi_o3,model_o3)

# interpolate model o3 column to omi grid
model_o3c_interp_omi <- interp(omi_o3,model_o3)

Plot a legend with the range of values

Description

Plot a legend with the range of values

Usage

legend_range(
  x,
  y,
  text.width = NULL,
  dig = c(2, 2, 2),
  xjust = 0.005,
  yjust = 0.95,
  horiz = TRUE,
  y.intersp = 0.5,
  x.intersp = 0.5,
  show.mean = TRUE,
  unit = "",
  label_mean = "ALL:",
  ...
)

Arguments

Value

No return value

Note

for use with rast use before any change of projection

text.width can vary depending on map dimensions

Examples

x <- 1:10 + rnorm(10,sd = .4)
plot(x,ty='l')
legend_range(x)

Calculate 8-hour moving average

Description

function to calculate Ozone 8-hour moving average for a data.frame

Usage

ma8h(data, time = "date", var, verbose = TRUE, ...)

Arguments

Value

data.frame with time and the 8-hour moving average

See Also

mda8 for Maximum Daily 8-hour moving average

Examples

model_file <- paste(system.file("extdata", package = "eva3dm"),
                                "/model_o3_ugm3_36km.Rds", sep="")
model      <- readRDS(model_file)
model_8h   <- ma8h(model)
plot(model$date,model$Campinas, pch = 19,
     main = expression(O[3]~~'['*mu*g*m^-3*']'))
points(model_8h$date,model_8h$Campinas, col = 'blue', pch = 19)
legend('topleft',bty = 'n',
       pch = 19,
       legend = c('hourly','8h-mov average'),
       col = c('black','blue'))

Maximum Daily 8-hr Average

Description

function to calculate Ozone Maximum Daily 8-hr Average or 8-hr moving Average for a data.frame

Usage

mda8(data, time = "date", var, verbose = TRUE)

Arguments

Value

data.frame with time and the maximum daily 8-hr average

See Also

ma8h for 8-hour Moving Average

Examples

model_file <- paste(system.file("extdata", package = "eva3dm"),
                    "/model_o3_ugm3_36km.Rds", sep="")
model      <- readRDS(model_file)
model_mda8 <- mda8(model)
model_8h   <- ma8h(model)
plot(model$date,model$Campinas, pch = 19,
     main = expression(O[3]~~'['*mu*g*m^-3*']'))
points(model_8h$date,model_8h$Campinas, col = 'blue', pch = 19)
points(model_mda8$date + 17*60*60,model_mda8$Campinas,
       col = 'red', pch = 4, cex = 2)
legend('topleft',bty = 'n',
       pch = c(19,19,4),
       legend = c('hourly','8h-mov average','MD8A'),
       col = c('black','blue','red'))

Print a 'ncdump -h' command

Description

Read a NetCDF and print the medatada

Usage

ncdump(file = file.choose())

Arguments

Value

No return value, only display information

Examples

ncdump(file = paste0(system.file("extdata",package="eva3dm"),
                     '/wrfinput_d01'))

Plot or add points using a color scale

Description

Custon plot for SpatRaster (terra R-package) object based on terra package

Usage

overlay(
  p,
  z,
  col,
  col2,
  lim = range(z, na.rm = TRUE),
  symmetry = TRUE,
  pch = 19,
  pch2 = 1,
  cex = 1,
  cex2 = 1.2 * cex,
  outside = TRUE,
  add = FALSE,
  plg = list(tic = "none", shrink = 1),
  pax = list(),
  expand = 1.15,
  ...
)

Arguments

Value

No return value

Examples

sp<- terra::vect(paste0(system.file("extdata",package="eva3dm"),"/masp.shp"))
BR<- terra::vect(paste0(system.file("extdata",package="eva3dm"),"/BR.shp"))

p    <- readRDS(paste0(system.file("extdata",package="eva3dm"),"/sites_AQ_BR.Rds"))
p$id      <- row.names(p)
point     <- terra::vect(p)
point$NMB <- 1:45 - 20 # some values to plot

terra::plot(BR, main = 'add points',xlim = c(-52,-37),ylim = c(-25,-18))
terra::lines(BR)
terra::lines(sp, col = 'gray')
overlay(point,point$NMB,cex = 1.4, add = TRUE)

overlay(point,point$NMB,cex = 1.4, add = FALSE, main = 'new plot')
terra::lines(BR)
terra::lines(sp, col = 'gray')

Plot the difference from two SpatRaster objects

Description

Custom difference (x - y) plots for SpatRaster object (based on terra package)

Usage

plot_diff(
  x,
  y,
  col,
  absolute = TRUE,
  relative = TRUE,
  lim_a = NA,
  lim_r = NA,
  scale,
  unit = c(units(x), expression("%")),
  ...
)

Arguments

Value

No return value

Examples

folder <- system.file("extdata",package="eva3dm")
wrf    <- paste0(folder,"/wrfinput_d01")
A      <- wrf_rast(wrf,'XLAT')
terra::units(A) <- 'degrees'
B      <- wrf_rast(wrf,'XLONG')
plot_diff(A,B,int = 2)

Plot rast (SpatRaster) object

Description

Custon plot for SpatRaster (terra R-package) object based on terra package

Usage

plot_rast(
  r,
  color,
  ncolor = 21,
  proj = FALSE,
  plg = list(tic = "none", shrink = 1),
  pax = list(),
  latitude = TRUE,
  longitude = TRUE,
  int = 10,
  grid = FALSE,
  grid_int = int,
  grid_col = "#666666",
  add_range = FALSE,
  ndig = 2,
  log = FALSE,
  range,
  scale,
  min = -3,
  max,
  unit,
  ...
)

Arguments

Value

No return value

Note

color scales including: 'eva3', 'eva4', 'blues', 'diff', and 'rain'. Also reverse version with addition of a r ('eva3r' is the default).

Examples

wrf <- paste(system.file("extdata", package = "eva3dm"),
                         "/wrfinput_d01", sep="")

r <- wrf_rast(file=wrf, name='XLAT')

plot_rast(r)

Convert absolute humidity to relative humidity

Description

function to convert absolute humidity to relative humidity.

Usage

q2rh(q, t = 15, p = 101325)

Arguments

Value

vector or data.frame with time and the relative humidity, units are

Note

default values are from standard atmosphere (288.15 K (15C) / 101325 Pa)

if rh and temp arguments are data.frame, both need to have the same number of lines and columns, first column (time column) will be ignored.

Examples

# for a single value (or same length vectors)
q2rh(q = 0.0002038, t = 29.3, p = 100800)

# using all data.frames
times <- seq(as.POSIXct('2024-01-01',tz = 'UTC'),
             as.POSIXct('2024-01-02',tz = 'UTC'),
             by = 'hour')[1:5]
q2   <- data.frame(time = times, a = rep(0.0002038,5))
temp <- data.frame(time = times, a = rep(     29.3,5))
pres <- data.frame(time = times, a = rep(   100800,5))
q2rh(q = q2, t = temp, p = pres)

# using data.frame for q and t (p is cte.)
q2rh(q = q2, t = temp, p = 100000)

# using data.frame for q and p (t is cte.)
q2rh(q = q2, t = 26, p = pres)

# using data.frame only for q (p and t are cte.)
q2rh(q = q2, t = 26, p = 100000)

conversion of model precipitation to hourly precipitation

Description

function that converts model accumulated precipitation to hourly precipitation.

Usage

rain(rainc, rainnc, verbose = TRUE)

Arguments

Value

data.frame time and the hourly precipitation or SpatRaster hourly precipitation

Examples

times <- seq(as.POSIXct('2024-01-01',tz = 'UTC'),
             as.POSIXct('2024-01-01 04:00:00',tz = 'UTC'),
             by = 'hour')
RNC   <- data.frame(date = times, aa = c(0.149,0.149,0.149,0.149,0.149))
RNNC  <- data.frame(date = times, aa = c(0.919,1.0,1.1,1.1,2.919))
rain(rainc = RNC, rainnc = RNNC)

Function to convert/save a SpatRaster array/Netcdf

Description

Conversion of SpatRaster to array and optionally save on a existing Netcdf File.

Usage

rast_to_netcdf(r, file, name, unit = units(r), XY = FALSE, verbose = TRUE)

Arguments

Value

numerical array

Note

eva3dm::wrf_rast support 3d SpatRaster, in case of a 4d variable use other approach to save on file.

Examples

folder   <- system.file("extdata",package="eva3dm")
wrf_file <- paste0(folder,"/wrf.day1.o3.nc")

Rast     <- wrf_rast(wrf_file,'o3')
A        <- rast_to_netcdf(Rast)

Function to read stats and evaluation

Description

Function to read stats and evaluation output

Usage

read_stat(file, sep = ";", dec = ".", verbose = FALSE, ...)

Arguments

Value

No return value

Examples

sample <- read_stat(file    = paste0(system.file("extdata", package = "eva3dm"),"/sample.txt"),
                    verbose = TRUE)

sample <- read_stat(file    = paste0(system.file("extdata", package = "eva3dm"),"/sample.csv"),
                    verbose = TRUE)

Convert relative humidity to absolute humidity

Description

function to convert humidity to absolute humidity using Tetens formula, assuming standard atmosphere conditions.

Usage

rh2q(rh, temp = 15)

Arguments

Value

value of data.frame with time and the absolute humidity, units are g/g

Note

default values are from standard atmosphere (288.15 K / 15 C)

if rh and temp arguments are data.frame, both need to have the same number of lines and columns, first column (time column) will be ignored.

Examples

# for a singfle value
rh2q(rh = 99, temp = 25)

# vector of rh values
rh2q(rh = c(0,seq(1,100, by = 4)), temp = 25)

# vector of values for rh and temp
rh2q(rh = c(0,seq(1,100, by = 4)), temp = 10:35)

# rh is data.frame and temp is a value
times <- seq(as.POSIXct('2024-01-01',tz = 'UTC'),
             as.POSIXct('2024-01-02',tz = 'UTC'),
             by = 'hour')
rh2q(rh   = data.frame(time = times, a = seq(1,100, by = 4)),temp = 25)

# using both rh and temp are data.frames
rh2q(rh   = data.frame(time = times, a = seq(1,100, by = 4)),
     temp = data.frame(time = times, a = 11:35))

Functions to model evaluation using satellite

Description

functions to evaluate the spatial performance using satellite

Usage

sat(
  mo,
  ob,
  rname,
  table = NULL,
  n = 6,
  min = NA,
  max = NA,
  scale,
  method = "bilinear",
  eval_function = stat,
  mask,
  skip_interp = FALSE,
  verbose = TRUE,
  ...
)

Arguments

Value

a data.frame

Note

If a YOU DIED error message appears, means you are removing all the valid values using the arguments min or max.

If cate() is used for eval_function, the argument threshold must be included (see example).

Examples

model_no2 <- terra::rast(paste0(system.file("extdata",package="eva3dm"),
                              "/camx_no2.Rds"))
omi_no2   <- terra::rast(paste0(system.file("extdata",package="eva3dm"),
                              "/omi_no2.Rds"))

# generate the statistical indexes
sat(mo = model_no2,ob = omi_no2,rname = 'NO2_statistical')

# generate categorical evaluation using 3.0 as threshold
sat(mo = model_no2,ob = omi_no2,rname = 'NO2_categorical',
    eval_function = cate, threshold = 3.0)

# customizing the evaluation function: inclusion of p.value from stats::cor.test()
stat_p <- function(x, y, ...){
  table         <- eva3dm::stat(x, y, ...)
  cor.result    <- stats::cor.test(x, y, ... )
  table$p.value <- cor.result$p.value
  table         <- table[,c(1:4,12,5:11)]
  return(table)
}

sat(mo = model_no2,ob = omi_no2,rname = 'NO2_statistical_with_p',eval_function = stat_p)

Selection from data.frames with time-series

Description

Utility function to select periods from a data.frame. This function is inspired by openair::selectByDate.

Usage

select(
  data,
  year,
  month,
  day,
  hour,
  minutes,
  seconds,
  start,
  end,
  range,
  time = "date"
)

Arguments

Value

data.frame

See Also

See %IN% for selection based on position and model domains.

Examples

model <- readRDS(paste0(system.file("extdata",package="eva3dm"),
                        "/model.Rds"))
summary(model)
summary(select(data = model, start = '2012-01-09'))
summary(select(data = model, start = '2012-01-05', end  = '2012-01-09'))
summary(select(data = model, day  = 6))
summary(select(data = model, hour = 12))
summary(select(data = model, day = 6, hour = 12))
summary(select(data = model, day  = 'weekday'))
summary(select(data = model, day  = 'weekend'))
summary(select(data = model, day  = 'tue'))
summary(select(data = model, day  = 'jan'))

Calculate evaluation statistics from numerical vectors

Description

Calculate statistical indexes (Number of pairs, observation average, model average, correlation, Index Of Agreement, Factor of 2, Root Mean Square Error, Mean Bias, Mean error, Normalized Mean Bias, and Normalized Mean Bias) for model evaluation

Usage

stat(
  model,
  observation,
  wd = FALSE,
  cutoff = NA,
  cutoff_NME = NA,
  nobs = 8,
  rname,
  verbose = TRUE
)

Arguments

Value

data.frame with calculated Number of pairs, observation average, model average, correlation, Index Of Agreement, Factor of 2, Root Mean Square Error, Mean Bias, Mean error, Normalized Mean Bias, and Normalized Mean Bias

Note

the option wd = TRUE applies a rotation of 360 on model wind direction to minimize the angular difference.

References

Emery, C. and Tai., E. 2001. Enhanced Meteorological Modeling and Performance Evaluation for Two Texas Ozone Episodes.

Monk, K. et al. 2019. Evaluation of Regional Air Quality Models over Sydney and Australia: Part 1—Meteorological Model Comparison. Atmosphere 10(7), p. 374. doi: 10.3390/atmos10070374.

Ramboll. 2018. PacWest Newport Meteorological Performance Evaluation.

Zhang, Y. et al. 2019. Multiscale Applications of Two Online-Coupled Meteorology-Chemistry Models during Recent Field Campaigns in Australia, Part I: Model Description and WRF/Chem-ROMS Evaluation Using Surface and Satellite Data and Sensitivity to Spatial Grid Resolutions. Atmosphere 10(4), p. 189. doi: 10.3390/atmos10040189.

Emery, C., Liu, Z., Russell, A.G., Odman, M.T., Yarwood, G. and Kumar, N. 2017. Recommendations on statistics and benchmarks to assess photochemical model performance. Journal of the Air & Waste Management Association 67(5), pp. 582–598. doi: 10.1080/10962247.2016.1265027.

Zhai, H., Huang, L., Emery, C., Zhang, X., Wang, Y., Yarwood, G., ... & Li, L. (2024). Recommendations on benchmarks for photochemical air quality model applications in China—NO2, SO2, CO and PM10. Atmospheric Environment, 319, 120290.

Examples

model <- 1:100
data  <- model + rnorm(100,0.2)
stat(model = model, observation = data)

Create templates for model evaluation

Description

Create templates of code (r-scripts and bash job-submission script) to read, post-process and evaluate model results.

Usage

template(
  root,
  template = "WRF",
  case = "case",
  env = "rspatial",
  scheduler = "SBATCH",
  partition = "main",
  project = "PROJECT",
  verbose = TRUE
)

Arguments

Value

no value returned, create folders and other template scripts

Note

Templates types available:
- WRF (model post-process for METAR + INMET)
- WRF-3 (model post-process for METAR + INMET for triple nested domains)
- WRF-Chem (model post-process for METAR, AQS in Brazil and AERONET)
- EXP (model post-process for one experimental site including PBL variables)
- CAMx (post-process for triple tested domains) - METAR (download METAR observations from ASOS)
- MET (evaluation of meteorology)
- MET-3 (evaluation of meteorology for triple nested domains)
- AQ (evaluation of air quality)
- PSA (model post-processing with CDO for satellite evaluation)
- SAT (evaluation of precipitation using GPCP satellite)
- AQS_BR (download data from air quality stations at Sao Paulo and Rio de Janeiro)
- INMET (pre-processing of automatic and conventional meteorological data from INMET)
- merge (merge INMET data and merge METAR data)

Examples

temp <- file.path(tempdir(),"POST")
template(root = temp,template = 'WRF', case = 'WRF-only')

Function to calculate model wind direction

Description

Function to calculate model wind direction

Usage

uv2wd(u, v, verbose = TRUE)

Arguments

Value

vector or data.frame with time and the wind direction, units are degree north

Examples

times <- seq(as.POSIXct('2024-01-01',tz = 'UTC'),
             as.POSIXct('2024-01-02',tz = 'UTC'),
             by = 'hour')
U10 = data.frame(times = times,
                 test1 = c(3.29,2.07,1.96,2.82,3.73,
                           4.11,4.96,6.33,7.39,7.59,
                           7.51,7.22,6.81,6.43,5.81,
                           4.02,3.03,2.68,2.40,2.20,
                           2.09,1.95,1.66,1.39,1.4),
                 test2 = c(6.29,4.87,6.16,7.12,8.77,
                           10.16,10.85,11.45,11.21,11.04,
                           11.09,10.67,10.48,10.00,8.96,
                           6.36,5.62,5.83,5.83,5.25,
                           4.11,3.08,2.26,1.14,-0.10))
V10 = data.frame(times = times,
                 test1 = c(-8.87,-4.23,-2.81,-2.59,-4.58,
                           -4.80,-5.33,-5.86,-6.12,-6.13,
                           -6.11,-5.76,-5.91,-5.60,-5.09,
                           -3.33,-2.50,-2.29,-2.14,-2.07,
                           -1.95,-1.97,-2.04,-2.03,-1.9),
                 test2 = c(11.80,5.88,5.74,5.56,6.87,
                           8.39,8.68,8.33,7.90,7.42,
                           6.96,6.87,6.36,5.61,5.16,
                           4.16,4.25,4.59,4.51,3.90,
                           2.97,1.98,1.04,-0.08,-0.44))

uv2wd(u = U10, v = V10)

Function to calculate model wind speed

Description

Function to calculate model wind speed

Usage

uv2ws(u, v, verbose = TRUE)

Arguments

Value

vector or data.frame with time and the wind sped, units are m/s

Examples

times <- seq(as.POSIXct('2024-01-01',tz = 'UTC'),
             as.POSIXct('2024-01-02',tz = 'UTC'),
             by = 'hour')

U10 = data.frame(times = times,
                 test1 = c(3.29,2.07,1.96,2.82,3.73,
                           4.11,4.96,6.33,7.39,7.59,
                           7.51,7.22,6.81,6.43,5.81,
                           4.02,3.03,2.68,2.40,2.20,
                           2.09,1.95,1.66,1.39,1.4),
                 test2 = c(6.29,4.87,6.16,7.12,8.77,
                           10.16,10.85,11.45,11.21,11.04,
                           11.09,10.67,10.48,10.00,8.96,
                           6.36,5.62,5.83,5.83,5.25,
                           4.11,3.08,2.26,1.14,-0.10))
V10 = data.frame(times = times,
                 test1 = c(-8.87,-4.23,-2.81,-2.59,-4.58,
                           -4.80,-5.33,-5.86,-6.12,-6.13,
                           -6.11,-5.76,-5.91,-5.60,-5.09,
                           -3.33,-2.50,-2.29,-2.14,-2.07,
                           -1.95,-1.97,-2.04,-2.03,-1.9),
                 test2 = c(11.80,5.88,5.74,5.56,6.87,
                           8.39,8.68,8.33,7.90,7.42,
                           6.96,6.87,6.36,5.61,5.16,
                           4.16,4.25,4.59,4.51,3.90,
                           2.97,1.98,1.04,-0.08,-0.44))

uv2ws(u = U10, v = V10)

Function to return variable names

Description

Return variable names of a NetCDF

Usage

vars(file = NA, action = "get", verbose = FALSE)

Arguments

Value

string

Examples

vars(paste0(system.file("extdata",package="eva3dm"),'/wrfinput_d01'))

Creates SpatRaster object from wrf file

Description

Creates a SpatRaster (terra R-package) object from a variable from wrf file (or another compatible NetCDF)

Usage

wrf_rast(
  file = file.choose(),
  name = NA,
  map,
  level = 1,
  times,
  latlon = FALSE,
  method = "bilinear",
  as_polygons = FALSE,
  flip_h = FALSE,
  flip_v = FALSE,
  verbose = FALSE,
  ...
)

Arguments

Value

SpatRaster object (terra package)

Examples

{

wrf <- paste(system.file("extdata", package = "eva3dm"),
                         "/wrfinput_d01", sep="")

r <- wrf_rast(file=wrf, name='XLAT')

plot_rast(r)
}

Functions to write stats and evaluation

Description

Functions to write the output from evaluation functions. If the file name ends with .csv the function write.csv is used otherwise the function write.table is used.

Usage

write_stat(stat, file, sep = ";", dec = ".", verbose = FALSE, ...)

Arguments

Value

No return value

Examples

sample <- read_stat(paste0(system.file("extdata", package = "eva3dm"),"/sample.csv"),
                    verbose = TRUE)
dir.create(file.path(tempdir(), "stats"))

write_stat(file    = paste0(file.path(tempdir(), "stats"),'/sample.txt'),
           stat    = sample,
           verbose = TRUE)

write_stat(file    = paste0(file.path(tempdir(), "stats"),'/sample.csv'),
           stat    = sample,
           verbose = TRUE)