[Thread Prev][Thread Next][Index]

Re: ferret auto stats



Hi Cyndy and Joe,
	I started to answer this a few days ago but got side-tracked onto
something else.  Here is something that may still be timely and/or may
be of use to others in the group.

The hardest part of the problem is I think the automatic generation of a
histogram of a field that may be 1,2,3, or 4D, so at first put aside the
problem of 60 or more variables that you face. Using the script
autohist.jnl below the histogram of a chosen variable can be generated as
in the following demo

use coads_climatology
sho data
     currently SET data sets:
1> /opt/local/ferret/fer_dsets/descr/coads_climatology.des  (default)
 name     title                             I         J         K        L
SST      SEA SURFACE TEMPERATURE          1:180     1:90      ...       1:12
...
! pick SST and generate its histogram over all locations and months
let vhist=sst ; go autohist


After chugging along for a while (but not, I suspect, as long as some of
the other histogram tools which may involve a sort, or may not have 4D
capability) this produces a plot such as the attached GIF.

One could then pick another variable
   let vhist=sst ; go autohist
and repeat the process.

There may be a way I'm not seeing for automating a grab of all of the
variables present in a file, but a clunky way to make histograms that
would not be too labor intensive would be to cut and paste the results
of a "sho data" request to a file (say allvars.d) then do something
like
   awk '{print "let vhist="$1," ; go autohist"}' allvars.d > allvars.jnl
which in the case of the coads climatology would look like
let vhist=SST  ; go autohist
let vhist=AIRT  ; go autohist
let vhist=SPEH  ; go autohist
let vhist=WSPD  ; go autohist
let vhist=UWND  ; go autohist
let vhist=VWND  ; go autohist
let vhist=SLP  ; go autohist

Then a ferret session like
	use coads_climatology
	go allvars.jnl
would crunch away while you went bike-riding, lunch-eating, whatever ...

The autohist.jnl script is as follows, and its also attached.  Its not
fully ready for primetime - for example at present it assumes right now
that there is just one dataset loaded and is hardwired to bin the data
into about 20 classes.

The example shown is an unweighted histogram, with all the SST values
tossed into the hopper without regard to the area (or volume) each one
represents. It would not be too hard to modify autohist to produce a
weighted histogram I'd guess. Still, even as it stands, I think it may
be of use in the situation you described.

Good luck,
Mick

!-----------------------------------------------------------------------------
! autohist : automatically generates a histogram for a selected variable
!            over the current region in 4-d, without the need for sorting.
!
! usage : let vhist = selected variable
!         go autohist
!
def sym datafile `vhist,return=dset`
def sym histname `vhist,return=title`
!
let vhistn = vhist[i=@min,j=@min,k=@min,l=@min]
let vhistx = vhist[i=@max,j=@max,k=@max,l=@max]
let nvhpts = vhist[i=@ngd,j=@ngd,k=@ngd,l=@ngd]
!
! form ABOUT 20 evenly spaced histogram bins
!
ppl %range `vhistn` `vhistx` 20
!
! whose range and increment are given by the ppl symbol values
!     ($ppl$range_low),($ppl$range_high),($ppl$range_inc)
! the number of bins is given by
!
let nhistb = int((($ppl$range_high)-($ppl$range_low))/($ppl$range_inc))
!
! and an integer bin variable is
!
let binvar = int((vhist-($ppl$range_low))/($ppl$range_inc))+1
let binfreq = binsel[i=@ngd,j=@ngd,k=@ngd,l=@ngd]/nvhpts
let bin0 = ($ppl$range_low)-($ppl$range_inc)/2 ; let bind = ($ppl$range_inc)
! prepare polygon vertices
let boxx=ysequence({-0.5,0.5,0.5,-0.5}) ; let boxy=ysequence({0,0,1,1})
!
! write the bin frequencies to a temporary file ...
!
sp \/usr/bin/rm -f dohist.jnl ; sp touch dohist.jnl
sp \/usr/bin/rm -f histfile.tmp ; sp touch histfile.tmp
sp \/usr/bin/rm -f doit.tmpl
sp echo 'let binsel=if(binvar eq INDEX)then 1' > doit.tmpl
sp echo 'let bincent=bind*INDEX+bin0' >> doit.tmpl
sp echo 'list/nohead/app/form=(2e13.5)/file=histfile.tmp bincent,binfreq' >> doit.tmpl
repeat/i=1:`nhistb` (sp sed "s/INDEX/`i`/" doit.tmpl >> dohist.jnl)
go dohist
!
! ... then read them back ...
!
file/form=free/var=xhist,yhist histfile.tmp
!
! ... and plot them as a black polygon
!
plot/vs/sym=5/nolab/set xhist,yhist
ppl xaxis,($ppl$range_low),($ppl$range_high),($ppl$range_inc) ; ppl plot
polygon/nolab/pal=black/nokey/o boxx*bind+xhist,boxy*yhist
can data histfile.tmp
label/nouser 0 `($ppl$ylen)+0.3` -1 0 0.15 Histogram for ($histname)
label/nouser `($ppl$xlen)` `($ppl$ylen)+0.3` 1 0 0.15 Dataset: ($datafile)
label/nouser -0.7 `($ppl$ylen)/2` 0 90 0.15 Frequency
! show data extremes along the bottom axis
label/nouser 0 -0.5 -1 0 0.15 @sr`vhistn`
label/nouser `($ppl$xlen)` -0.5 1 0 0.15 @sr`vhistx`
!
! make a gif of the result
!
frame/file=($datafile)($histname).gif



|____Mick.Spillane@noaa.gov____|
|__Room 2070 Bldg#3 NOAA/PMEL__|
|____Phone_:_(206)526-6780_____|


Attachment: coads_climatologySST.gif
Description: coads_climatologySST.gif

! autohist : automatically generates a histogram for a selected variable
!            over the current region in 4-d, without the need for sorting.  
!
! usage : let vhist = selected variable
!         go autohist
!
def sym datafile `vhist,return=dset`
def sym histname `vhist,return=title`
!
let vhistn = vhist[i=@min,j=@min,k=@min,l=@min]
let vhistx = vhist[i=@max,j=@max,k=@max,l=@max]
let nvhpts = vhist[i=@ngd,j=@ngd,k=@ngd,l=@ngd]
!
! form ABOUT 20 evenly spaced histogram bins
!
ppl %range `vhistn` `vhistx` 20
!
! whose range and increment are given by the ppl symbol values
!     ($ppl$range_low),($ppl$range_high),($ppl$range_inc)
! the number of bins is given by
!
let nhistb = int((($ppl$range_high)-($ppl$range_low))/($ppl$range_inc))
!
! and an integer bin variable is 
!
let binvar = int((vhist-($ppl$range_low))/($ppl$range_inc))+1
let binfreq = binsel[i=@ngd,j=@ngd,k=@ngd,l=@ngd]/nvhpts
let bin0 = ($ppl$range_low)-($ppl$range_inc)/2 ; let bind = ($ppl$range_inc)
! prepare polygon vertices
let boxx=ysequence({-0.5,0.5,0.5,-0.5}) ; let boxy=ysequence({0,0,1,1})
! 
! write the bin frequencies to a temporary file ...
!
sp \/usr/bin/rm -f dohist.jnl ; sp touch dohist.jnl
sp \/usr/bin/rm -f histfile.tmp ; sp touch histfile.tmp
sp \/usr/bin/rm -f doit.tmpl
sp echo 'let binsel=if(binvar eq INDEX)then 1' > doit.tmpl
sp echo 'let bincent=bind*INDEX+bin0' >> doit.tmpl
sp echo 'list/nohead/app/form=(2e13.5)/file=histfile.tmp bincent,binfreq' >> doit.tmpl
repeat/i=1:`nhistb` (sp sed "s/INDEX/`i`/" doit.tmpl >> dohist.jnl)
go dohist
!
! ... then read them back ...
!
file/form=free/var=xhist,yhist histfile.tmp
!
! ... and plot them as a black polygon
!
plot/vs/sym=5/nolab/set xhist,yhist
ppl xaxis,($ppl$range_low),($ppl$range_high),($ppl$range_inc) ; ppl plot
polygon/nolab/pal=black/nokey/o boxx*bind+xhist,boxy*yhist
can data histfile.tmp
label/nouser 0 `($ppl$ylen)+0.3` -1 0 0.15 Histogram for ($histname)
label/nouser `($ppl$xlen)` `($ppl$ylen)+0.3` 1 0 0.15 Dataset: ($datafile)
label/nouser -0.7 `($ppl$ylen)/2` 0 90 0.15 Frequency
! show data extremes along the bottom axis
label/nouser 0 -0.5 -1 0 0.15 @sr`vhistn`
label/nouser `($ppl$xlen)` -0.5 1 0 0.15 @sr`vhistx`
!
! make a gif of the result
!
frame/file=($datafile)($histname).gif



[Thread Prev][Thread Next][Index]

Dept of Commerce / NOAA / OAR / PMEL / TMAP

Contact Us | Privacy Policy | Disclaimer | Accessibility Statement