# Re: Contouring "strangely" gridded data

```----------
X-Sun-Data-Type: text
X-Sun-Data-Description: text
X-Sun-Data-Name: text
X-Sun-Charset: us-ascii
X-Sun-Content-Lines: 59

Hi Mark,

I've never thought too much of the pplus gridding routine.
It is ok for a quick look, but its probably better to do
the gridding yourself, then feed ferret an already-gridded
field. The pplus routine is a subroutine called zgrid (I
have the code if you're interested). It was written in 1968,
so it's a nice example of very early fortran, but neither
very efficient or very intelligent. You will find that you get
*very* different results based on how you choose the input
parameters nx,ny.

Myself, I like gridding using running weighted gaussians.
The attached subroutines do this in either 2 or 3 dimensions.
The idea is that the gridded value at (x0,y0,t0) is the weighted
sum of nearby points, with weights given by exp{((xn-x0)/XX)**2},
where xn and x0 are the x-positions of data point n and the grid
point and XX is the chosen mapping scale in x. Similar calculations
are made for y and t. Each data point is part of the weighted sum
for several gridpoints, weighted according to its distance from each.
If V_n are the original data values in a list n=1,Np, then the gridded
value V' at gridpoint (x0,y0,t0) is:

Sum(n=1,Np) V_n*W_n
V'(x0,y0,t0) =   -------------------			(1)
Sum(n=1,Np) W_n

where Np is the number of data points and W_n are the weights in
3 dimensions described above. In practice one cuts off Np to some
"nearby" region (in these subroutines that is hardwired to be the
point where W_n=exp{-4}=0.018).

The big advantage of this method is that you can choose the mapping
scales XX,YY,TT appropriate for your data.

These routines also allow time-wrapping (to map data onto a
climatological year, say). This is flagged by subroutine argument
iwflag.

To use these routines, have a loop in your main program that loops
over the data values. Call the subroutine gausswt (or gausswt_2d
for 2-d mapping) once for each data value. These subroutines
accumulate the sums in (1). Then once all data points have been
summed, close that loop and call subroutine gaussfin to do the
final division in (1). I hope the routines are adequately documented.

Let me know if you have any problems. I hope there're no bugs in this
code.

This method is referenceable: Kessler, W.S. and J.P. McCreary, 1993.
The annual wind-driven Rossby wave in the subthermocline equatorial
Pacific. JPO, 23, 1192-1207.

Regards, Billy K

PS: if the data are very densely packed (as it appears from your message
to Ed), then it shouldn't make any difference what gridding method you
use.
----------
X-Sun-Data-Type: default
X-Sun-Data-Description: default
X-Sun-Data-Name: gaussmap-subs.f
X-Sun-Charset: us-ascii
X-Sun-Content-Lines: 161

c.............subroutines to do gaussian-weighted mapping onto grids.

c************************************************************************
c............sub gausswt forms the weight sums inside loop over all data.
c............--->>> 3-d mapping (x,y,t)
c............--->>> allows wraparound year (flagged by arg iwflag)
c............method is to call this sub for each data value
c............  sub loops on grid locations, maps each irregular data point to
c............  all possible gridpts, weighting by 3-d distance from gridpt.
c............  all calcs done in gridbox units
c............xx/x1/xf/xsc all in same units
c............yy/y1/yf/ysc all in same units
c............tt/t1/tf/tsc all in same units (but nm can be anything)
c............   note that, ie, t1 is the center of gridbox 1. So if dates
c............   are in months, Jan 15=1, Dec 15=12, and Jan 1=0.5, Dec 31=12.5.

c  i	xx,yy,tt=x/y/t location of data pt (data units)
c  i	val=data value at this (xx,yy,tt)
c  o	grid(nx,ny,nm)=sum of weighted values
c  o	wate(nx,ny,nm)=sum of weights
c  i	nx,ny,nm=size of grids
c  i	x1,y1,t1=west/south/date edge of grid (center of 1st box in data units)
c  i	xf,yf,tf=east/north/date edge of grid (center of final box)
c  i	xsc,ysc,tsc=mapping scales (data units)
c  i	iwflag=1 for time wrapping; 0 for no wrapping
c--------------------------------------------------------------------------

subroutine gausswt (xx,yy,tt,val,grid,wate,nx,ny,nm,x1,y1,t1,
.			xf,yf,tf,xsc,ysc,tsc,iwflag)

real grid(nx,ny,nm), wate(nx,ny,nm)

dx=(xf-x1)/real(nx-1)		  ! gridbox sizes in data units
dy=(yf-y1)/real(ny-1)
dt=(tf-t1)/real(nm-1)

xxg=(xx-x1)/dx+1.		  ! grid values of data location
yyg=(yy-y1)/dy+1.
ttg=(tt-t1)/dt+1.

cutoff=2.			  ! cutoff limits search (min wt=e**-4)
xcut=cutoff*xsc/dx		  ! cutoffs scaled to grid units
ycut=cutoff*ysc/dy		  ! look only twice the scale width
tcut=cutoff*tsc/dt		  !   from the gridbox center

do 100 i=1,nx			  ! loop on x gridpoints
xgp=real(i)			  ! center of gridbox
delx=abs(xgp-xxg)		  ! distance of data pt from grid ctr
if (delx.gt.xcut) go to 100  	  ! only do nearby points

do 100 j=1,ny			  ! loop on y gridpoints, same procedure
ygp=real(j)
dely=abs(ygp-yyg)
if (dely.gt.ycut) go to 100

do 100 m=1,nm			  ! loop on t gridpoints, same procedure
tgp=real(m)
delt=abs(tgp-ttg)
if (delt.gt.tcut .and. iwflag.eq.1)
.		delt=abs(delt-real(nm)) 	! allow flagged time wrapping
if (delt.gt.tcut) go to 100

xgas=(delx*dx/xsc)**2		  	! make gaussian exponents
ygas=(dely*dy/ysc)**2
tgas=(delt*dt/tsc)**2
expn=exp(-xgas-ygas-tgas)		! make the gaussian weight
wate(i,j,m)=wate(i,j,m)+expn		! sum the weights
grid(i,j,m)=grid(i,j,m)+val*expn	! sum the weighted values

100	continue

return
end

c************************************************************************
c............sub gaussfin divides weighted values by sum of weights
c............call this outside of summation loop.
c............for 2-d mapping, call gaussfin with nm=1.

subroutine gaussfin (nx,ny,nm,grid,wate)

real grid(nx,ny,nm),wate(nx,ny,nm)

do 100 i=1,nx
do 100 j=1,ny
do 100 m=1,nm
if (wate(i,j,m).gt.0.) then
grid(i,j,m)=grid(i,j,m)/wate(i,j,m)
else
endif
100	continue

return
end

c************************************************************************
c************************************************************************
c............sub gausswt-2d forms the weight sums inside loop over all data.
c............--->>> 2-d mapping (x,y)
c............--->>> allows wraparound year (flagged by arg iwflag) (1st index)
c............method is to call this sub for each data value
c............  sub loops on grid locations, maps each irregular data point to
c............  all possible gridpts, weighting by 3-d distance from gridpt.
c............  all calcs done in gridbox units
c............xx/x1/xf/xsc all in same units
c............yy/y1/yf/ysc all in same units
c............after sum loop, use gaussfin to finish. Call gaussfin with nm=1.

c  i	xx,yy=x/y location of data pt (data units)
c  i	val=data value at this (xx,yy)
c  o	grid(nx,ny)=sum of weighted values
c  o	wate(nx,ny)=sum of weights
c  i	nx,ny=size of grids
c  i	x1,y1=west/south edge of grid (center of 1st box in data units)
c  i	xf,yf=east/north edge of grid (center of final box)
c  i	xsc,ysc=mapping scales (data units)
c  i	iwflag=1 for time wrapping; 0 for no wrapping
c--------------------------------------------------------------------------

subroutine gausswt_2d (xx,yy,val,grid,wate,nx,ny,x1,y1,
.			xf,yf,xsc,ysc,iwflag)

real grid(nx,ny), wate(nx,ny)

dx=(xf-x1)/real(nx-1)		  ! x-grid size in data units
dy=(yf-y1)/real(ny-1)

xxg=(xx-x1)/dx+1.		  ! grid values of data location
yyg=(yy-y1)/dy+1.

cutoff=2.		 	  ! cutoff to limit search
xcut=cutoff*xsc/dx		  ! look only twice the scale width
ycut=cutoff*ysc/dy		  !   from the gridbox center

do 100 i=1,nx			  ! loop on x gridpoints
xgp=real(i)			  ! center of gridbox
delx=abs(xgp-xxg)		  ! distance of data pt from grid ctr
if (delx.gt.xcut .and. iwflag.eq.1)
.		delx=abs(delt-real(nx))   ! allow flagged time wrapping
if (delx.gt.xcut) go to 100  	  ! only do nearby points

do 100 j=1,ny			  ! loop on y gridpoints, same procedure
ygp=real(j)
dely=abs(ygp-yyg)
if (dely.gt.ycut) go to 100

xgas=(delx*dx/xsc)**2		  	! make gaussian exponents
ygas=(dely*dy/ysc)**2
expn=exp(-xgas-ygas)		! make the gaussian weight
wate(i,j)=wate(i,j)+expn	! sum the weights
grid(i,j)=grid(i,j)+val*expn	! sum the weighted values

100	continue

return
end

c************************************************************************

```