R/compHLLMLflux.R
In steingod/R-osisaf: Useful functions for handling EUMETSAT Ocean and Sea Ice SAF data
#
# NAME:
# statcompare
#
# PURPOSE:
# To compare OSISAF Radiative Flux data produced at HL and LML. Statistics
# are generated by C-code compstat which in turn reads the input generated
# by complmlhl.
#
# This software reads all the files generated and produces global
# statistics and illustrations.
#
# REQUIREMENTS:
#
# INPUT:
#
# OUTPUT:
#
# NOTES:
#
# BUGS:
#
# AUTHOR:
# Øystein Godøy, met.no/FOU, 15.07.2004
#
# MODIFIED:
# NA
#
compHLLMLflux <- function(basepath,experiment) {
##file <- "/home/steingod/software/osisaf/support/complmlhl/stathourlyssi.txt"
file <- paste(basepath,"/","stathourlyssi-",experiment,".txt",sep="")
t <- scan(file,what=list(d1=NULL,d2=NULL,d3=NULL,cl=integer(0)),
skip=1,nlines=1)
dsh <- scan(file,
what=list(n=integer(0),
year=integer(0),month=integer(0),day=integer(0),
hour=integer(0),minute=integer(0),
bias=double(0),biasc=double(0),biaso=double(0),
sd=double(0),sdc=double(0),sdo=double(0),
rms=double(0),rmsc=double(0),rmso=double(0),
biasrel=double(0),biasrelc=double(0),biasrelo=double(0),
sdrel=double(0),sdrelc=double(0),sdrelo=double(0),
rmsrel=double(0),rmsrelc=double(0),rmsrelo=double(0),
mhlfl=double(0),mhlflc=double(0),mhlflo=double(0),
mlmlfl=double(0),mlmlflc=double(0),mlmlflo=double(0),
sdhlfl=double(0),sdhlflc=double(0),sdhlflo=double(0),
sdlmlfl=double(0),sdlmlflc=double(0),sdlmlflo=double(0),
nobs=integer(0),nobsc=integer(0),nobso=integer(0),
mbiasbbao=double(0),sdbiasbbao=double(0),
relbiasbbao=double(0),relsdbiasbbao=double(0),
mhlbbao=double(0),mlmlbbao=double(0)
),
skip=3,na.strings="-99.00")
##file <- "/home/steingod/software/osisaf/support/complmlhl/stathourlydli.txt"
file <- paste(basepath,"/","stathourlydli-",experiment,".txt",sep="")
ddh <- scan(file,
what=list(n=integer(0),
year=integer(0),month=integer(0),day=integer(0),
hour=integer(0),minute=integer(0),
bias=double(0),biasc=double(0),biaso=double(0),
sd=double(0),sdc=double(0),sdo=double(0),
rms=double(0),rmsc=double(0),rmso=double(0),
biasrel=double(0),biasrelc=double(0),biasrelo=double(0),
sdrel=double(0),sdrelc=double(0),sdrelo=double(0),
rmsrel=double(0),rmsrelc=double(0),rmsrelo=double(0),
mhlfl=double(0),mhlflc=double(0),mhlflo=double(0),
mlmlfl=double(0),mlmlflc=double(0),mlmlflo=double(0),
sdhlfl=double(0),sdhlflc=double(0),sdhlflo=double(0),
sdlmlfl=double(0),sdlmlflc=double(0),sdlmlflo=double(0),
nobs=integer(0),nobsc=integer(0),nobso=integer(0)),
skip=3,na.strings="-99.00")
##file <- "/home/steingod/software/osisaf/support/complmlhl/statdaylissi.txt"
file <- paste(basepath,"/","statdailyssi-",experiment,".txt",sep="")
dsd <- scan(file,
what=list(n=integer(0),
year=integer(0),month=integer(0),day=integer(0),
hour=integer(0),minute=integer(0),
bias=double(0),sd=double(0),rms=double(0),
biasrel=double(0),sdrel=double(0),rmsrel=double(0),
mhlfl=double(0),mlmlfl=double(0),nobs=integer(0)),
skip=3,na.strings="-99.00")
##file <- "/home/steingod/software/osisaf/support/complmlhl/statdaylidli.txt"
file <- paste(basepath,"/","statdailydli-",experiment,".txt",sep="")
ddd <- scan(file,
what=list(n=integer(0),
year=integer(0),month=integer(0),day=integer(0),
hour=integer(0),minute=integer(0),
bias=double(0),sd=double(0),rms=double(0),
biasrel=double(0),sdrel=double(0),rmsrel=double(0),
mhlfl=double(0),mlmlfl=double(0),nobs=integer(0)),
skip=3,na.strings="-99.00")
dshtime <- ISOdatetime(dsh$year,dsh$month,dsh$day,dsh$hour,dsh$minute,0)
ddhtime <- ISOdatetime(ddh$year,ddh$month,ddh$day,ddh$hour,ddh$minute,0)
dsdtime <- ISOdatetime(dsd$year,dsd$month,dsd$day,dsd$hour,dsd$minute,0)
dddtime <- ISOdatetime(ddd$year,ddd$month,ddd$day,ddd$hour,ddd$minute,0)
layout(c(1,2,3))
par(pty="m")
# Daily plots
plot(dsdtime,dsd$bias,ylab=expression(plain(W)/plain(m)^2),
ylim=c(min(c(dsd$bias,dsd$sd,ddd$bias,ddd$sd),na.rm=T),
max(c(dsd$bias,dsd$sd,ddd$bias,ddd$sd),na.rm=T)),pch=1,col=4)
points(dsdtime,dsd$sd,pch=2,col=4)
points(dddtime,ddd$bias,pch=1,col=2)
points(dddtime,ddd$sd,pch=2,col=2)
abline(h=0)
title("Daily HL-LML bias",
sub=paste("SSI (blue), DLI (red), mean (circle), SD (triangle),",
"CL>=",t$cl))
plot(dsdtime,dsd$biasrel,ylab="%",
ylim=c(min(c(dsd$biasrel,dsd$sdrel,ddd$biasrel,ddd$sdrel),na.rm=T),
max(c(dsd$biasrel,dsd$sdrel,ddd$biasrel,ddd$sdrel),na.rm=T)),
pch=1,col=4)
points(dsdtime,dsd$sdrel,pch=2,col=4)
points(dddtime,ddd$biasrel,pch=1,col=2)
points(dddtime,ddd$sdrel,pch=2,col=2)
abline(h=0)
title("Daily HL-LML bias in percentage of mean",
sub=paste("SSI (blue), DLI (red), mean (circle), SD (triangle),",
"CL>=",t$cl))
plot(dsdtime,dsd$mhlfl,ylab=expression(plain(W)/plain(m)^2),
ylim=c(min(c(dsd$mhlfl,dsd$mlmlfl,ddd$mhlfl,ddd$mlmlfl),na.rm=T),
max(c(dsd$mhlfl,dsd$mlmlfl,ddd$mhlfl,ddd$mlmlfl),na.rm=T)),pch=1,col=4)
points(dsdtime,dsd$mlmlfl,pch=2,col=4)
points(dddtime,ddd$mhlfl,pch=1,col=2)
points(dddtime,ddd$mlmlfl,pch=2,col=2)
title("Daily mean HL and LML fluxes",
sub=paste("SSI (blue), DLI (red), HL (circle), LML (triangle),",
"CL>=",t$cl))
ans <- readline("Do you want to print this to a file (y/n): ")
if (substr(ans, 1, 1) == "y") {
dev.print(dev=jpeg,filename="osisaf-ic-daily-%03d.jpg",
width=600,height=600,bg="white")
}
# Hourly plots
plot(dshtime,dsh$bias,ylab=expression(plain(W)/plain(m)^2),
ylim=c(min(c(dsh$bias,dsh$sd,ddh$bias,ddh$sd),na.rm=T),
max(c(dsh$bias,dsh$sd,ddh$bias,ddh$sd),na.rm=T)),pch=1,col=4)
points(dshtime,dsh$sd,pch=2,col=4)
points(ddhtime,ddh$bias,pch=1,col=2)
points(ddhtime,ddh$sd,pch=2,col=2)
abline(h=0)
title("Hourly HL-LML bias",
sub=paste("SSI (blue), DLI (red), mean (circle), SD (triangle),",
"CL>=",t$cl))
plot(dshtime,dsh$biasrel,ylab="%",
ylim=c(min(c(dsh$biasrel,dsh$sdrel,ddh$biasrel,ddh$sdrel)),
max(c(dsh$biasrel,dsh$sdrel,ddh$biasrel,ddh$sdrel))),pch=1,col=4)
points(dshtime,dsh$sdrel,pch=2,col=4)
points(ddhtime,ddh$biasrel,pch=1,col=2)
points(ddhtime,ddh$sdrel,pch=2,col=2)
abline(h=0)
title("Hourly HL-LML bias in percentage of mean",
sub=paste("SSI (blue), DLI (red), mean (circle), SD (triangle),",
"CL>=",t$cl))
plot(dshtime,dsh$mhlfl,ylab=expression(plain(W)/plain(m)^2),
ylim=c(min(c(dsh$mhlfl,dsh$mlmlfl,ddh$mhlfl,ddh$mlmlfl),na.rm=T),
max(c(dsh$mhlfl,dsh$mlmlfl,ddh$mhlfl,ddh$mlmlfl),na.rm=T)),pch=1,col=4)
points(dshtime,dsh$mlmlfl,pch=2,col=4)
points(ddhtime,ddh$mhlfl,pch=1,col=2)
points(ddhtime,ddh$mlmlfl,pch=2,col=2)
title("Hourly mean HL and LML fluxes",
sub=paste("SSI (blue), DLI (red), HL (circle), LML (triangle),",
"CL>=",t$cl))
ans <- readline("Do you want to print this to a file (y/n): ")
if (substr(ans, 1, 1) == "y") {
dev.print(dev=jpeg,filename="osisaf-ic-hourly-%03d.jpg",
width=600,height=600,bg="white")
}
# Hourly clear sky plots
plot(dshtime,dsh$biasc,ylab=expression(plain(W)/plain(m)^2),
ylim=c(min(c(dsh$biasc,dsh$sdc,ddh$biasc,ddh$sdc),na.rm=T),
max(c(dsh$biasc,dsh$sdc,ddh$biasc,ddh$sdc),na.rm=T)),
pch=1,col=4)
points(dshtime,dsh$sdc,pch=2,col=4)
points(ddhtime,ddh$biasc,pch=1,col=2)
points(ddhtime,ddh$sdc,pch=2,col=2)
abline(h=0)
title("Hourly clear sky HL-LML bias",
sub=paste("SSI (blue), DLI (red), mean (circle), SD (triangle),",
"CL>=",t$cl))
plot(dshtime,dsh$biasrelc,ylab="%",
ylim=c(min(c(dsh$biasrelc,dsh$sdrelc,ddh$biasrelc,ddh$sdrelc),na.rm=T),
max(c(dsh$biasrelc,dsh$sdrelc,ddh$biasrelc,ddh$sdrelc),na.rm=T)),
pch=1,col=4)
points(dshtime,dsh$sdrelc,pch=2,col=4)
points(ddhtime,ddh$biasrelc,pch=1,col=2)
points(ddhtime,ddh$sdrelc,pch=2,col=2)
abline(h=0)
title("Hourly clear sky HL-LML bias in percentage of mean",
sub=paste("SSI (blue), DLI (red), mean (circle), SD (triangle),",
"CL>=",t$cl))
plot(dshtime,dsh$mhlflc,ylab=expression(plain(W)/plain(m)^2),
ylim=c(min(c(dsh$mhlflc,dsh$mlmlflc,ddh$mhlflc,ddh$mlmlflc),na.rm=T),
max(c(dsh$mhlflc,dsh$mlmlflc,ddh$mhlflc,ddh$mlmlflc),na.rm=T)),
pch=1,col=4)
points(dshtime,dsh$mlmlflc,pch=2,col=4)
points(ddhtime,ddh$mhlflc,pch=1,col=2)
points(ddhtime,ddh$mlmlflc,pch=2,col=2)
title("Hourly clear sky mean HL and LML fluxes",
sub=paste("SSI (blue), DLI (red), HL (circle), LML (triangle),",
"CL>=",t$cl))
ans <- readline("Do you want to print this to a file (y/n): ")
if (substr(ans, 1, 1) == "y") {
dev.print(dev=jpeg,filename="osisaf-ic-hourly-clear-%03d.jpg",
width=600,height=600,bg="white")
}
# Hourly overcast plots
plot(dshtime,dsh$biaso,ylab=expression(plain(W)/plain(m)^2),
ylim=c(min(c(dsh$biaso,dsh$sdo,ddh$biaso,ddh$sdo),na.rm=T),
max(c(dsh$biaso,dsh$sdo,ddh$biaso,ddh$sdo),na.rm=T)),
pch=1,col=4)
points(dshtime,dsh$sdo,pch=2,col=4)
points(ddhtime,ddh$biaso,pch=1,col=2)
points(ddhtime,ddh$sdo,pch=2,col=2)
abline(h=0)
title("Hourly overcast HL-LML bias",
sub=paste("SSI (blue), DLI (red), mean (circle), SD (triangle),",
"CL>=",t$cl))
plot(dshtime,dsh$biasrelo,ylab="%",
ylim=c(min(c(dsh$biasrelo,dsh$sdrelo,ddh$biasrelo,ddh$sdrelo),na.rm=T),
max(c(dsh$biasrelo,dsh$sdrelo,ddh$biasrelo,ddh$sdrelo),na.rm=T)),
pch=1,col=4)
points(dshtime,dsh$sdrelc,pch=2,col=4)
points(ddhtime,ddh$biasrelc,pch=1,col=2)
points(ddhtime,ddh$sdrelc,pch=2,col=2)
abline(h=0)
title("Hourly overcast HL-LML bias in percentage of mean",
sub=paste("SSI (blue), DLI (red), mean (circle), SD (triangle),",
"CL>=",t$cl))
plot(dshtime,dsh$mhlflo,ylab=expression(plain(W)/plain(m)^2),
ylim=c(min(c(dsh$mhlflo,dsh$mlmlflo,ddh$mhlflo,ddh$mlmlflo),na.rm=T),
max(c(dsh$mhlflo,dsh$mlmlflo,ddh$mhlflo,ddh$mlmlflo),na.rm=T)),
pch=1,col=4)
points(dshtime,dsh$mlmlflo,pch=2,col=4)
points(ddhtime,ddh$mhlflo,pch=1,col=2)
points(ddhtime,ddh$mlmlflo,pch=2,col=2)
title("Hourly overcast mean HL and LML fluxes",
sub=paste("SSI (blue), DLI (red), HL (circle), LML (triangle),",
"CL>=",t$cl))
ans <- readline("Do you want to print this to a file (y/n): ")
if (substr(ans, 1, 1) == "y") {
dev.print(dev=jpeg,filename="osisaf-ic-hourly-overcast-%03d.jpg",
width=600,height=600,bg="white")
}
plot(dshtime,dsh$mbiasbbao,ylab="%",
ylim=c(min(c(dsh$mbiasbbao,dsh$sdbiasbbao),na.rm=T),
max(c(dsh$mbiasbbao,dsh$sdbiasbbao),na.rm=T)),
pch=1,col=4)
points(dshtime,dsh$sdbiasbbao,pch=2,col=4)
title("Hourly overcast HL-LML BBA bias",
sub=paste("mean (circle), SD (triangle),",
"CL>=",t$cl))
abline(h=0)
plot(dshtime,dsh$relbiasbbao,ylab="%",
ylim=c(min(c(dsh$relbiasbbao,dsh$relsdbiasbbao),na.rm=T),
max(c(dsh$relbiasbbao,dsh$relsdbiasbbao),na.rm=T)),
pch=1,col=4)
points(dshtime,dsh$relsdbiasbbao,pch=2,col=4)
title("Hourly overcast HL-LML BBA bias in percentage of mean",
sub=paste("mean (circle), SD (triangle),",
"CL>=",t$cl))
abline(h=0)
plot(dshtime,dsh$mhlbbao,ylab="%",
ylim=c(min(c(dsh$mhlbbao,dsh$mlmlbbao),na.rm=T),
max(c(dsh$mhlbbao,dsh$mlmlbbao),na.rm=T)),
pch=1,col=4)
points(dshtime,dsh$mlmlbbao,pch=2,col=4)
title("Hourly overcast mean HL and LML BBA",
sub=paste("HL (circle), LML (triangle),",
"CL>=",t$cl))
ans <- readline("Do you want to print this to a file (y/n): ")
if (substr(ans, 1, 1) == "y") {
dev.print(dev=jpeg,filename="osisaf-ic-hourly-bbao-%03d.jpg",
width=600,height=600,bg="white")
}
md <- matrix(c(
"Ptype","bias","sd","relbias","relsd","MHL","MLML","N",
"SSI_24",
prettyNum(mean(dsd$bias,na.rm=T),digits=4),
prettyNum(mean(dsd$sd,na.rm=T),digits=4),
prettyNum(mean(dsd$biasrel,na.rm=T),digits=4),
prettyNum(mean(dsd$sdrel,na.rm=T),digits=4),
prettyNum(mean(dsd$mhlfl,na.rm=T),digits=4),
prettyNum(mean(dsd$mlmlfl,na.rm=T),digits=4),
length(dsd$bias),
"DLI_24",
prettyNum(mean(ddd$bias,na.rm=T),digits=4),
prettyNum(mean(ddd$sd,na.rm=T),digits=4),
prettyNum(mean(ddd$biasrel,na.rm=T),digits=4),
prettyNum(mean(ddd$sdrel,na.rm=T),digits=4),
prettyNum(mean(ddd$mhlfl,na.rm=T),digits=4),
prettyNum(mean(ddd$mlmlfl,na.rm=T),digits=4),
length(ddd$bias)
),
ncol=8,byrow=T)
print(md,justify="right")
mh <- matrix(c(
"Ptype","bias","sd","relbias","relsd","MHL","MLML","N",
"SSI_P",
prettyNum(mean(dsh$bias,na.rm=T),digits=4),
prettyNum(mean(dsh$sd,na.rm=T),digits=4),
prettyNum(mean(dsh$biasrel,na.rm=T),digits=4),
prettyNum(mean(dsh$sdrel,na.rm=T),digits=4),
prettyNum(mean(dsh$mhlfl,na.rm=T),digits=4),
prettyNum(mean(dsh$mlmlfl,na.rm=T),digits=4),
length(dsh$bias),
"SSI_PC",
prettyNum(mean(dsh$biasc,na.rm=T),digits=4),
prettyNum(mean(dsh$sdc,na.rm=T),digits=4),
prettyNum(mean(dsh$biasrelc,na.rm=T),digits=4),
prettyNum(mean(dsh$sdrelc,na.rm=T),digits=4),
prettyNum(mean(dsh$mhlflc,na.rm=T),digits=4),
prettyNum(mean(dsh$mlmlflc,na.rm=T),digits=4),
length(dsh$biasc),
"SSI_PO",
prettyNum(mean(dsh$biaso,na.rm=T),digits=4),
prettyNum(mean(dsh$sdo,na.rm=T),digits=4),
prettyNum(mean(dsh$biasrelo,na.rm=T),digits=4),
prettyNum(mean(dsh$sdrelo,na.rm=T),digits=4),
prettyNum(mean(dsh$mhlflo,na.rm=T),digits=4),
prettyNum(mean(dsh$mlmlflo,na.rm=T),digits=4),
length(dsh$biaso),
"DLI_P",
prettyNum(mean(ddh$bias,na.rm=T),digits=4),
prettyNum(mean(ddh$sd,na.rm=T),digits=4),
prettyNum(mean(ddh$biasrel,na.rm=T),digits=4),
prettyNum(mean(ddh$sdrel,na.rm=T),digits=4),
prettyNum(mean(ddh$mhlfl,na.rm=T),digits=4),
prettyNum(mean(ddh$mlmlfl,na.rm=T),digits=4),
length(ddh$bias),
"DLI_PC",
prettyNum(mean(ddh$biasc,na.rm=T),digits=4),
prettyNum(mean(ddh$sdc,na.rm=T),digits=4),
prettyNum(mean(ddh$biasrelc,na.rm=T),digits=4),
prettyNum(mean(ddh$sdrelc,na.rm=T),digits=4),
prettyNum(mean(ddh$mhlflc,na.rm=T),digits=4),
prettyNum(mean(ddh$mlmlflc,na.rm=T),digits=4),
length(ddh$biasc),
"DLI_PO",
prettyNum(mean(ddh$biaso,na.rm=T),digits=4),
prettyNum(mean(ddh$sdo,na.rm=T),digits=4),
prettyNum(mean(ddh$biasrelo,na.rm=T),digits=4),
prettyNum(mean(ddh$sdrelo,na.rm=T),digits=4),
prettyNum(mean(ddh$mhlflo,na.rm=T),digits=4),
prettyNum(mean(ddh$mlmlflo,na.rm=T),digits=4),
length(ddh$biaso),
"BBA_PO",
prettyNum(mean(dsh$mbiasbbao,na.rm=T),digits=4),
prettyNum(mean(dsh$sdbiasbbao,na.rm=T),digits=4),
prettyNum(mean(dsh$relbiasbbao,na.rm=T),digits=4),
prettyNum(mean(dsh$relsdbiasbbao,na.rm=T),digits=4),
prettyNum(mean(dsh$mhlbbao,na.rm=T),digits=4),
prettyNum(mean(dsh$mlmlbbao,na.rm=T),digits=4),
length(dsh$biaso)
),
ncol=8,byrow=T)
print(mh,justify="right")
layout(c(1))
}
steingod/R-osisaf documentation built on May 30, 2019, 2:32 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#
# NAME:
# statcompare
#
# PURPOSE:
# To compare OSISAF Radiative Flux data produced at HL and LML. Statistics
# are generated by C-code compstat which in turn reads the input generated
# by complmlhl.
#
# This software reads all the files generated and produces global
# statistics and illustrations.
#
# REQUIREMENTS:
#
# INPUT:
#
# OUTPUT:
#
# NOTES:
#
# BUGS:
#
# AUTHOR:
# Øystein Godøy, met.no/FOU, 15.07.2004
#
# MODIFIED:
# NA
#
compHLLMLflux <- function(basepath,experiment) {
##file <- "/home/steingod/software/osisaf/support/complmlhl/stathourlyssi.txt"
file <- paste(basepath,"/","stathourlyssi-",experiment,".txt",sep="")
t <- scan(file,what=list(d1=NULL,d2=NULL,d3=NULL,cl=integer(0)),
skip=1,nlines=1)
dsh <- scan(file,
what=list(n=integer(0),
year=integer(0),month=integer(0),day=integer(0),
hour=integer(0),minute=integer(0),
bias=double(0),biasc=double(0),biaso=double(0),
sd=double(0),sdc=double(0),sdo=double(0),
rms=double(0),rmsc=double(0),rmso=double(0),
biasrel=double(0),biasrelc=double(0),biasrelo=double(0),
sdrel=double(0),sdrelc=double(0),sdrelo=double(0),
rmsrel=double(0),rmsrelc=double(0),rmsrelo=double(0),
mhlfl=double(0),mhlflc=double(0),mhlflo=double(0),
mlmlfl=double(0),mlmlflc=double(0),mlmlflo=double(0),
sdhlfl=double(0),sdhlflc=double(0),sdhlflo=double(0),
sdlmlfl=double(0),sdlmlflc=double(0),sdlmlflo=double(0),
nobs=integer(0),nobsc=integer(0),nobso=integer(0),
mbiasbbao=double(0),sdbiasbbao=double(0),
relbiasbbao=double(0),relsdbiasbbao=double(0),
mhlbbao=double(0),mlmlbbao=double(0)
),
skip=3,na.strings="-99.00")
##file <- "/home/steingod/software/osisaf/support/complmlhl/stathourlydli.txt"
file <- paste(basepath,"/","stathourlydli-",experiment,".txt",sep="")
ddh <- scan(file,
what=list(n=integer(0),
year=integer(0),month=integer(0),day=integer(0),
hour=integer(0),minute=integer(0),
bias=double(0),biasc=double(0),biaso=double(0),
sd=double(0),sdc=double(0),sdo=double(0),
rms=double(0),rmsc=double(0),rmso=double(0),
biasrel=double(0),biasrelc=double(0),biasrelo=double(0),
sdrel=double(0),sdrelc=double(0),sdrelo=double(0),
rmsrel=double(0),rmsrelc=double(0),rmsrelo=double(0),
mhlfl=double(0),mhlflc=double(0),mhlflo=double(0),
mlmlfl=double(0),mlmlflc=double(0),mlmlflo=double(0),
sdhlfl=double(0),sdhlflc=double(0),sdhlflo=double(0),
sdlmlfl=double(0),sdlmlflc=double(0),sdlmlflo=double(0),
nobs=integer(0),nobsc=integer(0),nobso=integer(0)),
skip=3,na.strings="-99.00")
##file <- "/home/steingod/software/osisaf/support/complmlhl/statdaylissi.txt"
file <- paste(basepath,"/","statdailyssi-",experiment,".txt",sep="")
dsd <- scan(file,
what=list(n=integer(0),
year=integer(0),month=integer(0),day=integer(0),
hour=integer(0),minute=integer(0),
bias=double(0),sd=double(0),rms=double(0),
biasrel=double(0),sdrel=double(0),rmsrel=double(0),
mhlfl=double(0),mlmlfl=double(0),nobs=integer(0)),
skip=3,na.strings="-99.00")
##file <- "/home/steingod/software/osisaf/support/complmlhl/statdaylidli.txt"
file <- paste(basepath,"/","statdailydli-",experiment,".txt",sep="")
ddd <- scan(file,
what=list(n=integer(0),
year=integer(0),month=integer(0),day=integer(0),
hour=integer(0),minute=integer(0),
bias=double(0),sd=double(0),rms=double(0),
biasrel=double(0),sdrel=double(0),rmsrel=double(0),
mhlfl=double(0),mlmlfl=double(0),nobs=integer(0)),
skip=3,na.strings="-99.00")
dshtime <- ISOdatetime(dsh$year,dsh$month,dsh$day,dsh$hour,dsh$minute,0)
ddhtime <- ISOdatetime(ddh$year,ddh$month,ddh$day,ddh$hour,ddh$minute,0)
dsdtime <- ISOdatetime(dsd$year,dsd$month,dsd$day,dsd$hour,dsd$minute,0)
dddtime <- ISOdatetime(ddd$year,ddd$month,ddd$day,ddd$hour,ddd$minute,0)
layout(c(1,2,3))
par(pty="m")
# Daily plots
plot(dsdtime,dsd$bias,ylab=expression(plain(W)/plain(m)^2),
ylim=c(min(c(dsd$bias,dsd$sd,ddd$bias,ddd$sd),na.rm=T),
max(c(dsd$bias,dsd$sd,ddd$bias,ddd$sd),na.rm=T)),pch=1,col=4)
points(dsdtime,dsd$sd,pch=2,col=4)
points(dddtime,ddd$bias,pch=1,col=2)
points(dddtime,ddd$sd,pch=2,col=2)
abline(h=0)
title("Daily HL-LML bias",
sub=paste("SSI (blue), DLI (red), mean (circle), SD (triangle),",
"CL>=",t$cl))
plot(dsdtime,dsd$biasrel,ylab="%",
ylim=c(min(c(dsd$biasrel,dsd$sdrel,ddd$biasrel,ddd$sdrel),na.rm=T),
max(c(dsd$biasrel,dsd$sdrel,ddd$biasrel,ddd$sdrel),na.rm=T)),
pch=1,col=4)
points(dsdtime,dsd$sdrel,pch=2,col=4)
points(dddtime,ddd$biasrel,pch=1,col=2)
points(dddtime,ddd$sdrel,pch=2,col=2)
abline(h=0)
title("Daily HL-LML bias in percentage of mean",
sub=paste("SSI (blue), DLI (red), mean (circle), SD (triangle),",
"CL>=",t$cl))
plot(dsdtime,dsd$mhlfl,ylab=expression(plain(W)/plain(m)^2),
ylim=c(min(c(dsd$mhlfl,dsd$mlmlfl,ddd$mhlfl,ddd$mlmlfl),na.rm=T),
max(c(dsd$mhlfl,dsd$mlmlfl,ddd$mhlfl,ddd$mlmlfl),na.rm=T)),pch=1,col=4)
points(dsdtime,dsd$mlmlfl,pch=2,col=4)
points(dddtime,ddd$mhlfl,pch=1,col=2)
points(dddtime,ddd$mlmlfl,pch=2,col=2)
title("Daily mean HL and LML fluxes",
sub=paste("SSI (blue), DLI (red), HL (circle), LML (triangle),",
"CL>=",t$cl))
ans <- readline("Do you want to print this to a file (y/n): ")
if (substr(ans, 1, 1) == "y") {
dev.print(dev=jpeg,filename="osisaf-ic-daily-%03d.jpg",
width=600,height=600,bg="white")
}
# Hourly plots
plot(dshtime,dsh$bias,ylab=expression(plain(W)/plain(m)^2),
ylim=c(min(c(dsh$bias,dsh$sd,ddh$bias,ddh$sd),na.rm=T),
max(c(dsh$bias,dsh$sd,ddh$bias,ddh$sd),na.rm=T)),pch=1,col=4)
points(dshtime,dsh$sd,pch=2,col=4)
points(ddhtime,ddh$bias,pch=1,col=2)
points(ddhtime,ddh$sd,pch=2,col=2)
abline(h=0)
title("Hourly HL-LML bias",
sub=paste("SSI (blue), DLI (red), mean (circle), SD (triangle),",
"CL>=",t$cl))
plot(dshtime,dsh$biasrel,ylab="%",
ylim=c(min(c(dsh$biasrel,dsh$sdrel,ddh$biasrel,ddh$sdrel)),
max(c(dsh$biasrel,dsh$sdrel,ddh$biasrel,ddh$sdrel))),pch=1,col=4)
points(dshtime,dsh$sdrel,pch=2,col=4)
points(ddhtime,ddh$biasrel,pch=1,col=2)
points(ddhtime,ddh$sdrel,pch=2,col=2)
abline(h=0)
title("Hourly HL-LML bias in percentage of mean",
sub=paste("SSI (blue), DLI (red), mean (circle), SD (triangle),",
"CL>=",t$cl))
plot(dshtime,dsh$mhlfl,ylab=expression(plain(W)/plain(m)^2),
ylim=c(min(c(dsh$mhlfl,dsh$mlmlfl,ddh$mhlfl,ddh$mlmlfl),na.rm=T),
max(c(dsh$mhlfl,dsh$mlmlfl,ddh$mhlfl,ddh$mlmlfl),na.rm=T)),pch=1,col=4)
points(dshtime,dsh$mlmlfl,pch=2,col=4)
points(ddhtime,ddh$mhlfl,pch=1,col=2)
points(ddhtime,ddh$mlmlfl,pch=2,col=2)
title("Hourly mean HL and LML fluxes",
sub=paste("SSI (blue), DLI (red), HL (circle), LML (triangle),",
"CL>=",t$cl))
ans <- readline("Do you want to print this to a file (y/n): ")
if (substr(ans, 1, 1) == "y") {
dev.print(dev=jpeg,filename="osisaf-ic-hourly-%03d.jpg",
width=600,height=600,bg="white")
}
# Hourly clear sky plots
plot(dshtime,dsh$biasc,ylab=expression(plain(W)/plain(m)^2),
ylim=c(min(c(dsh$biasc,dsh$sdc,ddh$biasc,ddh$sdc),na.rm=T),
max(c(dsh$biasc,dsh$sdc,ddh$biasc,ddh$sdc),na.rm=T)),
pch=1,col=4)
points(dshtime,dsh$sdc,pch=2,col=4)
points(ddhtime,ddh$biasc,pch=1,col=2)
points(ddhtime,ddh$sdc,pch=2,col=2)
abline(h=0)
title("Hourly clear sky HL-LML bias",
sub=paste("SSI (blue), DLI (red), mean (circle), SD (triangle),",
"CL>=",t$cl))
plot(dshtime,dsh$biasrelc,ylab="%",
ylim=c(min(c(dsh$biasrelc,dsh$sdrelc,ddh$biasrelc,ddh$sdrelc),na.rm=T),
max(c(dsh$biasrelc,dsh$sdrelc,ddh$biasrelc,ddh$sdrelc),na.rm=T)),
pch=1,col=4)
points(dshtime,dsh$sdrelc,pch=2,col=4)
points(ddhtime,ddh$biasrelc,pch=1,col=2)
points(ddhtime,ddh$sdrelc,pch=2,col=2)
abline(h=0)
title("Hourly clear sky HL-LML bias in percentage of mean",
sub=paste("SSI (blue), DLI (red), mean (circle), SD (triangle),",
"CL>=",t$cl))
plot(dshtime,dsh$mhlflc,ylab=expression(plain(W)/plain(m)^2),
ylim=c(min(c(dsh$mhlflc,dsh$mlmlflc,ddh$mhlflc,ddh$mlmlflc),na.rm=T),
max(c(dsh$mhlflc,dsh$mlmlflc,ddh$mhlflc,ddh$mlmlflc),na.rm=T)),
pch=1,col=4)
points(dshtime,dsh$mlmlflc,pch=2,col=4)
points(ddhtime,ddh$mhlflc,pch=1,col=2)
points(ddhtime,ddh$mlmlflc,pch=2,col=2)
title("Hourly clear sky mean HL and LML fluxes",
sub=paste("SSI (blue), DLI (red), HL (circle), LML (triangle),",
"CL>=",t$cl))
ans <- readline("Do you want to print this to a file (y/n): ")
if (substr(ans, 1, 1) == "y") {
dev.print(dev=jpeg,filename="osisaf-ic-hourly-clear-%03d.jpg",
width=600,height=600,bg="white")
}
# Hourly overcast plots
plot(dshtime,dsh$biaso,ylab=expression(plain(W)/plain(m)^2),
ylim=c(min(c(dsh$biaso,dsh$sdo,ddh$biaso,ddh$sdo),na.rm=T),
max(c(dsh$biaso,dsh$sdo,ddh$biaso,ddh$sdo),na.rm=T)),
pch=1,col=4)
points(dshtime,dsh$sdo,pch=2,col=4)
points(ddhtime,ddh$biaso,pch=1,col=2)
points(ddhtime,ddh$sdo,pch=2,col=2)
abline(h=0)
title("Hourly overcast HL-LML bias",
sub=paste("SSI (blue), DLI (red), mean (circle), SD (triangle),",
"CL>=",t$cl))
plot(dshtime,dsh$biasrelo,ylab="%",
ylim=c(min(c(dsh$biasrelo,dsh$sdrelo,ddh$biasrelo,ddh$sdrelo),na.rm=T),
max(c(dsh$biasrelo,dsh$sdrelo,ddh$biasrelo,ddh$sdrelo),na.rm=T)),
pch=1,col=4)
points(dshtime,dsh$sdrelc,pch=2,col=4)
points(ddhtime,ddh$biasrelc,pch=1,col=2)
points(ddhtime,ddh$sdrelc,pch=2,col=2)
abline(h=0)
title("Hourly overcast HL-LML bias in percentage of mean",
sub=paste("SSI (blue), DLI (red), mean (circle), SD (triangle),",
"CL>=",t$cl))
plot(dshtime,dsh$mhlflo,ylab=expression(plain(W)/plain(m)^2),
ylim=c(min(c(dsh$mhlflo,dsh$mlmlflo,ddh$mhlflo,ddh$mlmlflo),na.rm=T),
max(c(dsh$mhlflo,dsh$mlmlflo,ddh$mhlflo,ddh$mlmlflo),na.rm=T)),
pch=1,col=4)
points(dshtime,dsh$mlmlflo,pch=2,col=4)
points(ddhtime,ddh$mhlflo,pch=1,col=2)
points(ddhtime,ddh$mlmlflo,pch=2,col=2)
title("Hourly overcast mean HL and LML fluxes",
sub=paste("SSI (blue), DLI (red), HL (circle), LML (triangle),",
"CL>=",t$cl))
ans <- readline("Do you want to print this to a file (y/n): ")
if (substr(ans, 1, 1) == "y") {
dev.print(dev=jpeg,filename="osisaf-ic-hourly-overcast-%03d.jpg",
width=600,height=600,bg="white")
}
plot(dshtime,dsh$mbiasbbao,ylab="%",
ylim=c(min(c(dsh$mbiasbbao,dsh$sdbiasbbao),na.rm=T),
max(c(dsh$mbiasbbao,dsh$sdbiasbbao),na.rm=T)),
pch=1,col=4)
points(dshtime,dsh$sdbiasbbao,pch=2,col=4)
title("Hourly overcast HL-LML BBA bias",
sub=paste("mean (circle), SD (triangle),",
"CL>=",t$cl))
abline(h=0)
plot(dshtime,dsh$relbiasbbao,ylab="%",
ylim=c(min(c(dsh$relbiasbbao,dsh$relsdbiasbbao),na.rm=T),
max(c(dsh$relbiasbbao,dsh$relsdbiasbbao),na.rm=T)),
pch=1,col=4)
points(dshtime,dsh$relsdbiasbbao,pch=2,col=4)
title("Hourly overcast HL-LML BBA bias in percentage of mean",
sub=paste("mean (circle), SD (triangle),",
"CL>=",t$cl))
abline(h=0)
plot(dshtime,dsh$mhlbbao,ylab="%",
ylim=c(min(c(dsh$mhlbbao,dsh$mlmlbbao),na.rm=T),
max(c(dsh$mhlbbao,dsh$mlmlbbao),na.rm=T)),
pch=1,col=4)
points(dshtime,dsh$mlmlbbao,pch=2,col=4)
title("Hourly overcast mean HL and LML BBA",
sub=paste("HL (circle), LML (triangle),",
"CL>=",t$cl))
ans <- readline("Do you want to print this to a file (y/n): ")
if (substr(ans, 1, 1) == "y") {
dev.print(dev=jpeg,filename="osisaf-ic-hourly-bbao-%03d.jpg",
width=600,height=600,bg="white")
}
md <- matrix(c(
"Ptype","bias","sd","relbias","relsd","MHL","MLML","N",
"SSI_24",
prettyNum(mean(dsd$bias,na.rm=T),digits=4),
prettyNum(mean(dsd$sd,na.rm=T),digits=4),
prettyNum(mean(dsd$biasrel,na.rm=T),digits=4),
prettyNum(mean(dsd$sdrel,na.rm=T),digits=4),
prettyNum(mean(dsd$mhlfl,na.rm=T),digits=4),
prettyNum(mean(dsd$mlmlfl,na.rm=T),digits=4),
length(dsd$bias),
"DLI_24",
prettyNum(mean(ddd$bias,na.rm=T),digits=4),
prettyNum(mean(ddd$sd,na.rm=T),digits=4),
prettyNum(mean(ddd$biasrel,na.rm=T),digits=4),
prettyNum(mean(ddd$sdrel,na.rm=T),digits=4),
prettyNum(mean(ddd$mhlfl,na.rm=T),digits=4),
prettyNum(mean(ddd$mlmlfl,na.rm=T),digits=4),
length(ddd$bias)
),
ncol=8,byrow=T)
print(md,justify="right")
mh <- matrix(c(
"Ptype","bias","sd","relbias","relsd","MHL","MLML","N",
"SSI_P",
prettyNum(mean(dsh$bias,na.rm=T),digits=4),
prettyNum(mean(dsh$sd,na.rm=T),digits=4),
prettyNum(mean(dsh$biasrel,na.rm=T),digits=4),
prettyNum(mean(dsh$sdrel,na.rm=T),digits=4),
prettyNum(mean(dsh$mhlfl,na.rm=T),digits=4),
prettyNum(mean(dsh$mlmlfl,na.rm=T),digits=4),
length(dsh$bias),
"SSI_PC",
prettyNum(mean(dsh$biasc,na.rm=T),digits=4),
prettyNum(mean(dsh$sdc,na.rm=T),digits=4),
prettyNum(mean(dsh$biasrelc,na.rm=T),digits=4),
prettyNum(mean(dsh$sdrelc,na.rm=T),digits=4),
prettyNum(mean(dsh$mhlflc,na.rm=T),digits=4),
prettyNum(mean(dsh$mlmlflc,na.rm=T),digits=4),
length(dsh$biasc),
"SSI_PO",
prettyNum(mean(dsh$biaso,na.rm=T),digits=4),
prettyNum(mean(dsh$sdo,na.rm=T),digits=4),
prettyNum(mean(dsh$biasrelo,na.rm=T),digits=4),
prettyNum(mean(dsh$sdrelo,na.rm=T),digits=4),
prettyNum(mean(dsh$mhlflo,na.rm=T),digits=4),
prettyNum(mean(dsh$mlmlflo,na.rm=T),digits=4),
length(dsh$biaso),
"DLI_P",
prettyNum(mean(ddh$bias,na.rm=T),digits=4),
prettyNum(mean(ddh$sd,na.rm=T),digits=4),
prettyNum(mean(ddh$biasrel,na.rm=T),digits=4),
prettyNum(mean(ddh$sdrel,na.rm=T),digits=4),
prettyNum(mean(ddh$mhlfl,na.rm=T),digits=4),
prettyNum(mean(ddh$mlmlfl,na.rm=T),digits=4),
length(ddh$bias),
"DLI_PC",
prettyNum(mean(ddh$biasc,na.rm=T),digits=4),
prettyNum(mean(ddh$sdc,na.rm=T),digits=4),
prettyNum(mean(ddh$biasrelc,na.rm=T),digits=4),
prettyNum(mean(ddh$sdrelc,na.rm=T),digits=4),
prettyNum(mean(ddh$mhlflc,na.rm=T),digits=4),
prettyNum(mean(ddh$mlmlflc,na.rm=T),digits=4),
length(ddh$biasc),
"DLI_PO",
prettyNum(mean(ddh$biaso,na.rm=T),digits=4),
prettyNum(mean(ddh$sdo,na.rm=T),digits=4),
prettyNum(mean(ddh$biasrelo,na.rm=T),digits=4),
prettyNum(mean(ddh$sdrelo,na.rm=T),digits=4),
prettyNum(mean(ddh$mhlflo,na.rm=T),digits=4),
prettyNum(mean(ddh$mlmlflo,na.rm=T),digits=4),
length(ddh$biaso),
"BBA_PO",
prettyNum(mean(dsh$mbiasbbao,na.rm=T),digits=4),
prettyNum(mean(dsh$sdbiasbbao,na.rm=T),digits=4),
prettyNum(mean(dsh$relbiasbbao,na.rm=T),digits=4),
prettyNum(mean(dsh$relsdbiasbbao,na.rm=T),digits=4),
prettyNum(mean(dsh$mhlbbao,na.rm=T),digits=4),
prettyNum(mean(dsh$mlmlbbao,na.rm=T),digits=4),
length(dsh$biaso)
),
ncol=8,byrow=T)
print(mh,justify="right")
layout(c(1))
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.