Re: Vectorize a XYT box or, convert a 3D var to 1D

[Andrew Wittenberg <Andrew.Wittenberg@xxxxxxxx>]

On Fri, 17 Sep 2004, Paulo B. Oliveira wrote:
> I need to compute some simple statistics (median, quantiles) on a XYT
> box which depend on sorting data values. Currently, I'm saving the box
> values on an ascii file, read them on a single axis variable then, using
> sortk, I compute the values I need. However this process is very slow.
> Is there another way to create a 1D variable from a 3D XYT box ?

Hi Paulo,

The ZSEQUENCE function does what you need.  Here's an example showing how
to compute quartiles (and more generally, quantiles/percentiles) from some
2-dimensional data:

use coads_climatology
set mode interp

! unwrap January tropical Pacific SSTs into a 1-dimensional array
let a = zsequence(sst[l=1,y=30s:30n,x=120e:80w])

! define a quantile axis
let ngood = `a[k=@ngd]`
def ax/z=0:1/np=`ngood` z_q

! sort the data (missing values will be last)
let a_z = samplek(a,sortk(a))
 
! place the sorted valid data onto the quantile axis
let/title="quantiles of `sst,r=title`"/unit="`sst,r=unit`" a_q = a_z[gz=z_q@asn]

! plot the quantile function (inverse cumulative distribution function)
plot/trans a_q

! list the min, 1st quartile, median, 3rd quartile, and max 
list a_q[z=0],a_q[z=.25],a_q[z=.5],a_q[z=.75],a_q[z=1]

Apparently these SST data are skewed cold, i.e. the coldest temperatures
are more extreme (farther from the median) than the warmest temperatures.

Cheers,

Andrew

--
Dr. Andrew T. Wittenberg
Physical Scientist
GFDL/NOAA, Princeton, NJ