Scripts and data for Lindzen and Choi 2011 are now online at CA here together with the original article.
It will take me a little while to get to this. The scripts are in IDL. Translations to R welcomed.
Scripts and data for Lindzen and Choi 2011 are now online at CA here together with the original article.
It will take me a little while to get to this. The scripts are in IDL. Translations to R welcomed.
53 Comments
There is also the free version, GDL:-
http://en.wikipedia.org/wiki/GNU_Data_Language
michel,
regrettably, IDL implements more functions than GDL currently does. “JULDAY” is probably only one of them.
GDL> main
% Compiled module: MAIN.
% Compiled module: OBS.
% Compiled module: FILE_LINES.
% Compiled module: STRSPLIT.
% Compiled module: UNIQ.
% Compiled module: MEAN.
% OBS: Function not found: JULDAY
% Execution halted at: OBS 45 obs.pro
% MAIN 15 main.pro
% $MAIN$
GDL>
You can get the source code for JULDAY at:
http://idlastro.gsfc.nasa.gov/ftp/ittvislib/julday.pro
Thanks for filling that gap. However the next whole is TS_SMOOTH. I guess it won’t be the last. GDL is clearly less well equipt in terms of its library of fucntions at this stage and probably would require more coding effort than porting to R.
Unfortunately neither approach allows identical replication of L&C but replication by similar means in another language should confirm the same result and if there’s a significant difference (not due to porting errors) that would be significant as well.
P. Solar,
I got the L&C-11 scripts running with gdl version 0.9.1 and closely reproduced at least “Fig.jpg”. Here is a list of urls for all needed IDL library routines:
h ttp://idlastro.gsfc.nasa.gov/ftp/ittvislib/julday.pro
h ttp://idlastro.gsfc.nasa.gov/ftp/ittvislib/label_date.pro
h ttp://idlastro.gsfc.nasa.gov/ftp/pro/plot/legend.pro
h ttp://idlastro.gsfc.nasa.gov/ftp/ittvislib/regress.pro
h ttp://idlastro.gsfc.nasa.gov/ftp/ittvislib/timegen.pro
h ttp://idlastro.gsfc.nasa.gov/ftp/ittvislib/ts_coef.pro
h ttp://idlastro.gsfc.nasa.gov/ftp/ittvislib/ts_fcast.pro
h ttp://idlastro.gsfc.nasa.gov/ftp/ittvislib/ts_smooth.pro
h ttp://astro.uni-tuebingen.de/software/idl/textoidl/textoidl.pro
h ttp://coco.atmos.washington.edu/ion_script/gamap2/graphics/screen2jpg.pro
h ttp://idl-coyote.googlecode.com/svn-history/r293/trunk/coyote/retired/tvread.pro
Put the above files into the L&C-2011 scripts directory; cd into this directory and put the patch below into a file named, say, “LC-11.diff”. Then apply the lc11.diff using the shell command
patch -p3 < LC-11.diff
One last thing: My gdl version behaves very strangely (spurious syntax error messages) when gdl was called with main.pro on the command line. This errors did not show up when “main” was called from the gdl prompt.
diff -ur a/lindzen_choi_2011/scripts/amip.pro b/lindzen_choi_2011/scripts/amip.pro
--- a/lindzen_choi_2011/scripts/amip.pro 2011-09-25 18:40:42.000000000 +0200
+++ b/lindzen_choi_2011/scripts/amip.pro 2011-09-27 20:06:10.000000000 +0200
@@ -3,7 +3,7 @@
ObsName='ERBE&CERES'
clim_period=59
tot_period=288
-path=['./amip_olr/','./amip_osr/']
+path=['../amip_olr/','../amip_osr/']
file=['ncar_ccsm3_0_20n20s.dat',$
'mpi_echam51_20n20s.dat',$
'iap_fgoals1_0_g1_20n20s.dat',$
diff -ur a/lindzen_choi_2011/scripts/legend.pro b/lindzen_choi_2011/scripts/legend.pro
--- a/lindzen_choi_2011/scripts/legend.pro 2011-08-21 16:37:24.000000000 +0200
+++ b/lindzen_choi_2011/scripts/legend.pro 2011-09-27 21:09:00.000000000 +0200
@@ -442,7 +442,7 @@
if vectorfont[i] ne '' then begin
; if (num eq 1) and vertical then xp = x + xt/2 ; IF 1, CENTERED.
xyouts,xp,yp,vectorfont[i],width=width,color=colors[i] $
- ,charsize=charsize,align=xalign,charthick = charthick,/norm,font=font
+ ,charsize=charsize,align=xalign,charthick = charthick,/norm;,font=font
xt = xt > width
xp = xp + width/2.
endif else begin
@@ -456,11 +456,11 @@
TEXT_ONLY:
if vertical and (vectorfont[i] eq '') and symline and (linestyle[i] eq -99) then x=x0 + xspacing
xyouts,x,y,delimiter,width=width,/norm,color=textcolors[i], $
- charsize=charsize,align=xalign,charthick = charthick,font=font
+ charsize=charsize,align=xalign,charthick = charthick;,font=font
x = x + width*xsign
if width ne 0 then x = x + 0.5*xspacing
xyouts,x,y,items[i],width=width,/norm,color=textcolors[i],charsize=charsize, $
- align=xalign,charthick=charthick,font=font
+ align=xalign,charthick=charthick;,font=font
x = x + width*xsign
if not vertical and (i lt (n-1)) then x = x+2*xspacing; ADD INTER-ITEM SPACE
xfinal = (x + xspacing*margin)
diff -ur a/lindzen_choi_2011/scripts/main.pro b/lindzen_choi_2011/scripts/main.pro
--- a/lindzen_choi_2011/scripts/main.pro 2011-09-25 18:40:40.000000000 +0200
+++ b/lindzen_choi_2011/scripts/main.pro 2011-09-27 22:50:22.000000000 +0200
@@ -1,4 +1,5 @@
pro main
+myctblack=0
; -----------------------------------------------------------------------------
; This is the main program calling all procedures, calculating feedbacks, and
@@ -24,7 +25,7 @@
; endfor
; ---------------
-; Display 1 - Timeseries
+print,"Display 1 - Timeseries"
; ---------------
window,1,xs=600,ys=500
loadct,5
@@ -32,19 +33,19 @@
!p.charsize=1.3
!P.background=255
dummy=label_date(date_format=('%y'))
- plot,jul,sst,color=!myct.black,$
+ plot,jul,sst,color=myctblack,$
xtitle='Year',ytitle='T (K)',$
xtickunits=['time'],xtickformat='label_date',xstyle=1,xminor=5
- plot,jul,replicate(0,258),color=!myct.black,yr=[-4,8],$
+ plot,jul,replicate(0,258),color=myctblack,yr=[-4,8],$
xtitle='Year',ytitle=textoidl('R (Wm^{-2})'),$
xtickunits=['time'],xtickformat='label_date',xstyle=1,xminor=5
oplot,jul,olr,color=112 ;red
oplot,jul,osr,color=47 ;blue
- legend,label=['LW','SW'],line=[0,0],lcolor=[112,47],$
- charsize=1.3,boxcolor=255,spacing=1.3
+ legend,['LW','SW'],line=[0,0],colors=[112,47],$
+ charsize=1.3,outline_color=255,spacing=1.3
; ---------------
-; Feedback analysis for ERBE and CERES;
+print, "Feedback analysis for ERBE and CERES"
; Call lc11.pro procedure, which include both LC11 method and simple regression.
; ---------------
; Call procedure for calculating climate feedback
@@ -55,8 +56,12 @@
LagDim,SlopDim,ConstDim,CorrDim,SlopDim2,ConstDim2,CorrDim2
; ---------------
-; Display 2 - Timeseries
- goto,jump
+print, "Display 2 - Timeseries"
+; goto,jump
+; goto not available in GDL outside of procedures, hence:
+if 0 ne 0 then begin
+
+
; ---------------
; SST timeseries graph with warming and cooling intervals
window,2,xs=600,ys=300
@@ -78,7 +83,7 @@
endfor
;
; ---------------
-; Display 3 - Feedback analysis for ERBE and CERES
+print,"Display 3 - Feedback analysis for ERBE and CERES"
; ---------------
window,3,xs=600,ys=600
!p.multi=[0,2,2,0,0]
@@ -90,12 +95,12 @@
xtitle='T (K)',ytitle='R (W/m2)' ;textoidl('R (Wm^{-2})')
oplot,findgen(10)-5,SlopDim2(5)*(findgen(10)-5)+ConstDim2(5)
if ConstDim2(5) ge 0 then begin
- legend,label=['y='+string(SlopDim2(5),format='(f4.1)')+'x+'$
- +strtrim(string(ConstDim2(5),format='(f4.1)'),1)],boxcolor=-5,charsize=1.3
+ legend,['y='+string(SlopDim2(5),format='(f4.1)')+'x+'$
+ +strtrim(string(ConstDim2(5),format='(f4.1)'),1)],outline_color=-5,charsize=1.3
endif
if ConstDim2(5) lt 0 then begin
- legend,label=['y='+string(SlopDim2(5),format='(f4.1)')+'x'$
- +strtrim(string(ConstDim2(5),format='(f4.1)'),1)],boxcolor=-5,charsize=1.3
+ legend,['y='+string(SlopDim2(5),format='(f4.1)')+'x'$
+ +strtrim(string(ConstDim2(5),format='(f4.1)'),1)],outline_color=-5,charsize=1.3
endif
f_sim= SlopDim2(5) ; Feedback from simple regression
; Scatter plots for Lindzen-Choi Method
@@ -105,14 +110,14 @@
xtitle='T (K)',ytitle='R (W/m2)';textoidl('R (Wm^{-2})')
oplot,findgen(10)-5,SlopDim(pos(0))*(findgen(10)-5)+ConstDim(pos(0))
if ConstDim(pos(0)) ge 0 then begin
- legend,label=['y='+string(SlopDim(pos(0)),format='(f4.1)')+'x+'$
+ legend,['y='+string(SlopDim(pos(0)),format='(f4.1)')+'x+'$
+strtrim(string(ConstDim(pos(0)),format='(f4.1)'),1)],$
- boxcolor=-5,charsize=1.3
+ outline_color=-5,charsize=1.3
endif
if ConstDim(pos(0)) lt 0 then begin
- legend,label=['y='+string(SlopDim(pos(0)),format='(f4.1)')+'x'$
+ legend,['y='+string(SlopDim(pos(0)),format='(f4.1)')+'x'$
+strtrim(string(ConstDim(pos(0)),format='(f4.1)'),1)],$
- boxcolor=-5,charsize=1.3
+ outline_color=-5,charsize=1.3
endif
f_lc=SlopDim(pos(0)) ; Feedback from LC11
; Lead-lag slopes
@@ -120,26 +125,28 @@
xtitle='Lagged month',xminor=2, $
ytitle='R/T (W/m2/K)';textoidl('R/T (Wm^{-2}K^{-1})')
oplot,LagDim,SlopDim2,linestyle=1
- legend,label=['LC','Simple'],line=[0,1],boxcolor=-5,charsize=1.3,spacing=2
+ legend,['LC','Simple'],line=[0,1],outline_color=-5,charsize=1.3,spacing=2
; Lead-lag correlations
plot,LagDim,CorrDim,xr=[-5,5],xs=1,yr=[-1,1],$
title='Lag Correlation',xtitle='Lagged month',xminor=2
oplot,LagDim,CorrDim2,linestyle=1
- legend,label=['LC','Simple'],line=[0,1],boxcolor=-5,charsize=1.3
+ legend,['LC','Simple'],line=[0,1],outline_color=-5,charsize=1.3
- jump: print,'Skip Display 2 and 3'
+;jump:
+ endif
+ print,'Skip Display 2 and 3'
; ---------------
-; Open AMIP data
+print, "Open AMIP data"
; Call amip.pro procedure, which include 11 AMIP model outputs
; ---------------
-; Call proceedure for reading AMIP model outputs
+print,"Call proceedure for reading AMIP model outputs"
amip,olr_amip,osr_amip
-; Generate julian day corresponding to AMIP outputs
+print,"Generate julian day corresponding to AMIP outputs"
julmo=timegen(start=julday(1,1,1985),final=julday(12,1,2008),units='M')
; ---------------
-; Display 4 - Timeseries for AMIP data
+print, "Display 4 - Timeseries for AMIP data -1"
; ---------------
title=['CCSM3','ECHAM5/MIP-OM','FGOALS-g1.0',$
'GFDL-CM2.1','GISS-ER','INM-CM3.0','IPSL-CM4','MRI-CGCM2.3.2',$
@@ -148,76 +155,86 @@
window,4,xs=800,ys=800
loadct,5
!p.multi=[0,2,7,0,0]
+
+ print, "Display 4 - Timeseries for AMIP data -2"
+
dummy=label_date(date_format=('%y'))
- plot,julmo,replicate(0,288),color=!myct.black,yr=[-8,8],$
+ plot,julmo,replicate(0,288),color=myctblack,yr=[-8,8],$
ytitle='R (W/m2)',$ ;textoidl('R (Wm^{-2})'),$
xticklen=0.1,xtickunits=['time'],xtickformat='label_date',$
xstyle=1,xminor=5,charsize=2,xmargin=[6,2],$
ystyle=1,yminor=2,yticks=2,ymargin=[1.9,0.3]
- legend,label=['ERBE&CERES'],$
- boxcolor=-5,charsize=1.3,spacing=1.4,halign=0.03,valign=0.91
+ legend,['ERBE&CERES'],$
+ outline_color=-5,charsize=1.3,spacing=1.4 ;,halign=0.03,valign=0.91
+
+print, "Display 4 - Timeseries for AMIP data -3"
oplot,jul,olr,color=112 ;red
oplot,jul,osr,color=47 ;blue
- legend,label=['LW','SW'],line=[0,0],lcolor=[112,47],$
- charsize=1,boxcolor=255,spacing=1
+ legend,['LW','SW'],line=[0,0],colors=[112,47],$
+ charsize=1,outline_color=255,spacing=1
+ print, "Display 4 - Timeseries for AMIP data -4"
for m=0,10 do begin
- plot,julmo,replicate(0,288),color=!myct.black,yr=[-8,8],$
+ print, "Display 4 - Timeseries for AMIP data -5"
+ plot,julmo,replicate(0,288),color=myctblack,yr=[-8,8],$
ytitle='R (W/m2)',$ ;textoidl('R (Wm^{-2})'),$
xticklen=0.1,xtickunits=['time'],xtickformat='label_date',$
xstyle=1,xminor=5,charsize=2,xmargin=[6,2],$
ystyle=1,yminor=2,yticks=2,ymargin=[1.9,0.3]
- legend,label=[title(m)],$
- boxcolor=-5,charsize=1.3,spacing=1.4,halign=0.03,valign=0.91
+ legend,[title(m)],$
+ outline_color=-5,charsize=1.3,spacing=1.4 ;,halign=0.03,valign=0.91
oplot,julmo,ts_smooth(olr_amip(*,m),3),color=112 ;red
oplot,julmo,ts_smooth(osr_amip(*,m),3),color=47 ;blue
endfor
; ---------------
-; Display 5 - Feedback analysis for AMIP 11 models
+print, "Display 5 - Feedback analysis for AMIP 11 models"
; ---------------
window,5,xs=800,ys=600
!p.multi=[0,4,3,0,0]
!x.margin=[7,2]
!y.margin=[3,2]
-; ERBE&CERES results
+print, "ERBE&CERES results"
+
pos=where(CorrDim eq max(CorrDim(5:10)))
plot,LagDim,SlopDim,xr=[-5,5],xs=1,yr=[-4,6],title ='ERBE&CERES',$
xtitle='Lagged month',xminor=2,charsize=2, $
ytitle='R/T (W/m2/K)';textoidl('R/T (Wm^{-2}K^{-1})')
oplot,LagDim,SlopDim2,linestyle=1
- legend,label=['LC','Simple'],line=[0,1],boxcolor=-5,charsize=1.3,spacing=1.4
- legend,label=['F_LC ='+string(SlopDim(pos(0)),format='(f4.1)'),$
+ legend,['LC','Simple'],line=[0,1],outline_color=-5,charsize=1.3,spacing=1.4
+ legend,['F_LC ='+string(SlopDim(pos(0)),format='(f4.1)'),$
'F_sim='+string(SlopDim2(5),format='(f4.1)')],$
- boxcolor=-5,charsize=1.3,spacing=1.4,halign=0.05,valign=0.95
-; 11 AMIP model results
+ outline_color=-5,charsize=1.3,spacing=1.4
+ ;,halign=0.05,valign=0.95
+
+print, '11 AMIP model results'
for m=0,10 do begin ; 11 models
olrmo=olr_amip(*,m) ; AMIP monthly LW flux from 1985
osrmo=osr_amip(*,m) ; AMIP monthly SW flux from 1985
-; Smoothing and scaling
+print, "Smoothing and scaling"
sstmo=ts_smooth(sst*0.78,3)
olrmo=ts_smooth(olrmo,3)
osrmo=ts_smooth(osrmo,3)
-; Call LC11.pro procedure
+print, "Call LC11.pro procedure"
; put 'olrmo' for LW, put 'osrmo' for SW, put 'olrmo+osrmo' for Net
fluxmo=olrmo
lc11,julmo,sstmo,fluxmo,int1,int2,dT,dTmax,dR,dRmax,$
LagDim,SlopDim,ConstDim,CorrDim,SlopDim2,ConstDim2,CorrDim2
;
- ; Lead-lag slopes
+print, "Lead-lag slopes"
pos=where(CorrDim eq max(CorrDim(5:10)))
plot,LagDim,SlopDim,xr=[-5,5],xs=1,yr=[-4,6],title =title(m),$
xtitle='Lagged month',xminor=2,charsize=2, $
ytitle='R/T (W/m2/K)';textoidl('R/T (Wm^{-2}K^{-1})')
oplot,LagDim,SlopDim2,linestyle=1
- ;legend,label=['LC','Simple'],line=[0,1],boxcolor=-5,charsize=1.3,spacing=1.4
- legend,label=['F_LC ='+string(SlopDim(pos(0)),format='(f4.1)'),$
+ ;legend,['LC','Simple'],line=[0,1],outline_color=-5,charsize=1.3,spacing=1.4
+ legend,['F_LC ='+string(SlopDim(pos(0)),format='(f4.1)'),$
'F_sim='+string(SlopDim2(5),format='(f4.1)')],$
- boxcolor=-5,charsize=1.3,spacing=1.4,halign=0.05,valign=0.95
+ outline_color=-5,charsize=1.3,spacing=1.4;,halign=0.05,valign=0.95
endfor
- screen2jpg,'Fig.jpg'
+ screen2jpg,'Fig-gdl.jpg'
end
diff -ur a/lindzen_choi_2011/scripts/obs.pro b/lindzen_choi_2011/scripts/obs.pro
--- a/lindzen_choi_2011/scripts/obs.pro 2011-09-25 18:40:42.000000000 +0200
+++ b/lindzen_choi_2011/scripts/obs.pro 2011-09-27 20:05:35.000000000 +0200
@@ -5,7 +5,7 @@
pro obs,jul,sst,sstmo,olr,osr
-path1='./obs/'
+path1='../obs/'
file1='oisst_36_20n20s_1985_1999.txt'
file2='amipsst_mon_20n20s_1985_2008.txt'
file3='erbe_s10n_36_20n20s_1985_1999.txt'
The last missing function is regress(). The L&C-2011 script seems to use an IDL-7.x version of regress(), which isn’t easily availabe on the net, but is contained as IDL source code within IDL installations. All others missing functions are available:
label_date.pro
legend.pro
textoidl.pro
ts_fcast.pro
julday.pro
ts_coef.pro
ts_smooth.pro
Since these files are mostly from older IDL versions, some minor editing of L&C scripts is needed (such as adapting named parameter to the older definitions or replacing color=”!myct.black” by color=0). If anyone understands the regress() call wrt the parameters “Correlation” and “Const” so far that the call could be adapted to use the older version available on the net, then a GDL version of the L&C scripts could become fully functional. It may make it easier to develop and test the R-version.
Actual 7.x – version of regress()
Result = REGRESS( X, Y, [, CHISQ=variable] [, CONST=variable] [, CORRELATION=variable] [, /DOUBLE] [, FTEST=variable] [, MCORRELATION=variable] [, MEASURE_ERRORS=vector] [, SIGMA=variable] [, STATUS=variable] [, YFIT=variable] )
http://star.pst.qub.ac.uk/idl/REGRESS.html
Obsoleted 5.4 Version of regress()
REGRESS(X, Y, Weights, Yfit, Const, Sigma, Ftest, R, Rmul, Chisq)
http://star.pst.qub.ac.uk/idl/Obs_in_54.html
Here is the IDL 8 version, which should be compatible with IDL 7x
Nice. I haven’t had to learn a new language this week , I was starting to get bored 😉
Good work on getting all this. You certainly seem to have influence.
And thanks of course to L&C for the innovative work and for making this available.
Two questions:
1.)How readily does R work with complex variables?
Because the EE’s need to keep track of magnitude and phase of our non-constant signals, we tend to default to phasor notation to keep both in a single (complex) variable.
2.) Is it correct to say that ARMA processes have delayed attenuated echoes of prior time series values? Is this what they refer to with the keyword ‘autocorrelated’ ?
I’m not seeing how that would be a great tool for diagnosing behaviour where the ‘feedback’ changes a lot depending on the frequency of the signal.
I am not a very mathy person, and so my worldview is not overly cluttered with equations, and flawless jargon to go with them.
Hopefully, something herein is sufficiently WRONGO to elicit a didactic response.
TIA
RR
You might find this summary of value.
geek
accusation or confession?
HMMM could that be said for your last name too 🙂
both
Is the paper/article meant to be there as well? LC11.pdf looks to be the readme converted to a pdf.
Translation to Mathematica is also welcome. 😉
If someone wants a zipped version of the folder, including the LC11 paper itself, you get the link (and a few comments) at:
http://motls.blogspot.com/2011/09/lindzen-choi-2011-paper-data-scripts.html
My small contribution:
Note that array subscripts in IDL start at 0, while in R is at 1. I’ve just added 1 to all array subscripts. This may be confussing for large pieces of code, if you have a better suggestion, please let me know. My commenst in the code start with 3 dashes ###
I have used “=” instead of “<-" it is the same and you have to type less.
amip.pro
Nice.. run it through formatR
Oops,
here he starts with the errors. in last for loop j is used only to get month number (j+1), not for subscripting so I should have left it as it was: (j in 0:11)
Note, if you hoover over the code, in firefox at least, you may copy it or see the code without line numbers
corrected:
The above code porting passes the parser but probably will not work as expected:
R> help (“>
x1 <- 0.5 – 0.3
x2 >
> x1 x2 x1 == x2
[1] FALSE
> x1
[1] 0.2
> x2
[1] 0.2
> x1 == x2
[1] FALSE
Yeah, figure that one out. 0.2 is not the same as 0.2 😕
Probably better to use
if (p<2) {
…
}else{
…
}
Sorry wordpress mangles that a bit.
Yeah, figure that one out. 0.2 is not the same as 0.2 😕
Probably better to use
OK workpress still mangles it when it’s in a protected. You’ve been warned , check out the help on relational operators for the full story.
It works fine here:
“R version 2.12.2 (2011-02-25)”
If wordpress doesn’t mangle it , this should allow grabbing the whole lot in one command.
Careful with the options or else you’ll get loads of other stuff too.
You got me a bit lost. I think it works fine for that case and for non-integers as well:
(use
for the code)
I mean [ sourcecode ] [ / sourcecode ] without spaces
sourcecode between squared brackets to open, same with slash before “s” to close
see the example and warning shown in the help:
I don’t understand the subtle difference in what this assignment is actually assigning that makes the result false.
Neither do I think that having the code path change as a result of the value of data makes for writing reliable code.
It seems to be one of those R specialities like the else clause parsing.
These are things that are likely to cause serious (and subtle) problems for those who do not some considerable air miles fly with R.
Apparently this is a “feature”.
It seems behaviour is stable for integers, so the “if p==1” etc should be OK. Apologies for the false alert.
For floating point equality, use all.equal(), as explained here.
Integer comparison does not present the same problem.
That’s quite true. Thanks for the advice.
It should be changed to
etc.
That criticism should be directed to Lindzen and Choi, I tried to keep the translation as close as possible to the original code.
I wrote too much in a hurry
It seems that the correct form should be
a bit winded R in this case I believe.
Solar ” Neither do I think that having the code path change as a result of the value of data makes for writing reliable code.”
jgc
“That criticism should be directed to Lindzen and Choi, I tried to keep the translation as close as possible to the original code.”
No that is an issue for R developers. It is a serious defect in a language if it shows that kind of instability wrt the values of the data. The pedantic run-around that is necessary to get stable results seems like a sticking plaster workaround.
L&C did not use R so they are not concerned.
See my earlier post. I think your original code was safe since p is only 1 or 2 and I think the problem only exists for (some) real values.
Having drawn attention to this potentially serious problem in porting this to R , I suggest we don’t get bogged down in discussing this further.
My bottom line is that your code was correct.
IDL actually should have the same issues when run on the same hardware. Any language will because the limitation is in the hardware. Not seeing the problems doesn’t necessarily mean they aren’t there.
NASA has a page devoted to IDL and floating point issues:
http://idlastro.gsfc.nasa.gov/idl_html_help/Accuracy_and_Floating_Point_Operations.html
It includes code for checking accuracy and errors.
Replacing TS_SMOOTH or examining what is does would seem critical here since it used some autoregression magic to end-fill in data window to avoid shortening the periods.
I’m instantly wary of this sort of thing , especially here where the data periods are relatively short due to the selection process.
How much of the “data” being analysed is end-fill projections and what do those projections look like in relation to the true data ?
The full description of the routine is here:
http://idlastro.gsfc.nasa.gov/ftp/ittvislib/ts_smooth.pro
; This function computes central, backward or forward moving-averages
; of an n-element time-series (X). Autoregressive forecasting and
; backcasting is used to extrapolate the time-series and compute a
; moving-average for each point of the time-series. The result is an
; n-element vector whose type is identical to X.
Thanks , I’ve seen what it does.
L&C are only smoothing over three points so there is one made up point at each end. Trouble is some of there data segments are only a few months long.
My approach with this is why use artificial data? If it does not make any significant difference, why use it? If it does make a difference DONT use it.
At the very least it will affect the R2 of the result in an artificial way.
Unfortunately I don’t have access to IDL so until it is simplified for GDL or ported to R I can’t investigate further.
Please see my reply above (https://climateaudit.org/2011/09/25/lindzen-choi-2011/#comment-305404)
Best
Thanks, sourcing those routines has been huge help.
However, I’m having trouble running the patch.
when it was in the same dir , patch -p3 could not find the files (typical patch run around). I took the diff file up one level and used -p2 and it now finds them. However, the patch fails straight away.
patching file scripts/amip.pro
patch: **** malformed patch at line 6: ObsName=’ERBE&CERES’
My guess is the diff options used to create the file.
I started to do it by hand but there’s quite a lot in there and you’ve already done the hard work so it seems silly to repeat it.
Can you see what is needed?
TIA.
OK, this seems to be spaces getting stripped out by WP at the beginning of lines.
I tried correcting by hand but I’m still getting errors. This is not a reliable way of doing it.
could you repost and wrap it with [ sourcecode ] [ / sourcecode ] tags (without spaces).
Thanks.
Ok, first testing sourcode tags
a b c d
e f g h
i j k l m >
Sourcecode tags don’t work too; leading spaces still appear to be stripped.
Try this link instead: h ttp://www.pastebin.ch/6669/txt
Hmm, Pastebin.ch adds an additional carriage return to each line, even in the plain text version. If this should become a problem with patch, check “h ttps://gist.github.com/1247625”, and download it either as tar.gz archive or use git clone git://gist.github.com/1247625.git . Besides, it would be very nice to have an upload directory at climateaudit.info, and even nicier if Mr. Pete could set up a git repository for such collaborative efforts.
Well I’m getting somewhere now. I still had a dozen or so hunks that failed in main.pro that I did by hand from the rejects. That looks good.
Some remaining issues in textoidl.
thanks for all your efforts.
BTW the tags worked on the last post there and no space gets lost.
😉
P. Solar,
seems like you have stale version of textoidl which uses textable. Maybe wget had produced something along “textoidl.pro.1” because another textoidl.pro version already existed on your disk? Anyway, simply remove all calls to textoidl and replace the call by plain strings (like ‘W/m^2’), since textoidl doesn’t work as expected anyway.
The rejects appear a bit strange. I tested the patch locally before uploading it. Did you use the patch from gist.github.com?
Re sourcecode tags, probably because I used angle brackets instead of squared brackets.
Btw, I’ve much beter version ready, which does also account for a severe GDL bug ( “PLOTS destroys plot layout with !P.MULTI – ID: 3323178”), fixes the color scheme and also prints the currently skipped plots. I’ll post this as soon as I’ve figured out why XTICKFORMAT isn’t working as expected.
Fine, I got clean copy of everything including github patch and it applied perfectly. Many thanks.
If you suss the other probs please post when you update.
P.Solar, I’ve posted a turnkey shell script, yet it landed in the moderation queue due to the URLs contained therein.
Well Steve, I don’t have any idea what you life victory conditions are, but you have started to people to deposit their data and mathematics on public databases.
That is one hell of a victory.
A day later, my shell script still does not show up here. Until it does, just use the copy at h ttps://gist.github.com/1252987 instead
I’ve ported the obs() procedure to R. obs() now returns a list of two time series, because sst,olr,osr on the one hand and sstmo on the other are of different length. When comparing these time serires to the LC11 GDL output, only small difference show up. The function can be found at h ttps://github.com/hm2/lc11/R. You may also use git clone git://github.com/hm2/lc11 to get the complete environment in a few seconds; the GDL port is also contained therein.
I’ve ported lc11.pro to R too. Rests to find an R equivalent for TS_SMOOTH(), or port ts_smooth.pro. See h ttps://github.com/hm2/lc11
Dear Steve,
Did you ever get around to examining the LC11 method? If it helps you may also be interested in Frankignoul, C. 1999: A Cautionary Note on the Use of Statistical Atmospheric Models in the Middle Latitudes: Comments on “Decadal Variability in the North Pacific as Simulated by a Hybrid Coupled Model”, J. Clim., 12:1871-1872. Also, papers that cite Frankignoul 1999.
One Trackback
[…] Steve just announced via Climate Audit that he was given all the files by Richard and his colleague and you may download the data files […]