comp.lang.idl-pvwave archive
Messages from Usenet group comp.lang.idl-pvwave, compiled by Paulo Penteado

Home » Public Forums » archive » Re: Gridding options
Show: Today's Messages :: Show Polls :: Message Navigator
E-mail to friend 
Switch to threaded view of this topic Create a new topic Submit Reply
Re: Gridding options [message #21478] Tue, 29 August 2000 00:00
Craig Markwardt is currently offline  Craig Markwardt
Messages: 1869
Registered: November 1996
Senior Member
Ben Tupper <btupper@bigelow.org> writes:

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

Okay now I understand. So in this case X would be the distance along
the cruise path, and Y would be the depth from the surface.

...
> 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?

The more appropriate question is probably, how broad should the
gaussian be in X and Y? This depends on how much smoothing you want
to acomplish, and the new sampling. For example, if your original
sampling was 10-20 km, then the interpolated image might have ~2 km
resolution. With minimal smoothing, the gaussian sigma would be
around 15 km (ie, comparable to your sampling). The response function
should have around +/- 2 sigmas = +/- 30 km, which is about 30 pixels.

Craig

--
------------------------------------------------------------ --------------
Craig B. Markwardt, Ph.D. EMAIL: craigmnet@cow.physics.wisc.edu
Astrophysics, IDL, Finance, Derivatives | Remove "net" for better response
------------------------------------------------------------ --------------
Re: Gridding options [message #21482 is a reply to message #21478] Tue, 29 August 2000 00:00 Go to previous message
Ben Tupper is currently offline  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
Re: Gridding options [message #21484 is a reply to message #21478] Tue, 29 August 2000 00:00 Go to previous message
Craig Markwardt is currently offline  Craig Markwardt
Messages: 1869
Registered: November 1996
Senior Member
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
------------------------------------------------------------ --------------
  Switch to threaded view of this topic Create a new topic Submit Reply
Previous Topic: Re: Keyword precedence
Next Topic: New and shiny web pages

-=] Back to Top [=-
[ Syndicate this forum (XML) ] [ RSS ] [ PDF ]

Current Time: Wed Oct 08 13:36:23 PDT 2025

Total time taken to generate the page: 0.00524 seconds