Hi All
I am creating EOF's using GPCP data and I have managed to create the same looking patterns for the eofs as the knmi explorer but the scale I have has very small values which do not look right as they should be around -1 and 1 and are currently -0.1 and 0.1. If someone could have a look through my code or just tell me what could potentially be creating such a small scale it would be very much appreciated :)
Many thanks!
Kind Regards
Melissa
Code:
; ==============================================================
; eof_1.ncl
;
; Concepts illustrated:
; - Calculating EOFs
; - Using coordinate subscripting to read a specified geographical region
; - Rearranging longitude data to span -180 to 180
; - Calculating symmetric contour intervals
; - Drawing filled bars above and below a given reference line
; - Drawing subtitles at the top of a plot
;
; ==============================================================
; Calculate EOFs of the Sea Level Pressure over the North Atlantic.
; ==============================================================
; The slp.mon.mean file can be downloaded from:
; http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.surface.html
; ==============================================================
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
begin
; ==============================================================
; User defined parameters that specify region of globe and
; ==============================================================
latS = -40.
latN = 0.
lonL = 10.
lonR = 60.
yrStrt = 1979
yrLast = 2012
season = "DJF" ; choose Dec-Jan-Feb seasonal mean
neof = 3 ; number of EOFs
optEOF = True
optEOF@jopt = 0 ; This is the default; most commonly used; no need to specify.
;;optEOF@jopt = 1 ; **only** if the correlation EOF is desired
optETS = False
; ==============================================================
; Open the file: Read only the user specified period
; ==============================================================
f = addfile ("/mnt/nfs2/geog/ml382/melphd/eof_sicz/gpcp_22.nc", "r")
;f = addfile ("/mnt/nfs2/geog/ml382/melphd/eof_sicz/cmap.prec.mon.mean.nc", "r")
lat = f->lat
TIME = f->time
YYYY = cd_calendar(TIME,-1)/100 ; entire file
iYYYY = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)
pr = f->prcp(iYYYY,0,:,:)
printVarSummary(pr) ; variable overview
;pr1 = flt2dble(pr)
; ==============================================================
; ;Remove annual cycle
; ==============================================================
PR1 = rmMonAnnCycTLL(pr)
; ==============================================================
; compute desired global seasonal mean: month_to_season (contributed.ncl)
; ==============================================================
PR = month_to_season (PR1, season)
nyrs = dimsizes(PR&time)
printVarSummary(PR)
; =================================================================
; normalize data at each gridpoint by local standard deviation at each grid pt
; =================================================================
PR = dim_standardize_n(PR,1,0)
PR = SqrtCosWeight(PR)
; =================================================================
; create weights: sqrt(cos(lat)) [or sqrt(gw) ]
; =================================================================
rad = 4.*atan(1.)/180.
clat = f->lat
clat = sqrt( cos(rad*clat) ) ; gw for gaussian grid
; =================================================================
; weight all observations
; =================================================================
wPR = PR ; copy meta data
wPR = PR*conform(PR, clat, 1)
wPR@long_name = "Wgt: "+wPR@long_name
; =================================================================
; Reorder (lat,lon,time) the *weighted* input data
; Access the area of interest via coordinate subscripting
; =================================================================
x = PR({lat|latS:latN},{lon|lonL:lonR},time|:)
eof = eofunc_Wrap(x, neof, optEOF)
eof_ts = eofunc_ts_Wrap (x, eof, optETS)
printVarSummary( eof ) ; examine EOF variables
printVarSummary( eof_ts )
;Correlations of EOF1 and nino3.4
y= (/0.313333, 0.433333, -0.37, 0.923333, 1.406667, -0.94333, -0.57667, -0.14333, 1.223333, -0.33, -1.20667, 0.183333, 0.816667, 1.303333, 0.186667, 0.4, 0.236667, -0.75, 0.573333, 0.963333, -1.55, -1.47667, -0.61333, 0.56, 0.76, 0.35, 0.03, -0.12667, -0.27333, -1.53667, 0, 0.366667, -1.31667, -0.62667/)
ccr = escorc(eof_ts, y)
print(ccr)
; =================================================================
; 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
; =================================================================
; 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
;============================================================
; PLOTS
;============================================================
wks = gsn_open_wks("X11","eof_gpcp")
gsn_define_colormap(wks,"gsltod") ; choose colormap
plot = new(neof,graphic) ; create graphic array
; only needed if paneling
; EOF patterns
res = True
res@gsnDraw = False ; don't draw yet
res@gsnFrame = False ; don't advance frame yet
;---This resource not needed in V6.1.0
res@gsnSpreadColors = True ; spread out color table
res@gsnAddCyclic = False ; plotted dataa are not cyclic
res@mpFillOn = False ; turn off map fill
res@mpMinLatF = latS ; zoom in on map
res@mpMaxLatF = latN
res@mpMinLonF = lonL
res@mpMaxLonF = lonR
res@cnFillOn = True ; turn on color fill
res@cnLinesOn = False ; True is default
;res@cnLineLabelsOn = False ; True is default
res@lbLabelBarOn = False ; turn off individual lb's
res@tmXBMode = "Explicit" ; Define own tick mark labels.
res@tmXBValues = (/ 10.,20.,30.,40.,50.,59 /)
res@tmXBLabels = (/ "10E","20E","30E","40E","50E","60E" /)
res@tmYLMode = "Explicit" ; Define own tick mark labels.
res@tmYLValues = (/ -39.,-30.,-20.,-10.,-1 /)
res@tmYLLabels = (/"40S","30S","20S","10S","0" /)
res@tmXBLabelFontHeightF = 0.03 ; resize tick labels
res@tmYLLabelFontHeightF = 0.03
;res@pmLabelBarOrthogonalPosF = .10 ; move whole thing down
; set symmetric plot min/max
symMinMaxPlt(eof, 16, False, res) ; contributed.ncl
; panel plot only resources
resP = True ; modify the panel plot
resP@gsnMaximize = True ; large format
resP@gsnPanelLabelBar = True ; add common colorbar
resP@lbLabelAutoStride = True ; auto stride on labels
resP@lbLabelFontHeightF = 0.02
resP@tmXBLabelFontHeightF = 0.03 ; resize tick labels
resP@tmYLLabelFontHeightF = 0.03
;resP@pmLabelBarOrthogonalPosF = .10 ; move whole thing down
res@cnLevelSelectionMode = "ManualLevels" ; manual levels
res@cnMinLevelValF = -0.15 ; min level
res@cnMaxLevelValF = 0.15 ; max level
res@cnLevelSpacingF = 0.025 ; contour spacing
yStrt = yyyymm(0)/100
yLast = yyyymm(nyrs-1)/100
;resP@txString = "pr: "+season+": "+yStrt+"-"+yLast
; reverse the first two colors
setvalues wks
"wkColorMap" : "gsltod"
"wkForegroundColor" : (/0.,0.,0./)
"wkBackgroundColor" : (/1.,1.,1./)
end setvalues
res@lbLabelFontHeightF = 0.01
;*******************************************
; first plot
;*******************************************
do n=0,neof-1
res@gsnLeftString = "EOF "+(n+1)
res@gsnRightString = sprintf("%5.1f", eof@pcvar(n)) +"%"
plot(n)=gsn_csm_contour_map_ce(wks,eof(n,:,:),res)
end do
gsn_panel(wks,plot,(/neof,1/),resP) ; now draw as one plot
;*******************************************
; second plot
;*******************************************
; EOF time series [bar form]
rts = True
rts@gsnDraw = False ; don't draw yet
rts@gsnFrame = False ; don't advance frame yet
rts@gsnScale = True ; force text scaling
; these four rtsources allow the user to stretch the plot size, and
; decide exactly where on the page to draw it.
rts@vpHeightF = 0.40 ; Changes the aspect ratio
rts@vpWidthF = 0.85
rts@vpXF = 0.10 ; change start locations
rts@vpYF = 0.75 ; the plot
rts@tiYAxisString = "mm/month" ; y-axis label
rts@gsnYRefLine = 0. ; reference line
rts@gsnXYBarChart = True ; create bar chart
rts@gsnAboveYRefLineColor = "red" ; above ref line fill red
rts@gsnBelowYRefLineColor = "blue" ; below ref line fill blue
; panel plot only resources
rtsP = True ; modify the panel plot
rtsP@gsnMaximize = True ; large format
;rtsP@txString = "pr: "+season+": "+yStrt+"-"+yLast
year = yyyymm/100
; create individual plots
do n=0,neof-1
rts@gsnLeftString = "EOF "+(n+1)
rts@gsnRightString = sprintf("%5.1f", eof@pcvar(n)) +"%"
plot(n) = gsn_csm_xy (wks,year,eof_ts(n,:),rts)
end do
gsn_panel(wks,plot,(/neof,1/),rtsP) ; now draw as one plot
end
Output:
Variable: pr
Type: float
Total Size: 16920576 bytes
4230144 values
Number of Dimensions: 3
Dimensions and sizes: [time | 408] x [lat | 72] x [lon | 144]
Coordinates:
time: [ 0..407]
lat: [-88.75..88.75]
lon: [1.25..358.75]
Number Of Attributes: 4
level : 0
long_name : precipitation
units : mm/day
_FillValue : 3e+33
Variable: PR
Type: float
Total Size: 1410048 bytes
352512 values
Number of Dimensions: 3
Dimensions and sizes: [time | 34] x [lat | 72] x [lon | 144]
Coordinates:
time: [ 0..396]
lat: [-88.75..88.75]
lon: [1.25..358.75]
Number Of Attributes: 6
anomaly_op_ncl : Annual Cycle Removed: rmMonAnnCycTLL: contributed.ncl
_FillValue : 3e+33
units : mm/day
long_name : DJF: precipitation
level : 0
NMO : 0
Variable: eof
Type: float
Total Size: 3840 bytes
960 values
Number of Dimensions: 3
Dimensions and sizes: [evn | 3] x [lat | 16] x [lon | 20]
Coordinates:
evn: [1..3]
lat: [-38.75..-1.25]
lon: [11.25..58.75]
Number Of Attributes: 7
eval_transpose : ( 5.749153, 2.830663, 2.75108 )
eval : ( 51.74238, 25.47597, 24.75972 )
pcvar : ( 18.30632, 9.013331, 8.759924 )
matrix : covariance
method : transpose
_FillValue : 3e+33
long_name : EOF: DJF: precipitation (sqrt cosine weighted)
Variable: eof_ts
Type: float
Total Size: 408 bytes
102 values
Number of Dimensions: 2
Dimensions and sizes: [evn | 3] x [time | 34]
Coordinates:
evn: [1..3]
time: [ 0..396]
Number Of Attributes: 4
ts_mean : ( -3.974707e-09, 2.176523e-09, 7.091753e-09 )
matrix : covariance
_FillValue : 3e+33
long_name : EOF: Amplitude: DJF: precipitation (sqrt cosine weighted)
Variable: ccr
Type: float
Total Size: 12 bytes
3 values
Number of Dimensions: 1
Dimensions and sizes: [3]
Coordinates:
Number Of Attributes: 1
_FillValue : 3e+33
(0) -0.4300042
(1) 0.5244811
(2) 0.2206682
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Dec 12 08:27:37 2013
This archive was generated by hypermail 2.1.8 : Fri Dec 13 2013 - 11:39:30 MST