[ncl-talk] Sub: error
Adv
advita6 at gmail.com
Mon Jun 6 16:54:03 MDT 2016
Hi ,
Could someone help me to fix this bug? I don't whats the fatal error says
about.
Thank you
;***************************************************
begin
;***************************************************
yrStrt = 196501
yrLast = 200512
season = "SON" ; choose Dec-Jan-Feb seasonal mean
;####################################################
............
f = addfile("hgt.mon.mean.nc", "r") ; note the "s" of addfile
u = f->hgt
printVarSummary (u)
YYYYMM = cd_calendar( f->time, -1)
iStrt = ind(YYYYMM.eq.yrStrt)
iLast = ind(YYYYMM.eq.yrLast)
n=u(level|:,lat|:,lon|:,time|iStrt:iLast)
anew=n({500},:,:,:)
printVarSummary (anew)
;############################################################
;anomalies-------------------------
temp= rmMonAnnCycTLL(anew(time|:,lat|:,lon|:))
printVarSummary(temp)
;dtrend the values/------------------
anew1Dtrend = dtrend(temp(lat|:,lon|:,time|:),True)
printVarSummary(anew1Dtrend)
;print(anew1Dtrend at slope)
;print(anew1Dtrend at y_intercept)
copy_VarCoords(anew,anew1Dtrend)
; ==============================================================
; compute desired global seasonal mean: month_to_season (contributed.ncl)
; ==============================================================
anew1a = month_to_season (anew1Dtrend(time|:,lat|:,lon|:), season)
nyrs = dimsizes(anew1a&time)
printVarSummary(anew1a)
printMinMax(anew1a,False)
;%%%%%%%%%%%%%%%EOF%%%%%%%%%%%%%%%%%%%%%%%%
;Dont change at all . Here is a trick..
;Globe
latS = 37
latN = 49.
lonL = -116
lonR = -90
;Center=0
;%%%%%%%%%%%%%%%EOF%%%%%%%%%%%%%%%%%%%%%%%%
neof = 3 ; number of EOFs
optEOF = True
optEOF at jopt = 0 ; This is the default; most commonly used; no need to
specify.
;;optEOF at jopt = 1 ; **only** if the correlation EOF is desired
; ==============================================================
; dataset longitudes span 0=>357.5
; Because EOFs of the North Atlantic Oscillation are desired
; use the "lonFlip" (contributed.ncl) to reorder
; longitudes to span -180 to 177.5: facilitate coordinate subscripting
; ==============================================================
anew1a = lonFlip(anew1a)
printVarSummary(anew1a) ; note the longitude
coord
orig=anew1a({lat|latS:latN},{lon|lonL:lonR},time|:)
printVarSummary(orig)
; =================================================================
; create weights: sqrt(cos(lat)) [or sqrt(gw) ]
; =================================================================
rad = 4.*atan(1.)/180.
clat = anew1a(0,:,0)
printVarSummary(clat)
clat = sqrt( cos(rad*clat) ) ; gw for gaussian grid
printVarSummary(clat)
; =================================================================
; weight all observations
; =================================================================
wSLP = anew1a ; copy meta data
printVarSummary(wSLP)
printMinMax(wSLP,False)
wSLP = anew1a*conform(anew1a, clat, 1)
wSLP at long_name = "Wgt: "+anew1a at long_name
printVarSummary(wSLP)
printMinMax(wSLP,False)
; =================================================================
; Reorder (lat,lon,time) the *weighted* input data
; Access the area of interest via coordinate subscripting
; =================================================================
x = wSLP({lat|latS:latN},{lon|lonL:lonR},time|:)
printVarSummary(x)
print(min(x))
print(max(x))
dimx = dimsizes( x )
mln = dimx(1)
sumWgt = mln*sum( clat({lat|latS:latN}) )
printMinMax(sumWgt,False)
eof = eofunc_Wrap(x, neof, optEOF)
printVarSummary(eof)
print(eof(0,:,:))
;return
;;;;EOF spatial patterns may be tested for orthogonality by using the dot
product:
d01 = sum(eof(0,:,:)*eof(1,:,:))
d12 = sum(eof(1,:,:)*eof(2,:,:))
d02 = sum(eof(0,:,:)*eof(2,:,:))
print("d01="+d01+" d12="+d12+" d02="+d02) ; may be +/- 1e-8
;return
printMinMax(eof(0,:,:),False)
;To retain same unit
do ne=0,neof-1
eof(ne,:,:) = eof(ne,:,:)*sqrt( eof at eval(ne) )
end do
;print(eof(:,:,:))
printMinMax(eof(0,:,:),False)
;return
;========================================================
;Print the Percentage variance;;;;;;
; print("% var="+ eof at pcvar(0))
;========================================================
eof_ts = eofunc_ts_Wrap (x, eof, False)
printVarSummary( eof_ts )
print( eof_ts )
;=======================================================
;To test the EOF time series for orthogonality, compute correlations;;;;;
r01 = escorc(eof_ts(0,:), eof_ts(1,:))
r12 = escorc(eof_ts(1,:), eof_ts(2,:))
r02 = escorc(eof_ts(0,:), eof_ts(2,:))
print("r01="+r01+" r12="+r12+" r02="+r02) ; numbers may be +/- 1e-8
; =================================================================
; Normalize time series: Sum spatial weights over the area of used
; =================================================================
; dimx = dimsizes( x )
; mln = dimx(1)
; sumWgt = mln*sum( clat({lat|latS:latN}) )
;eof_ts = eof_ts/sumWgt
;nlat = dimsizes(eof&lat) ; lat==> whatever the coordinate name is
; mlon= dimsizes(eof&lon)
;print(eof_ts)
eof_ts = dim_standardize_n( eof_ts, 0, 1) ; normalize
eof_regres = eof ; create an array w meta
data
do ne=0,neof-1
eof_regres(ne,:,:) = (/ regCoef(eof_ts(ne,:),orig(lat|:,lon|:,time|:))
/)
end do
printMinMax(eof_regres(0,:,:),False)
printVarSummary( eof_regres )
;print(eof_ts)
;print(eof(0,:,:))
;===================================================================
nlat = dimsizes(eof&lat) ; lat==> whatever the coordinate name is
mlon= dimsizes(eof&lon)
;------------------------------------------------------------
; =================================================================
; Extract the YYYYMM from the time coordinate
; associated with eof_ts [same as SLP&time]
; =================================================================
yyyymm = cd_calendar(eof_ts&time,-2)/100
;;yrfrac = yyyymm_to_yyyyfrac(yyyymm, 0.0); not used here
;*******************************************
; North significance test: any pcvar, eval_transpose, eval can be used
;*******************************************
dimp = dimsizes(x)
ntim = dimp(0) ; max #
eigenvalues possible
prinfo = True
sig = eofunc_north(eof at pcvar, ntim, prinfo)
;************************************************
; plotting parameters
;************************************************
;wks =
gsn_open_wks("x11","SlpReconst_Anom_eof_Reanal_"+season+"_temp1965-2005")
wks =
gsn_open_wks("x11","heightReconst_Anom_eof_Reanal_"+season+"_temp1965-2005")
gsn_define_colormap(wks,"NCV_jaisnd")
plot = new(neof,graphic) ; create graphic array
; only needed if paneling
; EOF patterns
res = True
res at gsnDraw = False ; don't draw yet
res at gsnFrame = False ; don't advance frame yet
;---This resource not needed in V6.1.0
res at gsnSpreadColors = True ; spread out color table
res at cnLevelSelectionMode = "ManualLevels" ; set manual contour levels
res at cnMinLevelValF = -10.0 ; set min contour level
res at cnMaxLevelValF = 10.0 ; set max contour level
res at cnLevelSpacingF = 0.5 ; set contour spacing
res at gsnAddCyclic = False ; plotted dataa are not cyclic
res at mpFillOn = False ; turn off map fill
res at mpMinLatF = eof&lat(0)+3 ; zoom in on map
res at mpMaxLatF = eof&lat(nlat-1)
res at mpMinLonF = eof&lon(0)
res at mpMaxLonF = eof&lon(mlon-1)
res at cnLineLabelsOn = False
res at cnFillOn = True ; turn on color fill
res at cnLinesOn = False ; True is default
res at lbLabelBarOn = False ; turn off individual lb's
; res at mpPerimOn = True ; draw box around map
res at mpGeophysicalLineThicknessF = 3.0
res at mpGeophysicalLineColor = "Black"; (/22/)
res at mpOutlineBoundarySets = "GeophysicalAndUSStates" ; add state boundaries
res at mpNationalLineColor = res at mpGeophysicalLineColor
res at mpUSStateLineThicknessF = 3.0
res at mpUSStateLineColor = res at mpGeophysicalLineColor
; set symmetric plot min/max
; symMinMaxPlt(eof, 16, False, res) ; contributed.ncl
; panel plot only resources
resP = True ; modify the panel plot
resP at gsnMaximize = True ; large format
resP at gsnPanelLabelBar = True ; add common colorbar
resP at lbLabelAutoStride = True ; auto stride on labels
resP at lbLabelFontHeightF = 0.017 ; label bar font
yStrt = yyyymm(0)/100
yLast = yyyymm(nyrs-1)/100
resP at txString = "Regres_Reanal_height(m): "+season+":
"+yStrt+"-"+yLast
;*******************************************
; first plot
;*******************************************
do n=0,neof-1
res at gsnLeftStringFontHeightF = 0.02 ; instead of using
txFontHeightF or gsnStringFontHeightF
res at gsnCenterStringFontHeightF = 0.02 ; to set the
gsnLeft/Center/RightString font heights,
res at gsnRightStringFontHeightF = 0.02 ; individually set each
string's font height.
res at gsnLeftString = "EOF "+(n+1)
res at gsnCenterString= "NORTH="+sig(n)
res at gsnRightString = sprintf("%5.1f", eof at pcvar(n)) +"%"
plot(n)=gsn_csm_contour_map_ce(wks,eof_regres(n,:,:),res)
end do
gsn_panel(wks,plot,(/neof,1/),resP) ; now draw as one plot
frame(wks)
end
ERROR printed:
fatal:Loop variable must be scalar, can't execute loop
fatal:["Execute.c":8575]:Execute: Error occurred at or near line 242 in
file heightReconstruct_Temp_eof_Reanal.ncl
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160606/4fe1c40b/attachment.html
More information about the ncl-talk
mailing list