[ncl-talk] sub: set lat issue..
Adv
advita6 at gmail.com
Sat Jun 4 21:06:31 MDT 2016
Hi,
If I set latS=40,latN=49, I didn't get the actual region corresponding to
those latitude. Could someone help me to fix this?
Thank you
yrStrt = 1965
yrLast = 2005
season = "DJF" ; choose Dec-Jan-Feb seasonal mean
;####################################################
f = addfile("air.2m.mon.mean.nc", "r") ; note the "s" of addfile
u = f->air
printVarSummary (u)
u=u-273.15
;anew=u(lat|:,lon|:,time|1380:)
anew=u(lat|:,lon|:,time|204:695)
printVarSummary(anew)
;print(anew)
;return
;############################################################
;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)
origanew1a = month_to_season (temp(time|:,lat|:,lon|:), season)
printVarSummary(origanew1a)
nyrs = dimsizes(anew1a&time)
printVarSummary(anew1a)
printMinMax(anew1a,False)
;Globe
latS = 37.12
latN = 49.70
lonL = 0-116.5
lonR = 0-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))
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)
;========================================================
;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
printVarSummary( eof_ts )
printVarSummary(origanew1a)
;print(eof_ts)
;print(eof(0,:,:))
; =================================================================
; Regress
; =================================================================
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
printVarSummary(eof_regres)
;===================================================================
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)
;**************************************************
; plot parameters
;**************************************************
;************************************************
; plotting parameters
;************************************************
wks = gsn_open_wks("x11","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 = -2.9 ; set min contour level
res at cnMaxLevelValF = 2.0 ; set max contour level
res at cnLevelSpacingF = 0.2 ; 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) ; 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 = "Reanal_Temp(~S~o~N~C): "+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
; txres = True
; txres at txFontHeightF = 0.015
; gsn_text_ndc(wks,"Figure 16: A smaller panel plot",0.5,0.16,txres)
frame(wks)
end
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160604/acbaadc7/attachment.html
More information about the ncl-talk
mailing list