Re: Interpolation Anomaly

Quite frankly, the 2nd code is hard to follow.
There are so many 'reorderings' it is confusing.

The "zon.ncl" creates a plot of level vs time of zonal averages.

The 2nd script does

==
; take zonal average
  uavg = dim_avg_Wrap(U2) ; (time,lev_p,lat)
==
Then you do a sequence of reorders, which make things obscure.

==
  ulatav = dim_avg_Wrap(uavg) ; (time,lev_p)
  ulatav_reorder = ulatav(lev_p | :, time | :) ; (lev_p,time)
  work = ulatav_reorder ;
(lev_p,time)
   uavg_reorder = uavg(lev_p | :, time | :, lat | :) ; (lev_p,time,lat)
====
later you
====
    work(lev_p|:,time|:) = uavg_reorder(:,:,nl) ;
(lev_p,time) ???????????
     plot = gsn_csm_pres_hgt(wks, work, res )
This is *very* confusing.
==============

Why not just the following
===
; take zonal average
  uavg = dim_avg_Wrap(U2) ; (time,lev_p,lat)

  wks = ...
  res = ...
  plot = gsn_csm_pres_hgt(wks, uavg(lev_p|:, time|:, lat|:nl) , res )

Good luck

Helen Parish wrote:
>
> I wanted to make a height versus time type plot for the zonally
> averaged zonal wind. Data for the 44 different times is stored in
> separate netcdf files (44 total).
>
> In the first case I simply plotted the zonally averaged zonal wind as
> a function of time in years on the x-axis, and model hybrid level on
> the y-axis. (Note that this is not the Earth - there are 50 levels
> total, and the pressure is 92 bars at the surface). (Also note that
> the y-axis label is produced automatically). Since the highest
> altitude level is level 0, and level 49 is at the surface, this simple
> plot is effectively upside down, with level 0, representing the
> highest altitude, at the bottom of the plot. The y-axis is also not
> directly scaled according to pressure.
>
> In the second case, I have tried to plot exactly the same results, as
> a function of time in years on the x-axis, but with the y-axis values
> given in terms of pressure / height. In order to do this I needed to
> interpolate the data from hybrid coordinates onto pressure
> coordinates, using NCL function vinth2p, for all 44 files. In this
> case the plot is the right way up, with the surface at the bottom of
> the plot.
>
> The two plots below show the non-interpolated, and the interpolated
> results respectively. The non-interpolated results look relatively
> smooth, but the interpolated results show jagged structures, which do
> not appear to have been in the original data. Also, I am not sure if
> corresponding features are still at the same times as they were in the
> non-interpolated version. This could either be due to something that
> happens in the interpolation process, or perhaps because I have done
> something incorrect in doing the interpolation.
>
> If it is a problem with how the interpolation works, perhaps there are
> some parameters I can adjust that may affect the way the results look
> ?. I am sure that these jagged features should not be there.
>
> I also include the two versions of the scripts I used. The script
> "zon.ncl" is the non-interpolated version. The script "zoninterp.ncl"
> it the interpolated version.
>
> Can anyone see what is going wrong here ?.
>
> Thanks,
> Helen.
>
>
> *a). NON-INTERPOLATED PLOT :*
>
> =
> ------------------------------------------------------------------------
>
>
> *b). INTERPOLATED PLOT :*
> =
> ------------------------------------------------------------------------
>
> *c). SCRIPT TO MAKE THE NON-INTERPOLATED PLOT :*
>
>
> ;***********************
> ; zon.ncl
> ;***********************
> 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
>
> ; Keep all files open when have a lot of files :
> setfileoption("nc","SuppressClose",False)
>
> diri = "./"
> fili = systemfunc("cd "+diri+" ; ls zeusrand09125*.cam2.h0*nc")
>
> f = addfiles (fili,"r")
> U = addfiles_GetVar(f,fili,"U")
> lat = addfiles_GetVar(f,fili,"lat")
> lon = addfiles_GetVar(f,fili,"lon")
> lev = addfiles_GetVar(f,fili,"lev")
> time = addfiles_GetVar(f,fili,"time")
> printVarSummary( U )
>
> dimu = dimsizes( U )
> ntim = dimu(0)
> klvl = dimu(1)
> nlon = dimu(2)
> mlon = dimu(3)
>
> uavg = dim_avg_Wrap(U)
>
> uavg_reverse = uavg(lev | :, time | :, lat | :)
>
> ;***********************
> ; Create Plot
> ;***********************
> wks = gsn_open_wks ("pdf", "zon" ) ; open workstation
> gsn_define_colormap(wks,"rainbow") ; choose colormap
>
> res = True
> res_at_cnFillOn = True
> res_at_lbLabelAutoStride = True
> res_at_gsnMaximize = True ; if [ps, eps, pdf] make large
> res_at_gsnSpreadColors = True ; span color map
>
> nl = 96 ; equator
>
> res_at_gsnLeftString = ""
> res_at_gsnCenterString = "Zon. av. Zonal Wind, 0 deg lat, years 155-198"
> res_at_gsnRightString = ""
>
> res_at_tiXAxisString = "Time (years)"
>
> plot = gsn_csm_contour(wks, uavg_reverse(:,:,nl), res ) ;
> (time,lev,lat)
>
> res_at_trYReverse = True
>
> end
>
>
>
> *d). SCRIPT TO MAKE INTERPOLATED PLOT :*
>
>
> ;***********************
> ; zoninterp.ncl
> ;***********************
> 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
>
> ; Keep all files open when have a lot of files :
> setfileoption("nc","SuppressClose",False)
>
> diri = "./"
> fili = systemfunc("cd "+diri+" ; ls zeusrand09125*.cam2.h0*nc")
> print(fili)
> nfili = dimsizes(fili)
>
> lev_p = (/ 91500., 90000., 87500., 80000., 67500., 50000., 30000.,
> 20000., 10000., 5000., 3000., 1500., 750., 390., 200., 100., 50.,25.,
> 14., 7., 3.5, 1.9, 0.95, 0.4, 0.1, 0.03 /)
>
> ; ***********************************************
> ; Read in multiple files, interpolate, and write to temporary output files
> ; **********************************************
> do nf = 0, nfili-1
>
> fi = addfile (diri+fili(nf), "r")
>
> diro = "./"
> filo = "uthelen"+nf+".nc"
> print(filo)
> system ("/bin/rm -f "+diro+filo) ; remove any pre-exist file
> fo = addfile (diro+filo, "c")
>
> hyam = fi->hyam ; read hybrid info
> hybm = fi->hybm
> hyai = fi->hyai ; read hybrid info
> hybi = fi->hybi
> PS = fi->PS
> P0mb = 0.01*fi->P0
> U = fi->U
> lat = fi->lat
>
> lev_p!0 = "lev_p" ; variable and dimension name
> the same
> lev_p&lev_p = lev_p ; create coordinate variable
> lev_p_at_long_name = "pressure" ; attach some attributes
> lev_p_at_units = "hPa"
> lev_p_at_positive = "down"
>
> Up = vinth2p (U, hyam, hybm, lev_p ,PS, 1, P0mb, 2, False)
> copy_VarAtts(U, Up)
> fo->U = Up ; write to netCDF file
>
> end do
> ; **************************
> ; calculate zonal average of zonal wind versus time for multiple
> interpolated files
> ; *************************
> diri = "./"
> fili2 = systemfunc("cd "+diri+" ; ls uthelen*nc")
> print(fili2)
> nfili = dimsizes(fili2)
>
> f2 = addfiles (fili2,"r")
> U2= addfiles_GetVar(f2,fili2,"U")
>
> printVarSummary( U2 ) ; (time, lev_p, lat,lon)
>
> dimt2 = dimsizes( U2 )
> ntim2 = dimt2(0)
> klvl2 = dimt2(1)
> nlat2 = dimt2(2)
> mlon2 = dimt2(3)
>
> ; take zonal average
>
> uavg = dim_avg_Wrap(U2) ; (time,lev_p,lat)
> printVarSummary( uavg )
>
> ; create variable "work" of appropriate size and with appropriate
> variable names
>
> ulatav = dim_avg_Wrap(uavg) ; (time,lev_p)
> printVarSummary( ulatav )
> ulatav_reorder = ulatav(lev_p | :, time | :) ; (lev_p,time)
> printVarSummary( ulatav_reorder )
> work = ulatav_reorder
> printVarSummary( work )
>
> ; reorder data to plot with gsn_csm_pres_hgt
>
> uavg_reorder = uavg(lev_p | :, time | :, lat | :) ; (lev_p,lat,time)
> printVarSummary( uavg_reorder )
>
> ;***********************
> ; Create Plot
> ;***********************
> wks = gsn_open_wks ("pdf", "zoninterp" ) ; open workstation
> gsn_define_colormap(wks,"rainbow") ; choose colormap
>
> res = True
> res_at_cnFillOn = True
> res_at_lbLabelAutoStride = True
> res_at_gsnMaximize = True ; if [ps, eps, pdf] make large
> res_at_gsnSpreadColors = True ; span color map
>
> nl = 96 ; equator
>
> work(lev_p|:,time|:) = uavg_reorder(:,:,nl) ; (lev_p,time)
> printVarSummary( work )
>
> res_at_gsnLeftString = ""
> res_at_gsnCenterString = "Zon av. ZW, 0 deg lat, years 155-198"
> res_at_gsnRightString = ""
>
> res_at_tiXAxisString = "Time (years)"
>
> plot = gsn_csm_pres_hgt(wks, work, res ) ; (lev,time)
>
> res_at_trYReverse = True
>
> end
>
> =
> ------------------------------------------------------------------------
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

-- 
======================================================
Dennis J. Shea                  tel: 303-497-1361    |
P.O. Box 3000                   fax: 303-497-1333    |
Climate Analysis Section                             |
Climate & Global Dynamics Div.                       |
National Center for Atmospheric Research             |
Boulder, CO  80307                                   |
USA                        email: shea 'at' ucar.edu |
======================================================
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk