[ncl-talk] model ensemble EOF analysis
Alessandra Giannini
alesall at iri.columbia.edu
Fri May 22 10:22:48 MDT 2020
hi again, Sri,
are you saying that you run into a problem when you compute EOFs in the line that says:
eof = eofunc_Wrap(xw, neof, optEOF)
?
Reading the function description <http://www.ncl.ucar.edu/Document/Functions/Contributed/eofunc_Wrap.shtml>
“eofunc_Wrap” expects the “time” dimension to be the rightmost. Which does not seem to be the case for you, because
work := reshape(wSLP ,(/nyrs*nens,nlat,mlon/))
should give you a “work” variable where time is the leftmost dimension.
If that is the issue, you could either reorder the dimensions before feeding them into “eofunc_Wrap”, to ensure that time is rightmost, or use “eofunc_n_Wrap”, which lets you specify where the “time” dimension is.
It is good practice to actually look into the structure of the data after these manipulations, to ensure that values are where you expect them to be, and that dimensions are read in the expected order…
warm regards, alessandra
> On May 22, 2020, at 10:32 AM, Sri nandini <bax8609 at uni-hamburg.de> wrote:
>
> thank you.
>
> I have tried both methods and they would archive different EOF patterns.
>
> The results from the first method (1) is not desirable as it doesnt include the signals from all the ensemble members.
>
> The second method (2) is what i tried, but got errors in my code when combining the time + ensemble dimension-see below. Would anyone point me in the right direction?
>
>
> ;==============================================================
> ; Open the file: Read only the user specified period
> ; =============================================================
> latS = 20.
> latN = 70.
> lonL = -80.
> lonR = 15.
>
> ;yrStrt = 2081
> ;yrLast = 2099
>
> season = "DJF" ; choose Dec-Jan-Feb seasonal mean
>
> neof = 3 ; number of EOFs
> optEOF = True
> optEOF at jopt = 0 ; This is the default; most commonly used; no need to specify.
> optETS = False
>
> ; =============================================================
> f = addfile ("anom_rcp45_zo_slr.nc", "r")
>
> slp=f->rcp45_anom ;[time | 240] x [ens | 100] x [lat | 45] x [lon | 90]
> printVarSummary(slp)
> printMinMax(slp,0)
>
> ; ==============================================================;
> dimx = dimsizes(x)
> ntim = dimx(0) ; 240
> nens = dimx(1) ; 100
> nlat = dimx(2) ; 45
> mlon = dimx(3) ; 90
>
> nmos = 12
> nyrs = ntim/nmos ; 20
>
> ; ==============================================================
> ; compute desired global seasonal mean: month_to_season (contributed.ncl)
> ; The first average (DJF=JF) and the last average (NDJ=ND) are actually two-month averages.
> ; So make climatology and extract the months needed????
> ; ==============================================================
> SLP = month_to_season (slp, season)
> nyrs = dimsizes(SLP&time)
> printVarSummary(SLP)
>
> ;==========================================================================
> ;Estimates and removes the least squares linear trend of the dimension /time/ from all grid points.The mean is also removed. Missing values are not allowed, use dtrend_msg_n if missing values exist.
> ;==========================================================================
>
> SLP = dtrend_n (SLP, False,0)
> printVarSummary(SLP)
>
> ; =================================================================
> ; create weights: sqrt(cos(lat)) [or sqrt(gw) ]
> ; =================================================================
> rad = 4.*atan(1.)/180.
> lat = f->lat
> 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, 1)
> printVarSummary(wSLP)
>
>
> ; ===============================================================
>
> ; compute EOFs on the “total variance”, in which case you would want your two dimensions to be space [lat x lon] and [time x ens member]. If so, then the "time series” would be the concatenation of the time series of each ensemble member. In order to get that, you need to combine time and ens member dimensions into one.
>
> ; ===============================================================
>
> ; work = new((/2000/), "integer")
> ; work!0 = “time”
> ; work&time = time
>
> x = new((/2000/),"integer")
> work = new((/dimsizes(x)),nlat,mlon/), "double", slp at _FillValue)
> printVarSummary(work)
>
> work := reshape(wSLP ,(/nyrs*nens,nlat,mlon/))
> printVarSummary(work) ; (2000,45,90)
>
> copy_VarCoords(slp(0,0,:,:), work)
> printVarSummary(work)
>
> ; =================================================================
> ; Reorder (lat,lon,time) the *weighted* input data
> ; Access the area of interest via coordinate subscripting
> ; =================================================================
>
> xw = work({lat|latS:latN},{lon|lonL:lonR},time|:)
> x = work(time|:,{lat|latS:latN},{lon|lonL:lonR})
>
> ; xw = wSLP({lat|latS:latN},{lon|lonL:lonR},time|:)
> ; x = wSLP(time|:,{lat|latS:latN},{lon|lonL:lonR})
>
> eof = eofunc_Wrap(xw, neof, optEOF)
> eof_ts = eofunc_ts_Wrap (xw, eof, optETS)
>
> printVarSummary( eof ) ; examine EOF variables
> printVarSummary( eof_ts )
>
> ; =================================================================
> ; Normalize time series: Sum spatial weights over the area of used
> ; =================================================================
> eof_ts = dim_standardize_n( eof_ts, 0, 1)
> printVarSummary(eof_ts)
>
> ;====================================================================
> ;my code: Regress
> ;======================================================================
>
> eof_regres=eof ;create an array with meta data
>
> do ne=0,neof-1
> eof_regres(ne,:,:)=(/ regCoef(eof_ts(ne,:), xw) /)
> end do
> printVarSummary(eof_regres)
>
>
> The problem:
>
> After reshaping the detrended winter sea surface height, i cannot perform EOF as the reshaped variable does not have a time dimension. I have tried to add it without success:
>
> work := reshape(wSLP ,(/nyrs*nens,nlat,mlon/))
> printVarSummary(work) ; (2000,45,90)
>
> Sri
>
>
>
>
>
>
>
>
>
> On 22.05.20 10:43, Alessandra Giannini wrote:
>>
>> hi Sri,
>>
>> in general terms, when you have an ensemble, you can do one of two things to compute EOFs:
>>
>> 1) compute EOFs on the ensemble mean, which is what you do here, since you used “dim_avg_n_Wrap” on your variable. In this case the two dimensions that the EOF routine sees are space [lat x lon] and time
>> 2) compute EOFs on the “total variance”, in which case you would want your two dimensions to be space [lat x lon] and [time x ens member]. If so, then the "time series” would be the concatenation of the time series of each ensemble member. In order to get that, you need to combine time and ens member dimensions into one.
>>
>> You may want to look up papers that discuss characterizations of total, forced and internal variance.
>> Here (1) gives you a characterization of “forced” variance, to the extent that the ensemble mean represents it.
>> (2) gives you a characterization of “total” variance.
>> And if you subtracted the ensemble mean from each ensemble member, and then proceeded as in (2), you would get the internal variance.
>>
>> Here are some examples from a google scholar search:
>>
>> Harzallah, A. and Sadourny, R., 1995. Internal versus SST-forced atmospheric variability as simulated by an atmospheric general circulation model. Journal of climate, 8(3), pp.474-495.
>> https://journals.ametsoc.org/doi/abs/10.1175/1520-0442(1995)008%3C0474:IVSFAV%3E2.0.CO;2
>>
>> Ting, M., Kushnir, Y., Seager, R. and Li, C., 2009. Forced and internal twentieth-century SST trends in the North Atlantic. Journal of Climate, 22(6), pp.1469-1481.
>> https://journals.ametsoc.org/doi/full/10.1175/2008JCLI2561.1
>>
>> Venzke, S., Allen, M.R., Sutton, R.T. and Rowell, D.P., 1999. The atmospheric response over the North Atlantic to decadal changes in sea surface temperature. Journal of Climate, 12(8), pp.2562-2584.
>> https://journals.ametsoc.org/doi/full/10.1175/1520-0442%281999%29012%3C2562%3ATAROTN%3E2.0.CO%3B2
>>
>> Reasons why things don’t match could also include whether you have subtracted the climatology or not
>>
>> Hope this helps!
>> alessandra
>>
>>
>>
>>
>> —
>> Alessandra Giannini
>> IRI for Climate and Society - The Earth Institute at Columbia University
>> P.O. Box 1000, Palisades NY 10964-8000, U.S.A.
>>
>> currently at:
>> LMD - École Normale Supérieure
>> 24, Rue Lhomond 75231 PARIS CEDEX 05, France
>>
>>> On May 22, 2020, at 1:01 AM, Sri nandini via ncl-talk <ncl-talk at ucar.edu> wrote:
>>>
>>> Dear all,
>>>
>>> I tried the below code for performing EOF on model ensemble data, with plot attached, but it doesn't match the EOF found in publications since i believe i averaged the ensemble dimension instead of including the signal from it. Can someone point me in the right direction?
>>>
>>>
>>> ;==============================================================
>>> ; Open the file: Read only the user specified period
>>> ; =============================================================
>>>
>>> latS = 20.
>>> latN = 80.
>>> lonL = -70.
>>> lonR = 40.
>>>
>>> ;yrStrt = 1986
>>> ;yrLast = 2005
>>>
>>> season = "DJF" ; choose Dec-Jan-Feb seasonal mean
>>> ; season = "JJA"
>>>
>>> neof = 3 ; number of EOFs
>>> optEOF = True
>>> optEOF at jopt = 0 ; This is the default; most commonly used; no need to specify.
>>>
>>> ; =============================================================
>>> f = addfile ("anom_hist_zo_slr.nc", "r")
>>> x=f->hist_anom ;[time | 240] x [ens | 100] x [lat | 45] x [lon | 90]
>>> printVarSummary(x)
>>> slp=dim_avg_n_Wrap(x,1)
>>> printVarSummary(slp)
>>> printMinMax(slp,0)
>>>
>>> ; ==============================================================
>>> ; compute desired global seasonal mean: month_to_season (contributed.ncl)
>>> ; The first average (DJF=JF) and the last average (NDJ=ND) are actually two-month averages.
>>> ; So make climatology and extract the months needed.
>>> ; ==============================================================
>>> SLP = month_to_season (slp, season)
>>> ; uClm = clmMonTLLL( u )
>>> nyrs = dimsizes(SLP&time)
>>> printVarSummary(SLP)
>>>
>>> ; =================================================================
>>> ; create weights: sqrt(cos(lat)) [or sqrt(gw) ]
>>> ; =================================================================
>>> rad = 4.*atan(1.)/180.
>>> lat = f->lat
>>> 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, 1)
>>> printVarSummary(wSLP)
>>>
>>> ; =================================================================
>>> ; Reorder (lat,lon,time) the *weighted* input data
>>> ; Access the area of interest via coordinate subscripting
>>> ; =================================================================
>>> xw = wSLP({lat|latS:latN},{lon|lonL:lonR},time|:)
>>> x = wSLP(time|:,{lat|latS:latN},{lon|lonL:lonR})
>>>
>>> xw= dtrend(xw(lat|:,lon|:,time|:),False)
>>> printVarSummary(xw)
>>>
>>>
>>> eof = eofunc_Wrap(xw, neof, optEOF)
>>> eof_ts = eofunc_ts_Wrap (xw, eof, optETS)
>>>
>>> printVarSummary( eof ) ; examine EOF variables
>>> printVarSummary( eof_ts )
>>> ; =================================================================
>>> ; Normalize time series: Sum spatial weights over the area of used
>>> ; =================================================================
>>> eof_ts = dim_standardize_n( eof_ts, 0, 1)
>>> printVarSummary(eof_ts)
>>>
>>> ;====================================================================
>>> ;my code: Regress
>>> ;======================================================================
>>>
>>> eof_regres=eof ;create an array with meta data
>>> ;eof_ts(0,:)=-eof_ts(0,:)
>>> do ne=0,neof-1
>>> eof_regres(ne,:,:)=(/ regCoef(eof_ts(ne,:), xw) /)
>>> end do
>>> printVarSummary(eof_regres)
>>>
>>> ; eof_regres=-eof_regres
>>>
>>> ;=================================================================
>>> yyyymm = cd_calendar(eof_ts&time,-2)/100
>>> ;============================================================
>>>
>>> wks = gsn_open_wks("pdf","hist_SSH_EOF2")
>>>
>>> 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 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 cnLineThicknessF = 2
>>> ;res at cnLineColor = 0
>>>
>>> res at mpMinLatF = latS ; zoom in on map
>>> res at mpMaxLatF = latN
>>> res at mpMinLonF = lonL
>>> res at mpMaxLonF = lonR
>>> ;res at mpFillOn = False ; turn on map fill
>>> res at mpOutlineOn = True ; turn the map outline on
>>>
>>> res at cnLevelSelectionMode = "ManualLevels" ; manual set levels
>>> res at cnMinLevelValF = -0.5
>>> res at cnMaxLevelValF = 0.5
>>> res at cnLevelSpacingF = .1 ; 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
>>> ; now change the size of the label bar labels
>>> resP at lbLabelFontHeightF = 0.017
>>>
>>> ;====================Create (but don't draw) both plots
>>> do n=0, neof -1
>>> res at gsnLeftString = "SSH DJF EOF "+(n+1)
>>> res at gsnRightString = sprintf("%5.1f", eof_regres 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
>>>
>>>
>>>
>>> The EOF1 should resemble EOF3 pattern instead, as per previous literature for sea surface height.
>>>
>>> Can someone help me with this?
>>>
>>> Sri
>>>
>>>
>>>
>>>
>>>
>>> On 21.05.20 17:18, Sri nandini via ncl-talk wrote:
>>>> Thank you.
>>>>
>>>> The idea is that SSH variability is computed over the ensemble dimension, rather than the traditional time-dimension- i.e computing EOF across individual ensemble member at each time step with looping, in which case the below would not make sense .i.e averaging the ensemble dimension?
>>>>
>>>> Perhaps i should perform the EOF on both the time and ensemble dimension as variations across time and ensemble are supposed to be the same after detrending, the more samples the more robust the EOF is?
>>>>
>>>> Best
>>>>
>>>> Sri
>>>>
>>>>
>>>>
>>>> On 21.05.20 05:09, Dennis Shea wrote:
>>>>> At each time stelp,and lat/lon location all 100 ensembel members
>>>>>
>>>>>
>>>>> ssh_ens_mean = dim_avg_n_Wrap(ssh,1) ; average
>>>>> printVarSummary(ssh_ens_mean)
>>>>> printMinMax(ssh_ens_mean,0)
>>>>>
>>>>> --
>>>>> Input 'ssh_ens_mean' as you would any other variable to the eof function
>>>>>
>>>>>
>>>>> On Wed, May 20, 2020 at 7:29 AM Sri nandini via ncl-talk <ncl-talk at ucar.edu> wrote:
>>>>> Hello dear fellow ncl users,
>>>>>
>>>>> I have been analyzing and plotting the standard EOF with success. Now i
>>>>> wish to proceed onto model ensemble EOF but having problems
>>>>> understanding this coding.
>>>>>
>>>>> My original data is in this format: SSH= [time | 240] x [ens | 100] x
>>>>> [lat | 45] x [lon | 90]
>>>>> How can i modify the standard EOF script on the NCL page to perform it
>>>>> on model ensemble of 100 members?e.g through looping it?
>>>>>
>>>>> Would be grateful for an example.
>>>>>
>>>>> best
>>>>>
>>>>> sri
>>>>>
>>>>>
>>>>>
>>>>> _______________________________________________
>>>>> ncl-talk mailing list
>>>>> ncl-talk at ucar.edu
>>>>> List instructions, subscriber options, unsubscribe:
>>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>>
>>>>
>>>> _______________________________________________
>>>> ncl-talk mailing list
>>>>
>>>> ncl-talk at ucar.edu
>>>>
>>>> List instructions, subscriber options, unsubscribe:
>>>>
>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>> <SSH_historical_DJF_EOF.pdf>_______________________________________________
>>> ncl-talk mailing list
>>> ncl-talk at ucar.edu
>>> List instructions, subscriber options, unsubscribe:
>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
More information about the ncl-talk
mailing list