fdbk_dt_interpolate: Bin a data.table column into user defined bins and replace it...

In rfxf/Rfdbk: Handling NetCDF feedback files

Description

Bin a data.table column into user defined bins and replace it with the bin center value. If breaks can be provided (e.g. no gaps between bins) try to use 'cut' instead.

Usage

fdbk_dt_interpolate(DT, varToInter = c("obs", "veri_data"),
  levelToInter = "plevel", interLevels = levels, varno = "varno")

Arguments

DT	data.table
varToBin	variable that should be binned (and will be replaced by the binned version)
mode	that will be used to defined the bin. Choices are "bin" or "level". In the first case the limtis of the bins have to be explicitly given in two vectors. The name given to the corresponding levels of the bin will be the mean of the lower and upper limit of the bin. In the second case a vector specifying the levels has to be given. The limits of the bins will be calculated by taking the mean between neighbouring levels. The two methods differ in the sense that the "bin" mode allow to have gaps between the bins, whereas the bins will be continuous in "level" mode. The "level" mode allow to have non-equally spaced levels without gaps between the bins, so that the level is not always at the center of the bin.
binLower	number/vector lower bins limits
binUpper	number/vector upper bins limits
levels	number/vector of levels on which the bins will be defined

Value

data.table with varToBin replaced by factorized mid-bin values (NA if variable falls in none of the bins)

Author(s)

Josue <[e-mail]>

Examples

# Example of linear interpolation based on an international standard atmosphere profile
require(ggplot2)
require(Rfdbk)
require(reshape2)
a1 = -6.5 # K/km standard atmosphere lapse rate, represents observations
a2 = -9 # K/km lapse rate obtained from a fictive model output
b1 = 288.15 # K standard atmosphere surface temperature
b2 = 295 # K surface temperature obtained from a fictive model output
Ho = 8.4 # km scale height
po = 1013.25 # standard atmosphere pressure in hPa
p = seq(250,1000,10) # pressure until the tropopause
T1 = a1*Ho*log(po/p)+b1 # Standard amtmosphere temperature profile
T2 =  a2*Ho*log(po/p)+b2 # Model output temperature profile 
Bias = T2-T1 # Bias = forecast - observation 

# Build a data table in feedback files format
obs = T1
veri_data = T2
veri_forecast_time = 24
veri_initial_date = 2015110900
time = -720
lat = 46.812
lon = 6.943
varno = 2
veri_model = "COSMO"
plevel = p
ident = 6610
levels = c(1000, 975, 950, 925, 900, 875, 850, 800, 750, 700, 650, 600, 550, 500, 450, 400, 350, 300, 250) # Levels on which to interpolate
DT = data.frame(obs,veri_data,veri_forecast_time,veri_initial_date,time,lat,lon,varno,veri_model,plevel,ident)
DT           = fdbk_dt_interpolate(DT,varToInter=c("obs","veri_data"), levelToInter = "plevel", interLevels = levels) # interpolate DT

data1 = melt(data.frame(T1,p),id="T1") # Data for the standard atmosphere temperature profile
data2 = melt(data.frame(T2=DT$obs,DT$plevel),id="T2") # Interpolation of data1 


plot = ggplot() + geom_point(data=data1,aes(x=T1,y=value,colour=variable)) + geom_point(data=data2,aes(x=T2,y=value,colour=variable))+scale_y_reverse()+
  xlab("T [K]") + ylab("pressure [hPa]")+  scale_colour_manual(name="Temperature",values=c("red","blue"),labels=c("Interpolation","Standard atmosphere"))+theme(legend.position = "top")
print(plot) # plot of the Standard atmosphere profile and its interpolation 

allscores = fdbk_dt_verif_continuous(DT,strat=c("varno","veri_model","plevel") ) # Data table with scores

data3 = melt(data.frame(Bias,p),id="Bias") # Bias calculated directly from the standard atmosphere and model output profiles
ME = allscores[allscores$scorename=="ME"]$scores # scores calculated with fdbk_dt_verif_continuous on interpolation levels
ME_levels = allscores[allscores$scorename=="ME"]$plevel # interpolation levels
data4 = melt(data.frame(ME,ME_levels),id="ME")
plot2 = ggplot() + geom_point(data=data3,aes(x=Bias,y=value,colour=variable)) + geom_point(data=data4,aes(x=ME,y=value,colour=variable))+scale_y_reverse()+
  xlab("T bias [K]") + ylab("pressure [hPa]")+scale_colour_manual(name="Bias",values=c("red","blue"),labels=c("Interpolation","Standard atmosphere"))+theme(legend.position = "top")
print(plot2) # plot of the bias calculated directly from the profiles and the bias from the interpolation

rfxf/Rfdbk documentation built on May 27, 2019, 7:22 a.m.