| Re: Gridding options [message #21482 is a reply to message #21478] |
Tue, 29 August 2000 00:00   |
Ben Tupper
Messages: 186
Registered: August 1999
|
Senior Member |
|
|
Craig Markwardt wrote:
>
> 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.
You have the right idea. The ship traveled along a long (mostly) straight
path. Every 10-20km the vessel stops and drops the CTD overboard, sampling
every 0.5 m over a total depth of 50m - 200m.
>
> 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
>
I do see what you are describing. This is quite similar (in methodology) to the
iterative gridding process used by a built in function GRID in PV-Wave (which I
am not using.)
How are NRX and NRY, for the response function, determined?
Thanks, Ben
--
Ben Tupper
Bigelow Laboratory for Ocean Science
West Boothbay Harbor, Maine
btupper@bigelow.org
note: email address new as of 25JULY2000
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
