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$