Re: Scale for EOF's

Hi Melissa,

the EOF routine in NCL returns normalized EOFs. Did you check if the
results you are comparing
against are normalized or do they have the same units as the original data?

If you look at the EOF routine description here
http://www.ncl.ucar.edu/Document/Functions/Built-in/eofunc.shtml
the last example tells you how to denormalize the EOF patterns (make the
EOFs have the same
units as the original data) by multiplying with the square root of the
eigenvalues. The first
paragraph in the description on the same page also mentions that.

I hope this helps,
Maria

> 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
>

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Dec 12 12:33:59 2013