filteroutreg.R 4.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128
  1. #!/usr/bin/env Rscript
  2. #
  3. # This script filters out from any 1d (time) to 3d (time-lat-lon) field the
  4. # effect of greenhouse gases, ENSO, or any mode of variability through a
  5. # multilinear regression over the corresponding index (greenhouse gas
  6. # concentration, ENSO index ..) at a set of lags you can choose.
  7. #
  8. # Usage : ./regressedts.R <ncfile> <ncvar> <indexfile> <indexvar> <outfile>
  9. # <lag1> <lag2> ... (as many lags as wanted)
  10. #
  11. # History : Virginie Guemas - Initial version - 2011
  12. # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  13. library(ncdf)
  14. args=commandArgs(TRUE)
  15. filein=args[1] # Input netcdf file
  16. varin=args[2] # Input netcdf variable name
  17. index=args[3] # Netcdf file containing the index which to filter the signal
  18. varindex=args[4] # Index variable name
  19. fileout=args[5] # Output file
  20. lags=c() # List of lags
  21. jind=6
  22. while ( jind <= length(args) ) {
  23. lags=c(lags,as.numeric(args[jind]))
  24. jind=jind+1
  25. }
  26. if ( is.null(lags) ) { lags=c(0) }
  27. fnc=open.ncdf(filein)
  28. absent=T
  29. for ( jvar in 1:fnc$nvars ) {
  30. if ( fnc$var[[jvar]]$name == varin ) {
  31. data=fnc$var[[jvar]]
  32. unit=data$units
  33. n0dims=data$ndims
  34. wdimsout=data$dim
  35. var=get.var.ncdf(fnc,varin)
  36. ndims=length(dim(var))
  37. miss=att.get.ncdf(fnc,varin,'_FillValue')
  38. if ( miss$hasatt == TRUE ) { var[var==miss$value]=NA }
  39. miss=att.get.ncdf(fnc,varin,'missing_value')
  40. if ( miss$hasatt == TRUE ) { var[var==miss$value]=NA }
  41. absent=F
  42. }
  43. }
  44. if (absent) { stop("Variable not in input file !") }
  45. close.ncdf(fnc)
  46. system(paste('cdo showmon ',filein,'> mons 2>/dev/null'))
  47. mons=read.table('mons')
  48. system(paste('cdo showyear ',filein,'> years 2>/dev/null'))
  49. years=read.table('years')
  50. system('rm -f mons years')
  51. if (is.null(ndims)) { ndims=1 }
  52. dimsvar=switch(ndims,length(var),dim(var),c(dim(var)[1]*dim(var)[2],dim(var)[3]))
  53. if ( ndims > 1 ) { tmpvar=t(array(var,dim=dimsvar)) } else { tmpvar=var }
  54. var_ts=ts(tmpvar,start=c(years[[1]],mons[[1]]),end=c(years[[length(years)]],mons[[length(mons)]]),frequency=12)
  55. fnc=open.ncdf(index)
  56. vind=get.var.ncdf(fnc,varindex)
  57. close.ncdf(fnc)
  58. system(paste('cdo showmon ',index,'> mons 2>/dev/null'))
  59. mons=read.table('mons')
  60. system(paste('cdo showyear ',index,'> years 2>/dev/null'))
  61. years=read.table('years')
  62. system('rm -f mons years')
  63. index_ts=ts(vind,start=c(years[[1]],mons[[1]]),end=c(years[[length(years)]],mons[[length(mons)]]),frequency=12)
  64. lstind=list()
  65. for (jlag in lags) { lstind=list(lstind,lag(index_ts,-jlag)) }
  66. def=TRUE
  67. if ( ndims > 1 ) {
  68. for (jpt in 1:dimsvar[1]) {
  69. varjpt=var_ts[,jpt]
  70. if (length(sort(varjpt)) > 1 ) {
  71. toto = varjpt
  72. f='dataf[,1] ~ '
  73. jind=2
  74. for (jlag in lags) {
  75. toto = ts.intersect(toto,lag(index_ts,-jlag), dframe=FALSE)
  76. f=paste(f,'dataf[,',as.character(jind),'] +',sep="")
  77. jind=jind+1
  78. }
  79. dataf=toto
  80. if (def==TRUE) {
  81. tmp=array(,dim=c(dimsvar[1],dim(toto)[1]))
  82. def=FALSE
  83. }
  84. f=substr(f,start=1,stop=(nchar(f)-1))
  85. dataf=toto
  86. lm.out=lm(as.formula(f),data=toto,na.action=NULL)
  87. #lm.out=lm(varjpt ~ index_ts,na.action=na.omit)
  88. #tmp[jpt,]=varjpt-(lm.out$coefficients[1]+index_ts*lm.out$coefficients[2])
  89. tmp[jpt,]=toto[,1]-lm.out$fitted.values
  90. }
  91. }
  92. }else{
  93. toto = var_ts
  94. f='dataf[,1] ~ '
  95. jind=2
  96. for (jlag in lags) {
  97. toto = ts.intersect(toto,lag(index_ts,-jlag), dframe=FALSE)
  98. f=paste(f,'dataf[,',as.character(jind),'] +',sep="")
  99. jind=jind+1
  100. }
  101. f=substr(f,start=1,stop=(nchar(f)-1))
  102. dataf=toto
  103. lm.out=lm(as.formula(f),data=dataf,na.action=NULL)
  104. tmp=toto[,1]-lm.out$fitted.values
  105. }
  106. tmp[is.na(tmp)]=1e20
  107. wtime=dim.def.ncdf("time",paste("months since ",as.character(start(toto)[1]),"-",as.character(start(toto)[2]),"-15 12:00:00",sep=""),seq(0,dim(toto)[1]-1))
  108. wdimsout[[n0dims]]=wtime
  109. dimsout=switch(ndims,length(tmp),dim(tmp),c(dim(var)[1],dim(var)[2],dim(tmp)[2]))
  110. wvarout=var.def.ncdf(varin,unit,wdimsout,1e20)
  111. f=create.ncdf(fileout,wvarout)
  112. put.var.ncdf(f,wvarout,array(tmp,dim=dimsout))
  113. close.ncdf(f)
  114. rm(list=ls())
  115. quit()