Re: Wind anomaly, MJO_PC_INDEX and MJO clivar (still different ?)

From: juki juki <juky_emc2_at_nyahnyahspammersnyahnyah>
Date: Tue Jun 12 2012 - 18:06:22 MDT

Mary, thanks for reply.  >We have updated the image on the web, and it should now match what you are getting. Yes, if I use the wind and olr anomalies provided by Dennis, I got the same result as your plot (web). However, if I try to use wind anomaly based on my calculation/code, the result is still different from the example (please see example). I am pretty sure that the discrepancy because of wind anomaly. I am not sure whether  my way to merge many files is correct or not. Probably you could suggest some inprovement or some possible errors. Please see the code below: ---------------------------------------------    ymdStrt = 19790101                         ; start yyyymmdd    ymdLast = 20111231                         ; last      yrStrt  = ymdStrt/10000    yrLast  = ymdLast/10000   nhar    = 4                                ; number of fourier comp          lev  = 200  ; desired level    var     = "uwnd"                            ; name of file                       vName   = "UWND"                            ; name for plots          ;========================   ; get list of all files and open as "one big file"   ;========================                                   all_files = systemfunc ("ls ./Data/uwnd*.nc")      f      = addfiles (all_files, "r")   ; note the "s" of addfile      print(all_files)   ;========================   ; choose how files are combined and read in variable across files   ;========================        ListSetType (f, "cat")             ; concatenate or "merge" (default)         if (NC) then        diro= "./Anom/"        ; output dir        filo= var+".day."+lev+".anomalies."+yrStrt+"-"+yrLast+".nc"  ; output file      end if    if (PLOT) then        wksType = "ps"        wksName = "mjoclivar"                  ; "mjo."+yrStrt+"_"+yrLast    end if ;*********************************************************** ; Read user specified time and create required yyyyddd                     ;***********************************************************                  time    = f[:]->time                          ; all times on file    ymd     = cd_calendar(time, -2)            ; yyyymmdd    iStrt   = ind(ymd.eq.ymdStrt)              ; index start    iLast   = ind(ymd.eq.ymdLast)              ; index last     delete(time)    delete(ymd) ;*********************************************************** ; Read user specified time and create required yyyyddd                     ;***********************************************************    time    = f[:]->time(iStrt:iLast)             ; time:units = "hours since"    TIME    = cd_calendar(time, 0)             ; type float     year    = floattointeger( TIME(:,0) )    month   = floattointeger( TIME(:,1) )    day     = floattointeger( TIME(:,2) )     ddd     = day_of_year(year, month, day)     yyyyddd = year*1000 + ddd                  ; needed for input ;*********************************************************** ; Read data: short2flt                                      ;***********************************************************     x       =  short2flt(f[:]->$var$(iStrt:iLast,{lev},:,:))    ; convert to float     printVarSummary(x)  .....after this part the code is the same as in the web. --------------------- Actually, I can use Dennis data, but it is also better to know/understand the problem of my calculation because sometimes I may need to analyze another data. thanks for help. Juki ________________________________ From: Mary Haley <haley_at_ucar.edu> To: juki juki <juky_emc2_at_yahoo.com> Cc: NCL Talk <ncl-talk_at_ucar.edu> Sent: Saturday, June 9, 2012 12:56 AM Subject: Re: Wind anomaly, MJO_PC_INDEX and MJO clivar (still different ?) Juki, We believe the plot on the web was simply wrong. It was flipped in the Y axis.  We have updated the image on the web, and it should now match what you are getting. --Mary On Jun 7, 2012, at 4:41 PM, juki juki wrote: > Thanks for reply Dennis; > > >If the daily anomalies are computed correctly, everything works as in > >the examples. > > I tried to plot the data that you shared at the ftp. However, the result is also very different from the example. Please see the plot for mjoclivar_15 (enclosed file) and compare with the example (http://www.ncl.ucar.edu/Applications/Images/mjoclivar_15_lg.png) > > You can see clearly the location of phase for each month is totally different, for instance: MJO phase of December (red color). > > You said that the result is only slightly different but I try in my PC it is totally different. Could you share your plot (mjoclivar_15_lg.png)  from  1979-2011 data anomaly ? I don't hope the error coming from my computer (cygwin). > > Thanks again, > > Best regards, > Jouke > > > > > > > From: Dennis Shea <shea_at_ucar.edu> > To: juki juki <juky_emc2_at_yahoo.com> > Cc: NCL Talk <ncl-talk_at_ucar.edu> > Sent: Thursday, June 7, 2012 11:15 PM > Subject: Re: [ncl-talk] Wind anomaly, MJO_PC_INDEX and MJO clivar > > I reran *all* the 'MJO Clivar' diagnostics. Everything works correctly. > > [1] downloaded the NCEP data from the website below. Period 1979-2011 > > [2] computed the anomalies. These are *slightly* different due to >    different climatological base periods (1980-2005 and 1979-2011). > > [3] redid all the computations and plots > > If the daily anomalies are computed correctly, everything works as in > the examples. > > --- > I have placed the anomalies and the scripts used in an ftp area. > > ftp ftp.cgd.ucar.edu > anonymous > email > cd pub/shea/MJO > prompt > mget * > quit > --- > > The files and sizes are: >        2636 Jun  7 07:52 mjoclicar_10.ncl >        4786 Jun  7 07:52 mjoclivar_11.ncl >        7275 Jun  7 07:52 mjoclivar_12.ncl >      11249 Jun  7 07:52 mjoclivar_13.ncl >      12992 Jun  7 07:52 mjoclivar_14.ncl >        6635 Jun  7 07:57 mjoclivar_15.ncl >        9302 Jun  7 08:03 mjoclivar_16.ncl >      10530 Jun  7 07:52 mjoclivar_2.ncl >      10374 Jun  7 07:52 mjoclivar_2.olr.ncl >      10530 Jun  7 07:52 mjoclivar_2.uwnd_200.ncl >      10530 Jun  7 07:52 mjoclivar_2.uwnd_850.ncl >      10530 Jun  7 07:52 mjoclivar_2.vwnd_200.ncl >      10530 Jun  7 07:52 mjoclivar_2.vwnd_850.ncl >        8361 Jun  7 07:52 mjoclivar_3.ncl >        3169 Jun  7 07:52 mjoclivar_4.ncl >        2876 Jun  7 07:52 mjoclivar_5.ncl >        5199 Jun  7 07:52 mjoclivar_6.ncl >        2750 Jun  7 07:52 mjoclivar_7.ncl >        8678 Jun  7 07:52 mjoclivar_8.ncl >        7137 Jun  7 07:52 mjoclivar_9.ncl > >      37400 Jun  7 07:43 MJO_PC_INDEX.nc >  1013804160 Jun  7 07:43 olr.day.anomalies.1979-2011.nc >  615513004 Jun  7 07:43 pregpcp19962008.daily.nc >  350655700 Jun  7 07:44 TRMM.200006-200109.thin.nc >  1013804216 Jun  7 07:44 uwnd.day.anomalies.1979-2011.200.nc >  1013804216 Jun  7 07:44 uwnd.day.anomalies.1979-2011.850.nc >  1013804216 Jun  7 07:45 vwnd.day.anomalies.1979-2011.200.nc >  1013804216 Jun  7 07:45 vwnd.day.anomalies.1979-2011.850.nc > > > > > On 6/6/12 3:15 AM, juki juki wrote: > > Thanks for reply Dennis. > > > >>What datasets were used? > > > > As in the website > > (http://www.ncl.ucar.edu/Applications/mjoclivar.shtml), I used > > http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.html, for > > example for wind I used daily wind data > > (http://www.esrl.noaa.gov/psd/cgi-bin/db_search/DBListFiles.pl?did=32&tid=34020&vid=666 > > <http://www.esrl.noaa.gov/psd/cgi-bin/db_search/DBListFiles.pl?did=32&tid=34020&vid=666>). > > > > For calculating the anomaly, I just simply modified mjoclivar_2.ncl. > > What I worried is about accumulating the daily data and then calculated > > anomaly from accumulated data. If anomaly is calculated from daily data > > (not from accumulated data) and then daily caclculated anomaly was > > accumulated, of course, the result will be different if it is calculated > > from accumulated data. > > > > All plot for OLR is the same as the example, only for wind is different > > and as we know olr data already accumulated by NOAA, but not the case > > for wind. Please see the code I used below. > > > > Again, thank you. > > > > > > Regards, > > Joeky > > > > -----the code > > > > 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 specifications > > ;------------------------------------------------------------- > > NC = True ; create netCDF? > > PLOT = False ; sample plots? > > > > ymdStrt = 19800101 ; start yyyymmdd > > ymdLast = 20120303 ; last > > yrStrt = ymdStrt/10000 > > yrLast = ymdLast/10000 > > lev = 200 ; desired level > > nhar = 4 ; number of fourier comp > > > > var = "uwnd" ; name of file > > vName = "U" ; name for plots > > fili = var+".day.mean.nc" ; input file > > > > ;======================== > > ; get list of all files and open as "one big file" > > ;======================== > > all_files = systemfunc ("ls ./Reana1/uwnd/uwnd*.nc") > > f = addfiles (all_files, "r") ; note the "s" of addfile > > print(all_files) > > ;======================== > > ; choose how files are combined and read in variable across files > > ;======================== > > ListSetType (f, "cat") ; concatenate or "merge" (default) > > > > if (NC) then > > diro= "./Anom/" ; output dir > > filo= var+".day."+lev+".anomalies."+yrStrt+"-"+yrLast+".nc" ; output file > > end if > > > > if (PLOT) then > > wksType = "ps" > > wksName = "mjoclivar" ; "mjo."+yrStrt+"_"+yrLast > > end if > > > > ;*********************************************************** > > ; Read user specified time and create required yyyyddd > > ;*********************************************************** > > ; f = addfile (diri+fili , "r") > > > > time = f[:]->time ; all times on file > > ymd = cd_calendar(time, -2) ; yyyymmdd > > iStrt = ind(ymd.eq.ymdStrt) ; index start > > iLast = ind(ymd.eq.ymdLast) ; index last > > delete(time) > > delete(ymd) > > > > ;*********************************************************** > > ; Read user specified time and create required yyyyddd > > ;*********************************************************** > > time = f[:]->time(iStrt:iLast) ; time:units = "hours since" > > TIME = cd_calendar(time, 0) ; type float > > year = floattointeger( TIME(:,0) ) > > month = floattointeger( TIME(:,1) ) > > day = floattointeger( TIME(:,2) ) > > ddd = day_of_year(year, month, day) > > yyyyddd = year*1000 + ddd ; needed for input > > > > ;*********************************************************** > > ; Read data: short2flt > > ;*********************************************************** > > x = short2flt(f[:]->$var$(iStrt:iLast,{lev},:,:)) ; convert to float > > printVarSummary(x) > > > > ;*********************************************************** > > ; Compute daily climatology: raw and then 'smoothed' > > ;*********************************************************** > > xClmDay = clmDayTLL(x, yyyyddd) ; daily climatology at each grid point > > printVarSummary(xClmDay) > > > > ;*********************************************************** > > ; Compute smoothed daily climatology using 'nhar' harmonics > > ;*********************************************************** > > xClmDay_sm = smthClmDayTLL(xClmDay, nhar) > > printVarSummary(xClmDay_sm) > > > > ;*********************************************************** > > ; Compute daily anomalies using raw and smoothed climatologies > > ;*********************************************************** > > xAnom = calcDayAnomTLL (x, yyyyddd, xClmDay) > > printVarSummary(xAnom) > > printMinMax(xAnom, True) > > > > xAnom_sm = calcDayAnomTLL (x, yyyyddd, xClmDay_sm) > > xAnom_sm_at_long_name = "Anomalies from Smooth Daily Climatology" > > printVarSummary(xAnom_sm) > > printMinMax(xAnom_sm, True) > > > > delete( x ) ; no longer needed > > > > > > ;*********************************************************** > > ; Create netCDF: convenience use 'simple' method > > ;*********************************************************** > > > > dimx = dimsizes(xAnom) > > ntim = dimx(0) > > nlat = dimx(1) > > mlon = dimx(2) > > > > if (NC) then > > > > lat = f[0]->lat > > lon = f[0]->lon > > > > system("/bin/rm -f "+diro+filo) ; rm any pre-exist file, if any > > fnc = addfile (diro+filo, "c") > > > > filAtt = 0 > > filAtt_at_title = vName+": Daily Anomalies: "+yrStrt+"-"+yrLast > > filAtt_at_source_file = fili > > filAtt_at_creation_date = systemfunc("date") > > fileattdef( fnc, filAtt ) ; copy file attributes > > > > setfileoption(fnc,"DefineMode",True) > > varNC = vName+"_anom" > > varNC_sm = vName+"_anom_sm" > > > > dimNames = (/"time", "lat", "lon"/) > > dimSizes = (/ -1 , nlat, mlon/) > > dimUnlim = (/ True , False, False/) > > filedimdef(fnc,dimNames,dimSizes,dimUnlim) > > > > filevardef(fnc, "time" ,typeof(time),getvardims(time)) > > filevardef(fnc, "lat" ,typeof(lat) ,getvardims(lat)) > > filevardef(fnc, "lon" ,typeof(lon) ,getvardims(lon)) > > filevardef(fnc, varNC ,typeof(xAnom) ,getvardims(xAnom)) > > filevardef(fnc, varNC_sm,typeof(xAnom) ,getvardims(xAnom)) > > > > filevarattdef(fnc,"time" ,time) ; copy time attributes > > filevarattdef(fnc,"lat" ,lat) ; copy lat attributes > > filevarattdef(fnc,"lon" ,lon) ; copy lon attributes > > filevarattdef(fnc,varNC ,xAnom) > > filevarattdef(fnc,varNC_sm,xAnom_sm) > > > > fnc->time = (/time/) > > fnc->lat = (/lat/) > > fnc->lon = (/lon/) > > fnc->$varNC$ = (/xAnom /) > > fnc->$varNC_sm$ = (/xAnom_sm/) > > end if > > > > ;************************************************ > > ; plotting parameters > > ;************************************************ > > if (PLOT) then > > LAT = (/ 60, 45, 5, -5, -45, 60 /) > > LON = (/270, 30, 90, 90, 180, 0 /) > > nPts = dimsizes( LAT ) > > > > plot = new ( nPts, graphic) > > data = new ( (/2,366/), typeof(xClmDay), getFillValue(xClmDay)) > > > > wks = gsn_open_wks (wksType,wksName) > > > > res = True ; plot mods desired > > res_at_gsnDraw = False > > res_at_gsnFrame = False > > res_at_trXMinF = 1 > > res_at_trXMaxF = 366 > > ;res_at_tiMainString = "" > > res_at_xyLineThicknesses = (/1.0, 2.0/) ; make 2nd lines thicker > > res_at_xyLineColors = (/"blue","red"/) ; change line color > > res_at_xyMonoDashPattern = True ; all solid > > > > do np=0,nPts-1 > > data(0,:) = xClmDay(:,{LAT(np)},{LON(np)}) > > data(1,:) = xClmDay_sm(:,{LAT(np)},{LON(np)}) > > res_at_gsnCenterString = "lat="+LAT(np)+" lon="+LON(np) > > plot(np) = gsn_csm_y (wks,data,res) ; create plot > > end do > > resP = True ; modify the panel plot > > resP_at_txString = vName+": Raw/Smooth Daily Climatology: "+yrStrt+"-"+yrLast > > resP_at_gsnMaximize = True > > resP_at_gsnPaperOrientation = "portrait" > > gsn_panel(wks,plot,(/(nPts/2),2/),resP) ; now draw as one plot > > > > ;========== > > ; Plot anomalies for an arbitrarily selected near equatorial location > > ; Time: Oct 1, 1996 to April 1,1997 [arbitrary selection] > > ;========== > > LATX = 0 > > LONX = 90 > > > > yyyymmdd = cd_calendar(time, -2) > > ;;yrfrac = yyyymmdd_to_yyyyfrac (yyyymmdd, 0) > > ;;delete(yrfrac_at_long_name) > > > > xAnom_at_long_name = "Anomalies from Raw" ; short labels for plot > > xAnom_sm_at_long_name = "Anomalies from Smooth" > > > > ntBegin = ind(yyyymmdd.eq.20011001) > > ntEnd = ind(yyyymmdd.eq.20020401) > > monthLabels = (/1,4,7,10/) > > monNam = (/"Jan","Feb","Mar","Apr","May","Jun" \ > > ,"Jul","Aug","Sep","Oct","Nov","Dec" /) > > xVal = new(ntim, typeof(xAnom&time) , "No_FillValue") ; bigger than > > xLab = new(ntim, "string", "No_FillValue") ; needed > > xValm = new(ntim, typeof(xAnom&time) , "No_FillValue") ; bigger than > > > > ntm = -1 > > cr = inttochar(10) ; carriage return > > do nt=ntBegin,ntEnd > > if (day(nt).eq.1) then > > ntm = ntm + 1 > > xVal(ntm) = xAnom&time(nt) > > xLab(ntm) = monNam(month(nt)-1) > > if (month(nt).eq.1) then > > xLab(ntm) = xLab(ntm) + cr +sprinti("%0.4i", year(nt)) > > end if > > end if > > end do > > > > rxy = True > > rxy_at_gsnDraw = False > > rxy_at_gsnFrame = False > > rxy_at_gsnYRefLine = 0.0 ; create a reference line > > rxy_at_gsnAboveYRefLineColor = "red" ; above ref line fill red > > rxy_at_gsnBelowYRefLineColor = "blue" ; below ref line fill blue > > rxy_at_xyLineThicknessF = 2.0 > > rxy_at_vpHeightF = 0.4 ; resize > > rxy_at_vpWidthF = 0.8 > > rxy_at_tmXBMode = "Explicit" > > rxy_at_tmXBValues = xVal(0:ntm) > > rxy_at_tmXBLabels = xLab(0:ntm) > > > > plot(0) = gsn_csm_xy (wks,time(ntBegin:ntEnd) \ > > ,xAnom(ntBegin:ntEnd,{0},{90}),rxy) > > plot(1) = gsn_csm_xy (wks,time(ntBegin:ntEnd) \ > > ,xAnom_sm(ntBegin:ntEnd,{0},{90}),rxy) > > resP_at_txString = vName+": Anomalies: (0,90E)" > > gsn_panel(wks,plot(0:1),(/2,1/),resP) > > > > end if > > end > > > > ------------------------------------------------------------------------ > > *From:* Dennis Shea <shea_at_ucar.edu> > > *To:* juki juki <juky_emc2_at_yahoo.com> > > *Cc:* NCL Talk <ncl-talk_at_ucar.edu> > > *Sent:* Wednesday, June 6, 2012 11:20 AM > > *Subject:* Re: [ncl-talk] Wind anomaly, MJO_PC_INDEX and MJO clivar > > > > The real key is the daily anomalies. If they are not correct, none > > of the MJO diagnostics are correct. You should look very carefully > > at how the anomalies were generated. What datasets were used? > > For debugging and testing, use a smaller temporal sample and > > plot the data. > > > > > > On 6/3/12 10:58 PM, juki juki wrote: > >  > Dear All; > >  > As my previous email that I would like to plot MJO phase diagram by > >  > usinghttp://www.ncl.ucar.edu/Applications/Scripts/mjoclivar_15.ncl. To > >  > see the performance of the code, first I plot same data period as in the > >  > website. The result is very different from the example > >  > (http://www.ncl.ucar.edu/Applications/Images/mjoclivar_15_lg.png) in > >  > which the location of all MJO phases for all months is different from > >  > that in the website. > >  > Thus, my MJO_PC_INDEX.nc <http://MJO_PC_INDEX.nc> > > <http://mjo_pc_index.nc/> may NOT be correct. > >  > As we know the input of mjoclivar_15.ncl is MJO_PC_INDEX.nc > >  > <http://MJO_PC_INDEX.nc> which is creted by > >  > http://www.ncl.ucar.edu/Applications/Scripts/mjoclivar_14.ncl. Moreover, > >  > the input of > >  > http://www.ncl.ucar.edu/Applications/Scripts/mjoclivar_14.ncl are: > >  > > >  > > >  > filolr = "olr.day.anomalies.1980-2012.nc > > <http://olr.day.anomalies.1980-2012.nc> > >  > <http://olr.day.anomalies.1980-2012.nc/>" > >  > filu200 = "uwnd.day.200.anomalies.1980-2012.nc > > <http://uwnd.day.200.anomalies.1980-2012.nc> > >  > <http://uwnd.day.200.anomalies.1980-2012.nc/>" > >  > filu850 = "uwnd.day.850.anomalies.1980-2012.nc > > <http://uwnd.day.850.anomalies.1980-2012.nc> > >  > <http://uwnd.day.850.anomalies.1980-2012.nc/>" > >  > Among these 3 inputs, it seems the wind anomalies are not correct because > >  > when I try to calculate Conventional (covariance) univariate EOF > >  > analysis <http://www.ncl.ucar.edu/Applications/Scripts/mjoclivar_12.ncl> > >  > fortime span 1995-1999 (as in the website) using > >  > http://www.ncl.ucar.edu/Applications/Scripts/mjoclivar_12.ncl, the > >  > result of OLR is close to the website, but for U850 and U200 is > >  > different and significant different for EOF-3 (enclosed). Thus, I think > >  > my wind anomalies is not correctly calculated. > >  > I calculated wind anomalies from daily data as the code below > >  > (enclosed). I don’t know whether the way to merge the files and then > >  > calculate anomaly from those files is correct or not. Probably you have > >  > suggestion and help for that. > >  > Thanks you, > >  > JoeKy > >  > > >  > > >  > _______________________________________________ > >  > 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 > > > <mjoclivar15.ps>_______________________________________________ > 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 Tue Jun 12 18:06:42 2012

This archive was generated by hypermail 2.1.8 : Fri Jun 15 2012 - 14:51:31 MDT