R Notes

kvalobs@pak:~/NETCDF$ cat R_runs 
Rscript plot_netcdf.R "Inverse Distance Square" "Paramid: 110 Typeid: All" IDSALL IDS_ALL.nc 
Rscript plot_netcdf.R "Inverse Distance Square <= 50 km" "Paramid: 110 Typeid: All" IDS50 IDS50.nc 
Rscript plot_netcdf.R "Inverse Distance Square <= 100 km" "Paramid: 110 Typeid: All" IDS100  IDS100.nc
Rscript plot_netcdf.R "Linear Interpolation of Triangulated Neighbours" "Paramid: 110 Typeid: All" TRI50 TRI50.nc
Rscript plot_netcdf.R "IDS <= 50 km & Height Correction" "Paramid: 110 Typeid: All" IDS50Ht IDS50Ht.nc

kvalobs@pak:~/NETCDF$ cat plot_netcdf.R 
args <- commandArgs(TRUE)
Method=args[1]
Extra=args[2]
Shortname=args[3]
FileNC=args[4]
library(ncdf)
nd = open.ncdf(FileNC)
x=get.var.ncdf(nd,"Original")
y=get.var.ncdf(nd,"Interpolations")
#y=rnorm(length(x),0,1)+x 
#y=rnorm(length(x),0,20)+x   ... for testing the statistics
#y=1.5*x
x[x==-1]=0
v <- array(x, length(x))
w <- array(y, length(y))
xx=x[!is.na(x) & x>=0 & x < 100]
yy=y[!is.na(x) & x>=0 & x < 100]
yyy=yy[!is.na(yy) & yy>=0 & yy < 100]
xxx=xx[!is.na(yy) & yy>=0 & yy < 100]
jpeg(filename=paste(Shortname,'jpg',sep='.'))
plot(xxx,yyy,pch=46,xlab='Original Value',ylab='Interpoltion',main=Method,sub=Extra,xlim=c(0,100),ylim=c(0,100))
XC=cor(xxx,yyy)
#ff <- lm(xxx ~ yyy)
ff <- lm(yyy ~ xxx)
ff$coefficients[1]
ff$coefficients[2]
z=ff$coefficients[2]*xxx + ff$coefficients[1]
abline(ff$coefficients)
text(0,95,"Linear Regression",pos=4,cex=0.8)
text(0,90,paste("Gradient =",format(ff$coefficients[2],digits=5)),pos=4,cex=0.8)
text(0,85,paste("Intercept =",format(ff$coefficients[1],digits=5)),pos=4,cex=0.8)
text(0,80,paste("Correlation coefficient =",format(XC,digits=5)),pos=4,cex=0.8)

meany=mean(yyy)
chi2=0
for (i in 1:length(yyy)) {
    OE=(yyy[i]-xxx[i])*(yyy[i]-xxx[i])
    Vy=(yyy[i]-meany)*(yyy[i]-meany)
    chi2=OE/Vy + chi2
  }
chi2_red=chi2/length(yyy)
text(0,70,"Comparison to y~x model",pos=4,cex=0.8)
text(0,65,expression("Reduced "*chi^2*" = "),pos=4,cex=0.8)
text(30,65,format(round(chi2_red, 1), nsmall = 1),pos=4,cex=0.8)
text(0,60,paste("N Points = ",length(yyy)),pos=4,cex=0.8)

dev.off()


kvalobs@pak:~/NETCDF$ 



h2d <- gplots::hist2d(xxx,yyy,show=FALSE, same.scale=TRUE, nbins=c(50,50))
contour( h2d$x, h2d$y, h2d$counts, nlevels=1000 )
lines(xxx,z)

yyyy=yy[!is.na(yy) & yy>=5 & yy < 100]
xxxx=xx[!is.na(yy) & yy>=5 & yy < 100]
yyyyy=yyyy[!is.na(xxxx) & xxxx>=5 & xxxx < 100]
xxxxx=xxxx[!is.na(xxxx) & xxxx>=5 & xxxx < 100]
XXC=cor(xxxxx,yyyyy)
fff <- lm(xxxxx ~ yyyyy)
fff$coefficients[1]
fff$coefficients[2]
zzz=fff$coefficients[1]*xxxxx + fff$coefficients[2]

h2d <- gplots::hist2d(xxxxx,yyyyy,show=FALSE, same.scale=TRUE, nbins=c(50,50))
contour( h2d$x, h2d$y, h2d$counts, nlevels=1000 )
lines(xxxxx,zzz)
plot(xxxxx,yyyyy,pch=46,xlab='Original Value',ylab='Interpoltion',main="...",sub="...")

kvalobs@pak:~/NETCDF$