Re: Locate good data point

[Jonathan Callahan <callahan@xxxxxxxxxxxxx>]

Jens-Olaf,

I've attached a script to do what you want.  It uses some functions
which are moderately recent additions so check to see if they exist for
you:

yes? show func compressi
COMPRESSI(DAT)
    Returns data, compressed along the I axis: Missing points moved to
the end
    DAT: variable to compress in I
yes? show func samplei
SAMPLEI(TO_BE_SAMPLED,X_INDICES)
    sample a field at a list of X indices
    TO_BE_SAMPLED: data to sample at list of X indices supplied
    X_INDICES: list of X indices at which to sample


Please note that this script will only work for datasets which contain a
single ocean basin.  You cannot take a global dataset, SET REGION, and
expect it to work.  As stated in the beginning of the file, I had to
create a subset of levitus_climatology in order to write this script.

Good luck,


-- Jonathan Callahan

Jens-Olaf Beismann wrote:
> 
> Dear Ferreters,
> 
> using output from an Atlantic ocean model, I'm trying to locate the
> i-index of the westernmost ocean point at a given latitude. Example: At,
> say, 30 degN, the westernmost ocean point would be roughly at 80 degW at
> the surface, and a bit more to the east at greater depths.
> I could of course LIST/y=30n/k=1/form=(e) TEMP to see at which index the
> values switch from special_value to real data at this level, but since I
> want to do it for all depths and all latitudes, I'm looking for something
> more automatic.
> 
> Any good ideas?
> 
> Cheers,
> 
> Jens-Olaf
> _____________________________________________________
>  Jens-Olaf Beismann
>  IfM Kiel - FB 1
>  Duesternbrooker Weg 20        Phone:__49-431-6004016
>  24105 Kiel, Germany           Fax  :__49-431-6001515
>
! Script to calculate the western boundary of an ocean model.
!
! The dataset used is derived from levitus_climatology with
!
! yes? use levitus_climatology
! yes? save/file=bop.nc temp[x=100w:20e,y=10n:60n]
!
! So it is only the North Atlantic and consists of 20 levels.

use bop.nc

let ix = if temp then i[gx=temp]    ! 'I' indices of good temp values.
let c_ix = compressi(ix)            ! Push all bad flags to end.
let coast_index = c_ix[i=1]         ! First value is westernmost index.

stat coast_index

message/quiet "Note that coast_index is on a YZ grid"

let mylon = x[gx=temp]              ! Longitude values from grid.
let mylat = y[gy=temp]              ! Latitude values from grid.

! Now we'll sample the longitude values
! based on the coast_index values.
!
! NOTE:  We need to create one of these for each 'K'
! because the index values passed to SAMPLEI need to be 1D.

let coast1 = samplei(mylon,coast_index[k=1])
let coast12 = samplei(mylon,coast_index[k=12])

! A couple of demonstration plots to show that it works

set view left
shade/lev=(0,30,1) temp[x=100w:40w,k=1]
plot/vs/over/symbol coast1,mylat

set view right
shade/lev=(0,30,1) temp[x=100w:40w,k=12]
plot/vs/over/symbol coast12,mylat