Re: Gridding options

[Craig Markwardt <craigmnet(at)cow.physics.wisc.edu>]

Ben Tupper <btupper@bigelow.org> writes:

> Hello,
> 
> I'm staring (again) at largish set of CTD casts from a recent cruise.
> The cast data is comprised of sample information from every 0.5 meters
> from the surface to the seafloor.   The 20 or so casts are separated
> from each other by about 10-20km and are nearly colinear.   I need to
> interpolate a 2d grid from these values.   In the past I have used the
> techniques described to grid the data.   I list them here in hopes that
> someone familiar with this kind of data can suggest alternatives.

I don't exactly understand what your data is like.  It sounds like you
have 0.5 m x 15000 m resolution, ie. extremely well sampled along one
axis and poorly sampled along another.  If that's the case, then the
following description may need to be modified.

I will describe my situation.  I have irregularly sampled data points,
which I wish to place on a regularly sampled 2D grid.  In my case the
resolution in X and Y is equal.  The measured data values are noisy,
so some form of averaging/smoothing is desireable.

My solution was to essentially convolve the measured points by a
spatial response function.  In my case it is the intrinsic spatial
response function of the measuring instrument, but a gaussian will
probably do fine for you.  Clearly you would want to tune the
parameters of your gaussian to be appropriate for your problem
(considering the spacing and noisiness of the data).  The trick is to
maintain the data and weighting functions separately, and divide them
at the end.  This provides a very natural weighting of nearby -- and
even overlapping -- data points.

Here is an example.  Suppose that your data is sampled at X and Y,
with value Z.  This example extends to more measurements trivially.
You are interested in making a MAP in the range [X0,X1] and [Y0,Y1],
in a NXBINS x NYBINS array.  The response function is RESP, an NRX x
NRY array: this is the gaussian, which should be centered at
RESP[NX/2,NY/2].  Here is my solution, with the real work being done
in the "drizzle" section.  Yes, a loop!

;; Discretize the positional values to IX And IY
xbinsize = (x1-x0)/nxbins 
ybinsize = (y1-y0)/nybins
ix = round((x-x0)/xbinsize) - nrx
iy = round((y-y0)/ybinsize) - nry

;; Make sure we keep all values in-bounds
wh = where(ix GE 0 AND ix LT nxbins-nrx AND iy GE 0 AND iy LT nybins-nry, ct)
if ct EQ 0 then $
  message, 'ERROR: no data within grid limits'
ix = ix(wh) & iy = iy(wh)
iz = z(wh)

;; Drizzle the points onto the map
map = dblarr(nxbins, nybins) & xmap = map & wmap = map
for i = 0L, ct-1 do begin
   map(ix(i),iy(i)) =  map(ix(i):ix(i)+nrx-1,iy(i):iy(i)+nry-1) + resp*iz(i)
  xmap(ix(i),iy(i)) = xmap(ix(i):ix(i)+nrx-1,iy(i):iy(i)+nry-1) + resp
endfor

;; Compute the weighted, convoluted map by dividing the data by the weighting
wh = where(xmap GT 0)
wmap(wh) = map(wh) / xmap(wh)

Maybe this helps!

Craig

-- 
--------------------------------------------------------------------------
Craig B. Markwardt, Ph.D.         EMAIL:    craigmnet@cow.physics.wisc.edu
Astrophysics, IDL, Finance, Derivatives | Remove "net" for better response
--------------------------------------------------------------------------