Re: Contouring "strangely" gridded data

[kessler@xxxxxxxxxxxxxxxxxxx (Billy Kessler)]

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

	realbad=1.e35

	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
		grid(i,j,m)=realbad
	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************************************************************************