Re: Interpolation Anomaly

Thank you Dennis, for your suggestion to improve the
comprehensibility of the code. I have made the adjustments you
suggested. However the more complicated approach must have worked,
as the results from running with the simplified code look exactly the
same, as shown below.

The code performs a relatively complicated process in order to read
in the data from multiple files, interpolate the data in each file,
write out to multiple output files, and then recalculate the zonal
average from those multiple output files. I am wondering if I have
made a mistake somewhere in this interpolation process that might be
causing the anomalous features in the interpolated results, which
were not there before the interpolation ?.

If I have not made a mistake in the interpolation process, then I
presume that these features appear due to the interpolation itself.
If so, I am wondering if there is anything that can be adjusted,
perhaps in the interpolation procedure, or in the spline fitting that
produces the graphical output, or whether there is some smoothing
algorithm that can be applied, to make the interpolated results look
more like the original un-interpolated results ?.

I also include the script I am currently using below.

Thanks,
Helen.

;***********************
; 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 Plot
;***********************
  wks = gsn_open_wks ("pdf", "zoninterp.
165to175.lat0" ) ; 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. ZW, 0 deg lat, years 165-175"
      res_at_gsnRightString = ""

      res_at_tiXAxisString = "Time (years)"

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

      res_at_trYReverse = True

end

a. Interpolated
results b.
Original un-interpolated results
                                                                         
                           (on y-axis levels grid, which is inverted
                                                                         
                                         top to bottom)

On Apr 16, 2009, at 4:20 PM, Dennis Shea wrote:

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

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk