| Title: | 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. |
| Authors: | Daniel Schuch [aut, cre]
|
| Maintainer: | Daniel Schuch <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.99.1 |
| Built: | 2024-11-25 19:31:13 UTC |
| Source: | CRAN |
combines the stats (from individual station evaluation) and site list in a SpatVector using row.names
stat %at% sitestat %at% site
stat |
data.frame with stats or other variable (containing row.names and other variables) |
site |
data.frame with site list (containing row.names, lat and lon) |
SpatVector (terra package)
sites<- readRDS(paste0(system.file("extdata",package="eva3dm"),"/sites_AQ_BR.Rds")) model<- readRDS(paste0(system.file("extdata",package="eva3dm"),"/model.Rds")) obs <- readRDS(paste0(system.file("extdata",package="eva3dm"),"/obs.Rds")) 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) print(stats) geo_stats <- stats %at% sites print(geo_stats)sites<- readRDS(paste0(system.file("extdata",package="eva3dm"),"/sites_AQ_BR.Rds")) model<- readRDS(paste0(system.file("extdata",package="eva3dm"),"/model.Rds")) obs <- readRDS(paste0(system.file("extdata",package="eva3dm"),"/obs.Rds")) 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) print(stats) geo_stats <- stats %at% sites print(geo_stats)
results of 'd01 in d02' style syntax
x %IN% yx %IN% y
x |
data.frame |
y |
data.frame or character string |
data.frame with common columns or a cropped SpatRaster
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
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')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')
Read and write metadata information of a NetCDF files
atr(file = NA, var = "?", att = NA, action = "get", value = NA, verbose = TRUE)atr(file = NA, var = "?", att = NA, action = "get", value = NA, verbose = TRUE)
file |
file name |
var |
variable name, 0 to global and "?" to show options |
att |
attribute names (NA for get all attnames) |
action |
"get" (default), "write" or "print" (return the value) of an attribute |
value |
value to write |
verbose |
display additional information |
string with the NetCDF attribute value
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')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 traditional statistics related to a threshold
cate( model, observation, threshold, cutoff = NA, nobs = 8, rname, to.plot = FALSE, col = "#4444bb", pch = 19, lty = 3, lcol = "#333333", lim, verbose = TRUE, ... )cate( model, observation, threshold, cutoff = NA, nobs = 8, rname, to.plot = FALSE, col = "#4444bb", pch = 19, lty = 3, lcol = "#333333", lim, verbose = TRUE, ... )
model |
numeric vector with paired model data |
observation |
numeric vector with paired observation data |
threshold |
reference value |
cutoff |
(optionally the maximum) valid value for observation |
nobs |
minimum number of observations |
rname |
row name |
to.plot |
TRUE to plot a scatter-plot |
col |
color for points |
pch |
pch of points |
lty |
lty of threshold lines |
lcol |
col of threshold lines |
lim |
limit for x and y |
verbose |
display additional information |
... |
arguments passed to plot |
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
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).
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)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)
function to calculate daily mean, min or max of a data.frame
daily( data, time = "date", var, stat = mean, min_offset = 0, hour_offset = 0, numerical = TRUE, verbose = TRUE )daily( data, time = "date", var, stat = mean, min_offset = 0, hour_offset = 0, numerical = TRUE, verbose = TRUE )
data |
data.frame with time column and variable columns to be processed |
time |
name of the time column (default is date) in POSIXct |
var |
name of the columns to be calculated |
stat |
function of the statistics to calculate (default is mean) |
min_offset |
minutes of observation from previous hour (default is 0) |
hour_offset |
hours of observation from previous day (default is 0) |
numerical |
TRUE (defoult) include only numerical columns |
verbose |
display additional information |
data.frame with time and the daily mean, min or max
sites <- c("OAKB") for(site in sites){ cat('downloading METAR from:',site,'...\n') DATA <- riem::riem_measures(station = sites, date_start = "2012-01-01", date_end = "2012-02-01") } 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)sites <- c("OAKB") for(site in sites){ cat('downloading METAR from:',site,'...\n') DATA <- riem::riem_measures(station = sites, date_start = "2012-01-01", date_end = "2012-02-01") } 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)
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.
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, time = "date", verbose = TRUE, ... )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, time = "date", verbose = TRUE, ... )
mo |
data.frame with model data |
ob |
data.frame with observation data |
rname |
row name of the output (default is site argument) |
table |
data.frame to append the results |
site |
name of the stations or "ALL" (default), see notes |
wd |
default is FALSE, see notes |
fair |
model data.frame (or list of names) to perform a fair comparison, see notes |
cutoff |
minimum (optionally the maximum) valid value for observation |
cutoff_NME |
minimum (optionally the maximum) valid value for observation for NME |
no_tz |
ignore tz from input (force GMT) |
nobs |
minimum number of valid observations, default is 8 |
eval_function |
evaluation function (default is stat) |
time |
name of the time column (containing time in POSIXct) |
verbose |
display additional information |
... |
arguments to be passing to stats and plot |
data.frame with statistical values from stat or cate functions.
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.
stat for additional information about the statistical evaluation and cate for categorical evaluation.
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)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)
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.
extract_max_8h( filelist, variable = "o3", field = "4d", prefix = "max_8h", units = "ug m-3", meta = TRUE, filename, verbose = TRUE )extract_max_8h( filelist, variable = "o3", field = "4d", prefix = "max_8h", units = "ug m-3", meta = TRUE, filename, verbose = TRUE )
filelist |
list of files to be read |
variable |
variable name |
field |
'4d' (default), '3d', '2d' or '2dz' see notes |
prefix |
to output file, default is serie |
units |
units on netcdf file (default is ug m-3), change to skip unit conversion |
meta |
use Times, XLONG and XLAT data (only works with 2d variable for file) |
filename |
name for the file, in this case prefix is not used |
verbose |
display additional information |
No return value
The field argument '4d' / '2dz' is used to read a 4d/3d variable droping the 3rd dimention (z).
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')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')
Read and calculate the mean value of a variable from a list of wrf output files.
extract_mean( filelist, variable = "o3", field = "4d", prefix = "mean", units = "ppmv", meta = TRUE, filename, verbose = TRUE )extract_mean( filelist, variable = "o3", field = "4d", prefix = "mean", units = "ppmv", meta = TRUE, filename, verbose = TRUE )
filelist |
list of files to be read |
variable |
variable name |
field |
'4d' (default), '3d', '2d' or '2dz' see notes |
prefix |
to output file, default is serie |
units |
units on netcdf file (default is ppmv) |
meta |
use Times, XLONG and XLAT data (only works with 2d variable for file) |
filename |
name for the file, in this case prefix is not used |
verbose |
display additional information |
No return value
The field argument '4d' / '2dz' is used to read a 4d/3d variable droping the 3rd dimention (z).
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'))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'))
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.
extract_serie( filelist, point, variable = "o3", field = "4d", prefix = "serie", new = "check", return.nearest = FALSE, fast = FALSE, use_ij = FALSE, latitude = "XLAT", longitude = "XLONG", use_TFLAG = FALSE, use_datesec = FALSE, id = "id", verbose = TRUE )extract_serie( filelist, point, variable = "o3", field = "4d", prefix = "serie", new = "check", return.nearest = FALSE, fast = FALSE, use_ij = FALSE, latitude = "XLAT", longitude = "XLONG", use_TFLAG = FALSE, use_datesec = FALSE, id = "id", verbose = TRUE )
filelist |
list of files to be read |
point |
data.frame with lat/lon |
variable |
variable name |
field |
'4d' (defoult), '3d', '2d' or '2dz' see notes |
prefix |
to output file, default is serie |
new |
TRUE, FALSE of 'check' see notes |
return.nearest |
return the data.frame of nearest points instead of extract the serie |
fast |
faster calculation of grid distances but less precise |
use_ij |
logical, use i and j from input instead of calculate |
latitude |
name of latitude coordinate variable in the netcdf |
longitude |
name of longitude coordinate variable in the netcdf |
use_TFLAG |
use the variable TFLAG (CMAQ / smoke) instead of Times (WRF) |
use_datesec |
use the variable date and datesec (WACCM / CAM-Chem) instead of Times (WRF) |
id |
name of the column with station names, if point is a SpatVector (points) from terra package |
verbose |
display additional information |
No return value
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).
cat('Example 1: METAR site list\n') sites <- readRDS(paste0(system.file("extdata",package="eva3dm"),"/sites_METAR.Rds")) cat('Example 2: AERONET site list\n') sites <- readRDS(paste0(system.file("extdata",package="eva3dm"),"/sites_AERONET.Rds")) cat('Example 3: list of INMET stations on Brazil\n') sites <- readRDS(paste0(system.file("extdata",package="eva3dm"),"/sites_INMET.Rds")) cat('Example 4: 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'))cat('Example 1: METAR site list\n') sites <- readRDS(paste0(system.file("extdata",package="eva3dm"),"/sites_METAR.Rds")) cat('Example 2: AERONET site list\n') sites <- readRDS(paste0(system.file("extdata",package="eva3dm"),"/sites_AERONET.Rds")) cat('Example 3: list of INMET stations on Brazil\n') sites <- readRDS(paste0(system.file("extdata",package="eva3dm"),"/sites_INMET.Rds")) cat('Example 4: 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
get_distances(lat1, long1, lat2, long2, R = 6371)get_distances(lat1, long1, lat2, long2, R = 6371)
lat1 |
Latitude in decimals |
long1 |
Longitude in decimals |
lat2 |
Latitude in decimals |
long2 |
Longitude in decimals |
R |
Radius of the earth in kmdescription (R=6371) |
A numeric vector with the distance in kilometers.
#' source: https://github.com/gustavobio/brclimate/blob/master/R/get_distances.R
function to calculate Ozone Maximum Daily 8-hr Average or 8-hr moving Average for a data.frame
hourly( data, time = "date", var, stat = mean, min_offset = 30, numerical = TRUE, verbose = TRUE )hourly( data, time = "date", var, stat = mean, min_offset = 30, numerical = TRUE, verbose = TRUE )
data |
data.frame with time column and variable columns to be processed |
time |
name of the time column (default is date) in POSIXct |
var |
name of the columns to be calculated |
stat |
function of the statistics to calculate (default is mean) |
min_offset |
minutes of observation from previous hour (default is 30) |
numerical |
TRUE (default) includes only numerical columns |
verbose |
display additional information |
data.frame including only numerical columns
data.frame with time and the hourly mean, min or max
sites <- c("OAHR") for(site in sites){ cat('downloading METAR from:',site,'...\n') DATA <- riem::riem_measures(station = sites, date_start = "2012-01-01", date_end = "2012-02-01") } 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)sites <- c("OAHR") for(site in sites){ cat('downloading METAR from:',site,'...\n') DATA <- riem::riem_measures(station = sites, date_start = "2012-01-01", date_end = "2012-02-01") } 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)
function to project and interpolate rast
interp(x, y, method = "bilinear", mask, verbose = FALSE)interp(x, y, method = "bilinear", mask, verbose = FALSE)
x |
rast to be interpolated |
y |
target rast of the interpolation |
method |
passed to terra::resample |
mask |
optional SpatVector to mask the results |
verbose |
display additional information (not used) |
SpatRaster (terra package)
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)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
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:", ... )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:", ... )
x |
rast or array |
y |
rast or array to mean (x is used only for the range in this case) |
text.width |
Longitude in decimals |
dig |
vector with number of digits for plot |
xjust |
passed to legend |
yjust |
passed to legend |
horiz |
passed to legend |
y.intersp |
passed to legend |
x.intersp |
passed to legend |
show.mean |
set TRUE to hide mean value |
unit |
a string for units |
label_mean |
label in case y is provided |
... |
extra arguments passed to legend |
No return value
for use with rast use before any change of projection
text.width can vary depending on map dimensions
x <- 1:10 + rnorm(10,sd = .4) plot(x,ty='l') legend_range(x)x <- 1:10 + rnorm(10,sd = .4) plot(x,ty='l') legend_range(x)
function to calculate Ozone 8-hour moving average for a data.frame
ma8h(data, time = "date", var, verbose = TRUE, ...)ma8h(data, time = "date", var, verbose = TRUE, ...)
data |
data.frame with time column and variable columns to be processed |
time |
name of the time column (default is date) in POSIXct |
var |
name of the columns to be calculated |
verbose |
display additional information |
... |
parameters passed to hourly |
data.frame with time and the 8-hour moving average
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'))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'))
function to calculate Ozone Maximum Daily 8-hr Average or 8-hr moving Average for a data.frame
mda8(data, time = "date", var, verbose = TRUE)mda8(data, time = "date", var, verbose = TRUE)
data |
data.frame with time column and variable columns to be processed |
time |
name of the time column (default is date) in POSIXct |
var |
name of the columns to be calculated |
verbose |
display additional information |
data.frame with time and the maximum daily 8-hr average
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'))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'))
Read a NetCDF and print the medatada
ncdump(file = file.choose())ncdump(file = file.choose())
file |
file name |
No return value, only display information
ncdump(file = paste0(system.file("extdata",package="eva3dm"), '/wrfinput_d01'))ncdump(file = paste0(system.file("extdata",package="eva3dm"), '/wrfinput_d01'))
Custon plot for SpatRaster (terra R-package) object based on terra package
overlay( p, z, col, lim = range(z, na.rm = TRUE), symmetry = TRUE, pch = 19, cex = 1, outside = TRUE, add = FALSE, plg = list(tic = "none", shrink = 1), pax = list(), expand = 1.15, ... )overlay( p, z, col, lim = range(z, na.rm = TRUE), symmetry = TRUE, pch = 19, cex = 1, outside = TRUE, add = FALSE, plg = list(tic = "none", shrink = 1), pax = list(), expand = 1.15, ... )
p |
SpatVector points |
z |
column name or a vector of values to plot |
col |
color |
lim |
range of values for scale |
symmetry |
calculate symmetrical scale |
pch |
type of point |
cex |
character expansion for the points |
outside |
to include values outside range |
add |
add to existing plot |
plg |
list of parameters passed to terra::add_legend |
pax |
list of parameters passed to graphics::axis |
expand |
to expand the plot region |
... |
arguments to be passing to terra::plot |
No return value
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')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')
Custom difference (x - y) plots for SpatRaster object (based on terra package)
plot_diff( x, y, col, absolute = TRUE, relative = TRUE, lim_1 = NA, lim_2 = NA, unit = c(units(x), expression("%")), ... )plot_diff( x, y, col, absolute = TRUE, relative = TRUE, lim_1 = NA, lim_2 = NA, unit = c(units(x), expression("%")), ... )
x |
SpatVector points |
y |
values to plot |
col |
color |
absolute |
to plot absolute difference |
relative |
to plot relative difference |
lim_1 |
range of values for scale |
lim_2 |
calculate symmetrical scale |
unit |
annotation for units |
... |
arguments to be passing to plot_raster |
No return value
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)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)
Custon plot for SpatRaster (terra R-package) object based on terra package
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, min = -3, max, unit, ... )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, min = -3, max, unit, ... )
r |
raster |
color |
color scale, or name of a custom color scale (see notes) |
ncolor |
number of colors |
proj |
TRUE to project the raster to lat-lon |
plg |
list of parameters passed to terra::add_legend |
pax |
list of parameters passed to graphics::axis |
latitude |
add a latitude axis |
longitude |
add a longitude axis |
int |
interval of latitude and longitude lines |
grid |
add grid (graticule style) |
grid_int |
interval of grid lines |
grid_col |
color for grid lines |
add_range |
add legend with max, average and min r values |
ndig |
number of digits for legend_range |
log |
TRUE to plot in log-scale |
range |
range of original values to plot |
min |
minimum log value for log scale (defoul is -3) |
max |
maximum log value for log scale |
unit |
title for color bar (defoult is ) |
... |
arguments to be passing to terra::plot |
No return value
color scale includes: 'eva3r' (default), 'eva4', 'blues' and 'diff'
wrf <- paste(system.file("extdata", package = "eva3dm"), "/wrfinput_d01", sep="") r <- wrf_rast(file=wrf, name='XLAT') plot_rast(r)wrf <- paste(system.file("extdata", package = "eva3dm"), "/wrfinput_d01", sep="") r <- wrf_rast(file=wrf, name='XLAT') plot_rast(r)
function to convert absolute humidity to relative humidity.
q2rh(q, t = 15, p = 101325)q2rh(q, t = 15, p = 101325)
q |
vector (or data.frame) of absolute humidity (in g/Kg) |
t |
vector (or data.frame) of temperature (in Celcius) |
p |
vector (or data.frame) of pressure (in Pa) |
vector or data.frame with time and the relative humidity, units are
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.
# 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)# 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)
function that converts model accumulated precipitation to hourly precipitation.
rain(rainc, rainnc, verbose = TRUE)rain(rainc, rainnc, verbose = TRUE)
rainc |
data.frame or SpatRaster with RAINC variable |
rainnc |
data.frame or SpatRaster with RAINNC variable |
verbose |
set TRUE to display additional information |
data.frame time and the hourly precipitation or SpatRaster hourly precipitation
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)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)
Conversion of SpatRaster to array and optionally save on a Netcdf File.
rast_to_netcdf(r, file, name, unit = units(r), XY = FALSE, verbose = TRUE)rast_to_netcdf(r, file, name, unit = units(r), XY = FALSE, verbose = TRUE)
r |
SpatRaster object |
file |
Netcdf file name |
name |
variable name on a Netcdf file |
unit |
unit of the variable (set to NA to don't change unit) |
XY |
set to true if MemoryOrder is XY (only if file is missing) |
verbose |
display additional information |
numerical array
eva3dm::wrf_rast support 3d SpatRaster, in case of a 4d variable use other approach to save on file.
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)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 output
read_stat(file, sep = ";", dec = ".", verbose = FALSE, ...)read_stat(file, sep = ";", dec = ".", verbose = FALSE, ...)
file |
model data.frame |
sep |
the field separator string, passed to read.table function |
dec |
he string to use for decimal points, passed to read.table function |
verbose |
display additional information |
... |
arguments passed to read.table functions |
No return value
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)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)
function to convert humidity to absolute humidity using Tetens formula, assuming standard atmosphere conditions.
rh2q(rh, temp = 15)rh2q(rh, temp = 15)
rh |
vector (or data.frame) of relative humidity (in percentage) |
temp |
vector (or data.frame) of temperature (in Celsius) |
value of data.frame with time and the absolute humidity, units are g/g
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.
# 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))# 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 evaluate the spatial performance using satellite
sat( mo, ob, rname, table = NULL, n = 6, min = NA, max = NA, method = "bilinear", eval_function = stat, mask, verbose = TRUE, ... )sat( mo, ob, rname, table = NULL, n = 6, min = NA, max = NA, method = "bilinear", eval_function = stat, mask, verbose = TRUE, ... )
mo |
SpatRaster or raster with model |
ob |
SpatRaster or raster with observations |
rname |
passed to stat |
table |
data.frame to append the results |
n |
number of points from the boundary removed, default is 5 |
min |
minimum value cutoff |
max |
maximum value cutoff |
method |
passed to terra::resample |
eval_function |
evaluation function (default is stat) |
mask |
optional SpatVector to mask the results |
verbose |
set TRUE to display additional information |
... |
other arguments passed to stat |
a data.frame
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).
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")) # generate the statistical indexes sat(mo = model_o3,ob = omi_o3,rname = 'NO2_statistical') # generate categorical evaluation using 3.0 as threshold sat(mo = model_o3,ob = omi_o3,rname = 'NO2_categorical', eval_function = cate, threshold = 3.0)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")) # generate the statistical indexes sat(mo = model_o3,ob = omi_o3,rname = 'NO2_statistical') # generate categorical evaluation using 3.0 as threshold sat(mo = model_o3,ob = omi_o3,rname = 'NO2_categorical', eval_function = cate, threshold = 3.0)
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
stat( model, observation, wd = FALSE, cutoff = NA, cutoff_NME = NA, nobs = 8, rname, verbose = TRUE )stat( model, observation, wd = FALSE, cutoff = NA, cutoff_NME = NA, nobs = 8, rname, verbose = TRUE )
model |
numeric vector with paired model data |
observation |
numeric vector with paired observation data |
wd |
logical, set true to apply a rotation on wind direction, see notes |
cutoff |
(optionally the maximum) valid value for observation |
cutoff_NME |
(optionally the maximum) valid value for observation for NME, MFB and MFE |
nobs |
minimum number of observations |
rname |
row name |
verbose |
display additional information |
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
the option wd = TRUE applies a rotation of 360 on model wind direction to minimize the angular difference.
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.
model <- 1:100 data <- model + rnorm(100,0.2) stat(model = model, observation = data)model <- 1:100 data <- model + rnorm(100,0.2) stat(model = model, observation = data)
Create templates of code (r-scripts and bash job-submission script) to read, post-process and evaluate model results.
template( root, template = "WRF", case = "case", env = "rspatial", scheduler = "SBATCH", partition = "main", project = "PROJECT", verbose = TRUE )template( root, template = "WRF", case = "case", env = "rspatial", scheduler = "SBATCH", partition = "main", project = "PROJECT", verbose = TRUE )
root |
directory to create the template |
template |
template type (see notes) |
case |
case to be evaluated |
env |
name of the conda environment |
scheduler |
job scheduler used (SBATCH or PBS) |
partition |
partition name |
project |
project name |
verbose |
display additional information |
no value returned, create folders and other template scripts
Templates types available:
- WRF (model post-process for METAR + INMET)
- WRF-Chem (model post-process for METAR, AQS in Brazil and AERONET)
- EXP (model post-process for one experimental site including PBL variables)
- METAR (download observations)
- MET (evaluation of meteorology)
- AQ (evaluation of air quality)
- PSA (model post-processing with CDO for satellite evaluation)
- SAT (evaluation of precipitation using GPCP satellite)
temp <- file.path(tempdir(),"POST") template(root = temp,template = 'WRF', case = 'WRF-only')temp <- file.path(tempdir(),"POST") template(root = temp,template = 'WRF', case = 'WRF-only')
Function to calculate model wind direction
uv2wd(u, v, verbose = TRUE)uv2wd(u, v, verbose = TRUE)
u |
data.frame with model time-series of U10 |
v |
data.frame with model time-series of V10 |
verbose |
display additional information |
vector or data.frame with time and the wind direction, units are degree north
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)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
uv2ws(u, v, verbose = TRUE)uv2ws(u, v, verbose = TRUE)
u |
data.frame with model time-series of U10 |
v |
data.frame with model time-series of V10 |
verbose |
display additional information |
vector or data.frame with time and the wind sped, units are m/s
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)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)
Return variable names of a NetCDF
vars(file = NA, action = "get", verbose = FALSE)vars(file = NA, action = "get", verbose = FALSE)
file |
file name |
action |
'get' to return variable names or 'print' to print |
verbose |
display additional information |
string
vars(paste0(system.file("extdata",package="eva3dm"),'/wrfinput_d01'))vars(paste0(system.file("extdata",package="eva3dm"),'/wrfinput_d01'))
Creates a SpatRaster (terra R-package) object from a variable from wrf file (or another compatible NetCDF)
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, ... )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, ... )
file |
wrf file |
name |
variable name |
map |
(optional) file with lat-lon variables and grid information |
level |
only for 4d data, numeric, default is 1 for surface (include all times) |
times |
only for 4d data, numeric, set to select time instead of levels (include all levels) |
latlon |
logical (default is FALSE), set TRUE project the output to "+proj=longlat +datum=WGS84 +no_defs" |
method |
method passed to terra::projection, default is bilinear |
as_polygons |
logical, true to return a SpatVector instead of SpatRaster |
flip_h |
horizontal flip (by rows) |
flip_v |
vertical flip (by cols) |
verbose |
display additional information |
... |
extra arguments passed to ncdf4::ncvar_get |
SpatRaster object (terra package)
{ wrf <- paste(system.file("extdata", package = "eva3dm"), "/wrfinput_d01", sep="") r <- wrf_rast(file=wrf, name='XLAT') plot_rast(r) }{ wrf <- paste(system.file("extdata", package = "eva3dm"), "/wrfinput_d01", sep="") r <- wrf_rast(file=wrf, name='XLAT') plot_rast(r) }
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.
write_stat(stat, file, sep = ";", dec = ".", verbose = FALSE, ...)write_stat(stat, file, sep = ";", dec = ".", verbose = FALSE, ...)
stat |
observed data.frame |
file |
model data.frame |
sep |
the field separator string, passed to write.table function |
dec |
he string to use for decimal points, passed to write.table function |
verbose |
display additional information |
... |
arguments passed to write.table and write.csv functions |
No return value
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)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)