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