| zonal means [message #11082] |
Mon, 16 March 1998 00:00  |
Martin Schultz
Messages: 515
Registered: August 1997
|
Senior Member |
|
|
Hi everyone,
this is part question, part answer - I just want to make sure there
is nothing wrong with this:
Q: How do you compute zonal means from a 3D data cube ?
(example: A(72,46,14) is a data array with longitude, latitude, altitude
as dimensions, and I want to compute the averages over longitude for
each latitude and altitude)
A: well, you can do it in a loop (buuuuuhh!!)
for j=0,13 do begin
for i=0,45 do begin
b(i,j) = total(a(*,i,j)) / 72.
endfor
endfor
Q: but I hate loops !!
A: hmmm, we could try REBIN or CONGRID. How about that:
c = reform(rebin(a,1,46,14),46,14)
d = reform(congrid(a,1,46,14),46,14)
Q: That looks much more like IDL! But does it work ?
A: (after running a test with random data) Seems like REBIN is
just what you are looking for. CONGRID somehow screws it up. BTW: I am
sure you want to know the timing of these routines:
loop: 1.3769450 seconds
rebin: 0.15734899 seconds
congrid: 0.085584044 seconds
These results are obtained by repeating each method 100 times.
Q: But can I be ascertained that REBIN does the job correctly ? And does
anyone have an idea why CONGRID screws up ?
Thanks,
Martin.
------------------------------------------------------------ -------
Dr. Martin Schultz
Department for Earth&Planetary Sciences, Harvard University
186 Pierce Hall, 29 Oxford St., Cambridge, MA-02138, USA
phone: (617)-496-8318
fax : (617)-495-4551
e-mail: mgs@io.harvard.edu
IDL-homepage: http://www-as.harvard.edu/people/staff/mgs/idl/
------------------------------------------------------------ -------
|
|
|
|
| Re: zonal means [message #11244 is a reply to message #11082] |
Wed, 18 March 1998 00:00  |
Craig Markwardt
Messages: 1869
Registered: November 1996
|
Senior Member |
|
|
Martin Schultz <mgs@io.harvard.edu> writes:
>
> Q: ... And does anyone have an idea why CONGRID screws up ?
First, CONGRID uses a different form of interpolation compared to
REBIN. REBIN will find the simple average of array cells that are
collapsed, while CONGRID can be set to perform an interpolation. It's
not clear to me that these are the same.
Also, CONGRID() has a subtle bug in its implementation, at least from
the point of view of a physicist :-) Essentially, CONGRID attempts to
interpolate to the pixel *corners*, not their centers. This becomes a
problem at the edges of the array; one "extra" row and column appears
at the far edges.
At some point the IDL folks tried to fix CONGRID by adding the
MINUS_ONE keyword, which removes those extra pieces. Unfortunately
that doesn't solve the problem, if you are trying to be precise, since
you lose one row and column. RSI support seems to have confirmed my
report of this problem.
The proper solution is to interpolate to the pixel centers. I have a
modified version of CONGRID on my IDL web page called CMCONGRID,
http://astrog.physics.wisc.edu/~craigm/idl/idl.html
which does just that. CMCONGRID is a drop in replacement for CONGRID.
Use it just like CONGRID, except that you also need to set the keyword
HALF_HALF to enable the pixel-center interpolation.
Craig
--
------------------------------------------------------------ --------------
Craig B. Markwardt, Ph.D. EMAIL: craigmnet@astrog.physics.wisc.edu
Astrophysics, IDL, Finance, Derivatives | Remove "net" for better response
------------------------------------------------------------ --------------
|
|
|
|
|
|
|
|
|
|
|
|
|
|
| Re: zonal means [message #11339 is a reply to message #11082] |
Mon, 23 March 1998 00:00 |
Martin Schultz
Messages: 515
Registered: August 1997
|
Senior Member |
|
|
James Tappin wrote:
>
>
> Now for the $64,000 question -- how do you do the same thing for
> medians? What we need here is an index argument to SORT analogous to the
> second argument of TOTAL. Methinks that the only efficient solution is
> currently to write the confounded thing as a CALL_EXTERNAL piece of C
> code.
>
while we are at it: a "dimension" parameter for MIN and MAX would also
be useful. But is it really worth 64k ?
Martin.
--
------------------------------------------------------------ -------
Dr. Martin Schultz
Department for Earth&Planetary Sciences, Harvard University
186 Pierce Hall, 29 Oxford St., Cambridge, MA-02138, USA
phone: (617)-496-8318
fax : (617)-495-4551
e-mail: mgs@io.harvard.edu
IDL-homepage: http://www-as.harvard.edu/people/staff/mgs/idl/
------------------------------------------------------------ -------
|
|
|
|
| Re: zonal means [message #11342 is a reply to message #11082] |
Mon, 23 March 1998 00:00 |
James Tappin
Messages: 54
Registered: December 1995
|
Member |
|
|
Evilio del Rio wrote:
>
> Martin Schultz wrote:
>> ...
>> Q: How do you compute zonal means from a 3D data cube ?
>> (example: A(72,46,14) is a data array with longitude, latitude, altitude
>> as dimensions, and I want to compute the averages over longitude for
>> each latitude and altitude)
> In my opinion you should use the TOTAL function with a second argument:
>
> IDL> help,a
> A FLOAT = Array[72, 46, 14]
> IDL> b = TOTAL(A,1) ; The argument 1 tells TOTAL to sum just in the
> first dim.
Now for the $64,000 question -- how do you do the same thing for
medians? What we need here is an index argument to SORT analogous to the
second argument of TOTAL. Methinks that the only efficient solution is
currently to write the confounded thing as a CALL_EXTERNAL piece of C
code.
--
+------------------------+-------------------------------+-- -------+
| James Tappin, | School of Physics & Astronomy | O__ |
| sjt@star.sr.bham.ac.uk | University of Birmingham | -- \/` |
| Ph: 0121-414-6462. Fax: 0121-414-3722 | |
+--------------------------------------------------------+-- -------+
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
