[ncl-talk] EOF PC timeseries

[Sri nandini]

In particular i wanted to understand how to decompose or "reverse the 
reshape" to get back the original dimension data for plotting the EOF 
timeseries.

The original data is [time|240,ens|100,lat,lon]

The reshaped data is [24000,lat,lon]

I wish to get back the original data from my reshaped data. does NCL 
have a function for this?

Sri

On 13.04.21 14:26, Sri nandini via ncl-talk wrote:

> Hello ncl-users,
>
> I have an EOF script, which im applying to my data over the 
> time+ensemble level. The EOF spatial pattern comes out correctly, 
> however, the timeseries does not get plotted. Could someone help me 
> out with my code below?
>
>   begin
>
>   latS   =  20.
>   latN   =  68.
>   lonL   = -78.
>   lonR   =  15.
>
>   neof   = 2       ; number of EOFs
>   optEOF = True
>   optEOF at jopt = 0   ;uses correlation matrix to compute EOFs. The 
> default is to use a covariance matrix (optEOF at jopt = 0).
>   optETS = False
>   yrStrt                                   = 1986
>   yrLast                                   = 2006
>
>   f     = addfile ("data.nc", "r")
>    time   = f->time
>    YYYY   = cd_calendar(time,-1)/100
>    iYYYY  = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)
>    slp=f->hist_trend(iYYYY,:,{latS:latN},{lonL:lonR})
>    printVarSummary(slp)    ;[240,100,13,24]
>
>   slp = dtrend_msg_n(ispan(0,dimsizes(slp&time)-1,1),slp,True,False,0) 
> ; detrend the data
> ; =================================================================
> ; create weights:
> ; =================================================================
>   rad    = 4.*atan(1.)/180.
>   lat = f->lat({latS:latN})
>    if (typeof(lat).eq."double") then
>        clat = sqrt( cos(rad*tofloat(lat)) )
>    else
>        clat = sqrt( cos(rad*lat) )
>    end if
>    copy_VarCoords(lat, clat) ; contributed
>    printVarSummary(clat)
> ; =================================================================
> ; weight all observations
> ; =================================================================
>   wSLP   = slp  ; type float
>   wSLP   = slp*conform(slp, clat, 2)
>   printVarSummary(wSLP)
> ; =================================================================
> ; Reorder (lat,lon,time) the *weighted* input data Access the area of 
> interest via coordinate subscripting
> ; =================================================================
>
>   work := reshape(wSLP ,(/ntim*nens,nlat,mlon/))    ;concatanate the 
> time+ensemeble dimension together to get total variance
>   copy_VarCoords(wSLP(0,0,:,:), work(0,:,:))      ; trick to copy 
> coord information
>   printVarSummary( work )
>
>   eof    = eofunc_n_Wrap(work,neof,optEOF,0 )
>   eof_ts = eofunc_ts_n_Wrap(work, eof, optETS,0)
>   printVarSummary( eof )                         ; should be [3,13,24]
>   printVarSummary( eof_ts )                      ; should be [3,2400]
>
>    eof_ts = dim_standardize_n( eof_ts, 0, 1)
>   ;printVarSummary(eof_ts)
>
> ;====================================================================
> ;Regressnig data onto standardized PC1,2
> ;======================================================================
>
>   eof_regres=eof ;create an array with meta data
>   printVarSummary(eof_regres)
>
>   do ne=0,neof-1
>    eof_regres(ne,:)=(/ regCoef_n(eof_ts(ne,:), work,0,0) /)
>     end do
>   printVarSummary(eof_regres)
>
>   plotfile="EOF_figure"
>   wks_type = "png"
>   wks = gsn_open_wks(wks_type,plotfile)
>   plot                    = new(neof,graphic) ; create graphic array; 
> only needed if paneling
>   res                     = True
>   res at mpProjection        = "LambertConformal"; choose projection
>   res at gsnDraw             = False           ; don't draw
>   res at gsnFrame            = False           ; don't advance frame
>
>   res at cnFillOn            = True            ; turn on color
>   res at cnLinesOn   = False                       ; turn off the contour 
> lines
>   res at cnLineThicknessF = 0.9
>   res at cnLineColor = 0
>   res at mpDataBaseVersion   = "MediumRes"
>   res at cnLineLabelsOn      = False    ; turn off contour line labels
>
>   res at cnFillDrawOrder     = "PreDraw"       ; draw contours before 
> continents
>   res at cnFillPalette = "BlRe"
>   res at mpMinLatF            = latS           ; zoom in on map
>   res at mpMaxLatF            = latN
>   res at mpMinLonF            = lonL
>   res at mpMaxLonF            = lonR
>
>   res at mpOutlineOn = True                        ; turn the map outline on
>
>   res at cnLevelSelectionMode = "ManualLevels"     ; manual set levels
>   res at cnMinLevelValF       = -0.4
>   res at cnMaxLevelValF       = .4
>   res at cnLevelSpacingF      = 0.05                 ; 20 contour levels
> res at lbLabelBarOn         = False              ; turn off individual lb's
>   res at lbBoxEndCapStyle     = "TriangleBothEnds" ; Added in NCL V6.4.0
>   res at lbLabelAutoStride    = True               ; Control labelbar 
> spacing
>   res at gsnMaximize          = True               ; large format in 
> landscape
>   res at gsnAddCyclic         = False              ; plotted dataa are 
> not cyclic
>   res at gsnMaskLambertConformal = True            ; turn on lc masking
>
> ;=============================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 lbLabelFontHeightF = 0.02
>   res at gsnCenterStringFontHeightF = 0.02              ; ditto
> ;====================Create (but don't draw) both plots
>
>      do n=0,neof-1
>      res at gsnStringFontHeightF = 0.02
>      res at gsnLeftString  = "EOF "+(n+1)
>      res at gsnCenterString= "a) Hist SSH"
>      res at gsnRightString = sprintf("%5.1f", eof at pcvar(n)) +"%"
>      plot(n)=gsn_csm_contour_map(wks,eof_regres(n,:,:),res)
>    end do
>   gsn_panel(wks,plot,(/neof,1/),resP)     ; now draw as one plot
>
> ;============================================================
> ; time series (principal component) plot
> ;============================================================
>    eof_ts_ens_avg at long_name = "Standardized"
>
>    yyyymm = cd_calendar(slp&time, -1)         ; create a 
> yyyy.fraction_of_year
>    eof_ts_ens_avg&time = yyyymm_to_yyyyfrac(yyyymm, 0.0) ; better time 
> axis
>    ;eof_ts&time at long_name = "YEAR"
>    ntim   = dimsizes(yyyymm)
>
>    ;printVarSummary(frac)
>    rts           = True
>    rts at gsnDraw   = False       ; don't draw yet
>    rts at gsnFrame  = False       ; don't advance frame yet
>    rts at tmXBFormat= "f"         ; no unnessary decimal pts
>
>    rts at vpHeightF = 0.40        ; Changes the aspect ratio
>    rts at vpWidthF  = 0.85
>    rts at vpXF      = 0.10        ; change start locations
>    rts at vpYF      = 0.75        ; the plo
>    rts at tiYAxisString = "Standardized"          ; y-axis label
>    rts at trXMinF = 1986
>    rts at trXMaxF = 2006
>    rts at trYMinF = -3
>    rts at trYMaxF =  3
>    rts at gsnYRefLine           = 0.              ; reference line
>    rts at gsnXYBarChart         = True            ; create bar chart
>    rts at gsnAboveYRefLineColor = "red"           ; above ref line fill red
>    rts at gsnBelowYRefLineColor = "blue"          ; below ref line fill blue
>
>  ; panel plot only resources
>    rtsP                      = True            ; modify the panel plot
>    rtsP at gsnMaximize          = True            ; large format
>  ; create individual plots
>   do n=0,neof-1
>      rts at gsnLeftString  = "EOF "+(n+1)
>      rts at gsnRightString = sprintf("%5.1f", eof_regres at pcvar(n)) +"%"
>      plot(n) = gsn_csm_xy (wks,work(:,0,0),eof_ts(n,:),rts)
>   end do
>   gsn_panel(wks,plot,(/neof,1/),rtsP)     ; now draw as one plot
> end
>
>
-- 
Dr. Sri, Nandini-Weiss
Research Associate

Universität Hamburg
Center for Earth System Research and Sustainability (CEN)
Cluster of Excellence 'Climate, Climatic Change, and Society' (CLICCS)

Bundesstrasse 53, 20146 Hamburg