Re: About Station data to gridded data

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Fri, 26 Dec 2008 08:57:34 -0700

Hello

This is offline.

I am attaching a Barnes/Cressman objective analysis
written in NCL.

It is not well tested. It should be similar
to GrADS OACRES. Just like OACRES, it can be slow.

There is no documentaion. You must define the grid onto which
you want the data. The lat/lon of the grids must be monotonically
increasing.

The function you call is:

    zGrid = tripleObjAnalysis(zlon,zlat,z,glon,glat,dcrit,opt)

where zlon[*],zlat[*],z[*] are the observation triplet.

glon[*], glat[*] is the grid

dcrit[*] are the search radii.
; dcrit - 1D array containing successive radii of influence.
; Must be expressed in degrees latitude and should be
; monotonically decreasing. eg: dcrit = (/10, 5, 3/)

I am also attaching a sample test script which you can model
your script after.

Let me know if it works for you.

===
Note also: in the next release of NCL [5.1.0, January 2009],
there is a "bin_avg.ncl" which does binning.

http://www.ncl.ucar.edu/future_release.shtml

Good luck
D

daxiawj wrote:
> Hello,
>
> Does anybody knows how to interpolate station data to gridded data
> in NCL? I've tried
> natgrid and ccsgrid, but it looks like the results are a little strange
> and doesn't like the ones
> interpolated by GrADS. Are there any Functions to do that like in GrADS(
> OACRES or so)?
>
> There is another question, when regridding temperature, the
> influence of topography
> must take in to consideration, and are there any procedure in NCL can
> do that?
>
> Thank you! Merry Christmas and Happy New Year!
>
> Best Wishes!
> Wang Jun
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

;-----------------------------------------------------------------------------------
function tripleObjAnalysis1D(x[*],y[*],z[*],lon[*],lat[*],dcrit[*]:numeric,opt:logical)
; This should *not* be invoked directly by the user.
; It is called by "function objAnalysis"
;
; Nomenclature
; x,y,z - lon,lat,observation triplets
; lat - lat of returned grid. Need not be equally spaced
; but must be monotonically increasing. Should have the
; units attribute assigned: lat_at_units="degrees_north"
; lon - lon of returned grid. Need not be equally spaced
; but must be monotonically increasing. Should have the
; units attribute assigned: lon_at_units="degrees_east"
; dcrit - 1D array containing successive radii of influence.
; Must be expressed in degrees latitude and should be
; monotonically de. eg: dcrit = (/10, 5, 3/)
; opt - variable to which optional attributes are attached
; @guess = 2D guess array [input]
; @timing = True ; print times
; @count => Return number of observations used in each scan
; @nObs ==> (nscan,:,:)
;
local nlat, mlon, nScan, G, dimG, rankG, zonavg, i, j, ij, ns, nl, ml \
    , gcdist, diff, cf, nObs, dc2, nd
begin
  nlat = dimsizes(lat)
  mlon = dimsizes(lon)

  nScan = dimsizes(dcrit)
  nObs = new( (/nScan,nlat,mlon/), "integer", "No_FillValue" )
  nObs = 0

  if (opt .and. isatt(opt,"guess")) then

      G = opt_at_guess ; 1st guess
      dimG = dimsizes(G)
      rankG = dimsizes(dimG)
      if (.not.(rankG.eq.2)) then
          print("objAnalysis1D: rankG="+rankG+" expecting 2D")
          exit
      end if
      if (.not.(nlat*mlon.ne.prod(dimG))) then
          print("objAnalysis1D: dimension sizes of G and nlat*mlon must match")
          print("objAnalysis1D: dimG="+dimG)
          print("objAnalysis1D: nlat="+nlat+" mlon="+mlon)
          exit
      end if

  else

      G = new( (/nlat,mlon/), typeof(z), getFillValue(z) ) ; 1st guess

      if (isatt(opt,"zonal") .and. opt_at_zonal) then ; create zonal avg
          dlat = max(abs(lat(1:nlat-1)-lat(0:nlat-2)) ) ; nominal
         ;dlat = 2*dlat ; expand to get more data for zonal average
                                                         ; bigger range
          zonavg = new( nlat, typeof(z), z@_FillValue)
          do nl=0,nlat-1
             i = ind(y.le.(lat(nl)+dlat) .and. y.ge.(lat(nl)-dlat))
             if (.not.all(ismissing(i))) then
                 zonavg(nl) = avg( z(i) ) ; zonal avg of all observations
             end if
             delete(i)
          end do

          if (any(ismissing(zonavg))) then
              zonavg = linmsg(zonavg, -1) ; linearly interpolate
          end if
         print("objAnalysis1D: lat="+lat+" zonavg="+zonavg)
                                               ; arbitrary smooth
        ;;zonavg = wgt_runave(zonavg, (/0.25, 0.50, 0.25/), 1)
          zonavg = wgt_runave(zonavg, filwgts_normal (7, 1.0, 0) , 1)

          do nl=0,nlat-1
             G(nl,:) = zonavg(nl)
          end do

          delete(zonavg)
      else
          G = 0.0 ; direct ... no 1st guess
      end if
  end if

  G!0 = "lat"
  G!1 = "lon"
  G&lat = lat
  G&lon = lon

  wcStrt = systemfunc("date")
  do ns=0,nScan-1
     nsStrt = systemfunc("date")
     dc2 = dcrit(ns)^2

    do nl=0,nlat-1
       i = ind( abs(y-lat(nl)).le.dcrit(ns))
       if (.not.any(ismissing(i)) ) then
           do ml=0,mlon-1
              if (ns.eq.0 .or. (ns.gt.0 .and. nObs(ns-1,nl,ml).gt.0)) then
                  gcdist = gc_latlon(lat(nl),lon(ml), y(i),x(i), 0,2)
                  nd = num(gcdist.le.dcrit(ns))
                  nObs(ns,nl,ml) = nd ; # observations within radius
     
                  if (nd.gt.0) then
                      j = ind(gcdist.le.dcrit(ns))
                      ij = i(j)
                      diff = z(ij)-G(nl,ml) ; normally interpolate G to z but .....
                      wgt = exp(-4*gcdist(j)^2/dc2)
                      cf = sum(wgt*diff)/sum(wgt) ; correction factor
                     ;print("ns="+ns+" nl="+nl+" ml="+ml+" cf="+cf)
         
                      G(nl,ml) = G(nl,ml) + cf ; update Guess
                      
                      delete(j)
                      delete(ij)
                      delete(cf)
                      delete(wgt)
                      delete(diff)
                  end if ; nd
                  delete(gcdist)
              end if
           end do ; ml
       end if
       delete(i)
    end do ; nl
                                     ; default is to smooth
    if (.not.isatt(opt,"smooth") .or. opt_at_smooth) then
        if (ns.lt.(nScan-1)) then
            G = smth9(G, 0.50,-0.25, 0) ; light local smoother
        else
            G = smth9(G, 0.50, 0.25, 0) ; heavy local smoother
        end if
    end if

    if (isatt(opt,"timing")) then
        wallClockElapseTime(nsStrt,"objAnalysis1D: ns="+ns, 0)
    end if
  end do ; ns`

  if (isatt(opt,"timing")) then
      wallClockElapseTime(wcStrt,"Total time: objAnalysis1D: nScan="+nScan , 0)
  end if

  if (isatt(opt,"count")) then
      G_at_nObs = nObs
  end if

  return(G)
end
;-----------------------------------------------------------------------------------
function tripleObjAnalysis2D(x[*],y[*],z[*],lon2d[*][*],lat2d[*][*],dcrit[*]:numeric,opt:logical)
;
; This should not be invoked directly by the user.
; It is called by "function objAnalysis"
;
; Nomenclature
; x,y,z - lon,lat,observation triplets
; lat - lat of returned grid. Need not be equally spaced
; but must be monotonically increasing. Should have the
; units attribute assigned: lat_at_units="degrees_north"
; lon - lon of returned grid. Need not be equally spaced
; but must be monotonically increasing. Should have the
; units attribute assigned: lon_at_units="degrees_east"
; dcrit - 1D array containing successive radii of influence.
; Must be expressed in degrees latitude and should be
; monotonically de. eg: dcrit = (/10, 5, 3/)
; opt - variable to which optional attributes are attached
; @guess = 2D guess array
; @timing = print times
; @count = number of observation used in each scan
; Return @nObs ==> (nscan,:,:)
;
local nlat, mlon, nScan, G, dimG, rankG, zonavg, i, j, ij, ns, nl, ml \
    , gcdist, diff, cf, nObs, dc2, nd, dimlat, ranklt, LAT, LON, G1D
begin
  dimlat = dimsizes(lat2d)
  nlat = dimlat(0)
  mlon = dimlat(1)

  nScan = dimsizes(dcrit)
  nObs = new( (/nScan,nlat,mlon/), "integer", "No_FillValue" )
  nObs = 0

  if (opt .and. isatt(opt,"guess")) then

      G = opt_at_guess ; 1st guess
      dimG = dimsizes(G)
      rankG = dimsizes(dimG)
      ranklt = dimsizes(dimlat)
      if (.not.(rankG.eq.ranklt)) then
          print("objAnalysis2D: rankG="+rankG+" ranklt="+ranklt)
          exit
      end if
      if (.not.all(dimlat.eq.dimG)) then
          print("objAnalysis2D: all dimension sizes must be the same")
          print("objAnalysis2D: dimltt="+dimlat+" dimG="+dimG)
          exit
      end if

  else

      G = new( dimlat, typeof(z), getFillValue(z) ) ; 1st guess

      if (isatt(opt,"zonal") .and. opt_at_zonal) then ; create zonal avg
          mnlat = min(lat2d)
          mxlat = max(lat2d)
          lat1d = fspan(mnlat,mxlat,nlat) ; nominal
          dlat = (mxlat-mnlat)/(nlat-1) ; nominal
          dlat = 4*dlat ; expand to get more data for zonal average
                                                         ; bigger range
          zonavg = new( nlat, typeof(z), z@_FillValue)
          do nl=0,nlat-1
             i = ind(y.le.(lat1d(nl)+dlat) .and. y.ge.(lat1d(nl)-dlat))
             if (.not.all(ismissing(i))) then
                 zonavg(nl) = avg( z(i) ) ; zonal avg of all observations
             end if
             delete(i)
          end do

          if (any(ismissing(zonavg))) then
              zonavg = linmsg(zonavg, -1) ; linearly interpolate
          end if

        ;;zonavg = wgt_runave(zonavg, (/0.25, 0.50, 0.25/), 1) ;smooth
          zonavg = wgt_runave(zonavg, filwgts_normal (7, 1.0, 0) , 1) ;smooth

          LAT = ndtooned(lat2d)
          G1D = ndtooned(G)

          do nl=0,nlat-1
             i = ind(LAT.le.(lat1d(nl)+dlat) .and. LAT.ge.(lat1d(nl)-dlat))
             G1D(i) = zonavg(nl)
             delete(i)
          end do

          delete(zonavg)
          delete(lat1d)
          delete(LAT)

          G = onedtond(G1D, dimlat)
          delete(G1D)
      else
          G = 0.0 ; direct ... no 1st guess
      end if
  end if

  G_at_lat2d = lat2d
  G_at_lon2d = lon2d

 ;plotGrid("ps", "GUESS", "1st GUESS: ZONAL", G, lat2d, lon2d)

  LAT = ndtooned( lat2d )
  LON = ndtooned( lon2d )

  wcStrt = systemfunc("date")
  do ns=0,nScan-1
     nsStrt = systemfunc("date")
     dc2 = dcrit(ns)^2

    do nl=0,nlat-1
      do ml=0,mlon-1
         if (ns.eq.0 .or. (ns.gt.0 .and. nObs(ns-1,nl,ml).gt.0)) then
             i = ind( abs(y-lat2d(nl,ml)).le.dcrit(ns) )
             if (.not.any(ismissing(i)) ) then
                 gcdist = gc_latlon(lat2d(nl,ml),lon2d(nl,ml), y(i),x(i), 0,2)
    
                 nd = num(gcdist.le.dcrit(ns))
                 nObs(ns,nl,ml) = nd ; # observations within radius
    
                 if (nd.gt.0) then
                     j = ind(gcdist.le.dcrit(ns))
                     ij = i(j)
                     diff = z(ij)-G(nl,ml) ; normally interpolate G to z but .....
                     wgt = exp(-4*gcdist(j)^2/dc2)
                     cf = sum(wgt*diff)/sum(wgt)
                    ;print("ns="+ns+" nl="+nl+" ml="+ml+" cf="+cf)
        
                     G(nl,ml) = G(nl,ml) + cf
                     
                     delete(j)
                     delete(ij)
                     delete(cf)
                     delete(wgt)
                     delete(diff)
                 end if
                 delete(gcdist)
             end if
             delete(i)
         end if
      end do ; ml
    end do ; nl
                                     ; default is to smooth
    if (.not.isatt(opt,"smooth") .or. opt_at_smooth) then
        if (ns.lt.(nScan-1)) then
            G = smth9(G, 0.50,-0.25, 0) ; light local smoother
        else
            G = smth9(G, 0.50, 0.25, 0) ; heavy local smoother
        end if
    end if

    if (isatt(opt,"timing")) then
        wallClockElapseTime(nsStrt,"objAnalysis2D: ns="+ns , 0)
    end if
  end do ; ns`

  if (isatt(opt,"count")) then
      G_at_nObs = nObs
  end if

  if (isatt(opt,"timing")) then
      wallClockElapseTime(wcStrt,"Total time: objAnalysis2D: nScan="+nScan , 0)
  end if

  return(G)
end
;----------------------------------------------------------------------------
function tripleObjAnalysis(X[*],Y[*],Z[*],lon,lat,dcrit[*]:numeric,opt:logical)
; Perform Barnes [ Cressman ] type iterative correction objective analysis
;
; Nomenclature
; x,y,z - lon,lat,observation triplets
; lat - lat of returned grid. Need not be equally spaced
; but must be monotonically increasing. Should have the
; units attribute assigned: lat_at_units="degrees_north"
; lon - lon of returned grid. Need not be equally spaced
; but must be monotonically increasing. Should have the
; units attribute assigned: lon_at_units="degrees_east"
; dcrit - 1D array containing successive radii of influence.
; Must be expressed in degrees latitude and should be
; monotonically de. eg: dcrit = (/10, 5, 3/)
; dims = dimsizes(dcrit) , nscan = dims(0)
; opt - variable to which optional attributes are attached
; @guess = 2D guess array [input]
; @timing = print times
; @count = number of observation used in each scan
; Return @nObs ==> (nscan,:,:)
;
local wcStrt, i, XX, YY, ZZ, k, x, y, z, dimLat, dimLon, rankLat, rankLon
begin
  wcStrt = systemfunc("date")

  if (.not.isatt(Z,"_FillValue")) then
      Z@_FillVlaue = 1e20
  end if
                            ; eliminate missing values
  i = ind(.not.ismissing(Z))
  XX = X(i)
  YY = Y(i)
  ZZ = Z(i)
  delete(i)
                           
  k = dim_pqsort(YY, 1) ; sort obs in ascending latitude order
  x = XX(k) ; not used here ... too lazy to change code
  y = YY(k)
  z = ZZ(k)
  delete(k)

  delete(XX)
  delete(YY)
  delete(ZZ)

 ;print("objAnalysis: z="+z+" lon="+x+" lat="+y )

  dimLat = dimsizes(lat)
  dimLon = dimsizes(lon)
  rankLat = dimsizes(dimLat)
  rankLon = dimsizes(dimLon)
  if (rankLat.ne.rankLon) then
      print("objAnalysis: ranks of lat and lon must match")
      print("objAnalysis: rankLat="+rankLat)
      print("objAnalysis: rankLon="+rankLon)
  end if
  if (rankLat.gt.2) then
      print("objAnalysis: ranks of lat and lon must be 1 or 2")
      print("objAnalysis: rankLat="+rankLat)
  end if

  if (rankLat.eq.1) then
      zGrid = tripleObjAnalysis1D(X,Y,Z,lon,lat,dcrit,opt)
  else
      zGrid = tripleObjAnalysis2D(X,Y,Z,lon,lat,dcrit,opt)
  end if
  if (isatt(opt,"timing")) then
      wallClockElapseTime(wcStrt,"objAnalysis", 0)
  end if

  if (isatt(Z,"long_name")) then
      zGrid_at_long_name = Z_at_long_name
  end if
  if (isatt(Z,"units")) then
      zGrid_at_units = Z_at_units
  end if
 
  return(zGrid)
end

;
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "./CAS_generic.ncl"

begin
;
; Data is stored in four columns: station_name lat lon pwv
;
  ncol = 4
  fname = "temp.USA.tst_objAnalysis_1.txt"
  npts = numAsciiRow (fname )
  table = asciiread(fname, (/npts,ncol/), "float")

  id = floattointeger( table(:,0) )
  LON = table(:,1)
  LAT = table(:,2)
  z = table(:,3)
;;print(LON+" "+LAT+" "+z)

  latS = 0
  latN = 74
  lonL = -164
  lonR = -60

  pltType = "x11"
  pltName = "tst_objAnalysis_1D_temp"
  wks = gsn_open_wks(pltType, pltName)
  gsn_define_colormap(wks,"so4_23")

  res = True
  res_at_gsnMaximize = True ; make as large as possible
  res_at_gsnSpreadColors = True
  res_at_gsnSpreadColorStart = 2
  res_at_gsnSpreadColorEnd = -1

  res_at_gsnPaperOrientation = "portrait"
  res_at_gsnAddCyclic = False ; data not cyclic

  res_at_cnFillOn = True ; color fill
;;res_at_cnFillMode = "RasterFill" ; Raster Mode
  res_at_cnLinesOn = False ; True is default
  res_at_cnLineLabelsOn = False ; True is default
  res_at_lbLabelAutoStride = True

 ;res_at_mpProjection = "CylindricalEquidistant"
  res_at_mpMinLatF = latS ; min((/min(y)-1 , 0/))
  res_at_mpMaxLatF = latN ; max((/max(y)+1 , 90/))
  res_at_mpMinLonF = lonL ; min(x)-1
  res_at_mpMaxLonF = lonR ; max(x)+1
  res_at_mpFillOn = False ; True is default for gsn_csm

  res_at_gsnFrame = False ; So we can draw markers
;
; Draw markers on the plot in the lat/lon locations.
;
SKIP = True
if (.not.SKIP) then
  plot = gsn_csm_map_ce(wks,res) ; draw map and contour
  mkres = True
  mkres_at_gsMarkerIndex = 17 ; Filled circle
  mkres_at_gsMarkerSizeF = 0.01
  gsn_polymarker(wks,plot,LON,LAT,mkres)
  frame(wks) ; Now advance the frame.
end if

  res_at_gsnFrame = True ; reset to default

; =============== OBJECTIVE ANALYSIS =======================

  dlat = 0.50
  dlon = 0.50
  nlat = floattointeger((latN-latS)/dlat + 1 )
  mlon = floattointeger((lonR-lonL)/dlon + 1 )
  lat = fspan(latS,latN,nlat)
  lon = fspan(lonL,lonR,mlon)
  lat!0 = "lat"
  lon!0 = "lon"
  lat_at_units = "degrees_north"
  lon_at_units = "degrees_east"
  lat&lat = lat
  lon&lon = lon

  dcrit = (/10,7,4,1/)
  nscan = dimsizes(dcrit)
  opt = True
 ;GUESS = new ( (/nlat,mlon/), typef(z), 1e20)
 ;GUESS = 0.0
 ;opt_at_guess = GUESS
  opt_at_zonal = True
  opt_at_count = True
  opt_at_timing = True
  opt_at_smooth = True

  zGrid = tripleObjAnalysis(LON,LAT,z,lon,lat,dcrit,opt) ; negative values
  printVarSummary(zGrid)
  printMinMax(zGrid, True)

  nObs = zGrid_at_nObs
  do ns=0,nscan-1
     zGrid = where(nObs(ns,:,:).eq.0, zGrid@_FillValue, zGrid)
     res_at_tiMainString = "BHALL data: ns="+ns+" dcrit="+dcrit(ns)
     plot = gsn_csm_contour_map_ce(wks,zGrid,res)
  end do

  
  option = True ; method=1 and domain=0 are set
  option_at_distmx = 222 ; kilometers
  ZGRID = triple2grid(LON,LAT,z, lon,lat, option)
  ZGRID!0 = "lat"
  ZGRID!1 = "lon"
  ZGRID&lat = lat
  ZGRID&lon = lon
  ZGRID@_FillValue = 0.0 ; make transparent
  ZGRID@_FillValue = 1.e20
  res_at_tiMainString = "USA: triple2grid"
  plot = gsn_csm_contour_map_ce(wks,ZGRID,res) ; draw map and contour

;;;;; USE poisson_grid_fill ***********

end
Received on Fri Dec 26 2008 - 08:57:34 MST

This archive was generated by hypermail 2.2.0 : Wed Dec 31 2008 - 10:52:54 MST