20041004: horizontal section in ocean-s-coordinate?

>From: Rich Signell <address@hidden>
>Organization: USGS
>Keywords: 200410042018.i94KIMUE018549 IDV horizontal slice

Hi Rich-

>Now that my vertical section issue in ocean-s coordinate is solved, I've 
>been trying to figure out if there is a "horizontal section" display or 
>formula such  that I could plot, say, the temperature at 5 m depth (at 
>constant Z).  This is something we want to do all the time. 

you users - always wanting something. ;-)

>There doesn't seem to be a tool built in, but is this something I could 
>craft with jython on my own?   Or do I have to ask you guys and wait a 
>whole couple of hours to get it?  ;-)

Well, now you've opened quite a can of worms.  Ideally, it should
just be a matter of modifying the existing System->Slice at Level
formula (and corresponding Jython) to create a level which is in units 
of altitude.  The existing formula assumes that the input level is in 
the native coordinates (i.e. pressure, ocean_s, etc).  So, in a perfect world,
you could write a new Jython method:

import sys;
sys.add_package('visad');
sys.add_package('visad.data.units');
def getSliceAtAltitude(fieldimpl, alt, unit="m") :
  """ extract a 2D horizontal slice from a 3D grid at the given altitude
      level is a real number; if unit is supplied, it must
      be compatible with meters (ft, fathoms, etc)
      param fieldimpl is a VisAD FieldImpl which may have
      one or more time steps.  """
  #import methods from
  import ucar.unidata.data.grid.GridUtil as gu
  from visad import RealType
  from visad import Real
  import visad.data.units.Parser as parser
  alt = float(alt)
  unit = parser.parse(str(unit))
  altitude = Real(RealType.Altitude, alt, unit)
  ff = gu.sliceAtLevel(fieldimpl, altitude)
  return ff

and then create a formula to call it:

Name: altitudeSlice
Formula: getSliceAtAltitude (Field, user_Altitude, user_Unit[default=m])
(As of 1.1, could also be written as:
getSliceAtAltitude (Field, Altitude[isuser=true], Unit[isuser=true,default=m])
)
Description: Slice at Altitude
Group: My Cool Formulas
Displays: Plan Views

Then, select the Formulas Data Source and for the Fields, select
"Slice at Altitude".  You will be prompted for the user_/isuser fields,
then the field you want to slice.


As I said, this would be in a perfect world.  Unfortunately, this won't
work with the version I just gave you.  The problem is that the
GridUtil.sliceAtLevel method always assumed there was a 1:1 mapping
between a particular altitude and a particular native coordinate. 
Obviously, this is not the case with these ocean_s grids.   I took
a stab at refactoring the code to do a slice on the altitude grid,
but this has some far reaching implications in the bowels of the
IDV that could cause me to back out of it.  I put a new jar file out at:

ftp://ftp.unidata.ucar.edu/pub/idv/untested/idv.jar

which works with the code above and your sample grid.  I'm not sure
what it will break. ;-).  So, definitely back up the last JAR I
put there.

Please note that the input altitudes are positive up.  So, if you
want to use fathoms, make sure you put in that minus sign to indicate
that they are below sea level if that's what you want.  Alternatively,
you could just modify the jython to be sliceAtDepth and do the 
negation as:

alt = -float(alt)

for those who are used to positive down.

Let me know if this causes anything to break.  I'll have to really
test this out before I make this part of the core.  But in principal,
it should work.

Don Murray
NOTE: All email exchanges with Unidata User Support are recorded in the
Unidata inquiry tracking system and then made publicly available
through the web.  If you do not want to have your interactions made
available in this way, you must let us know in each email you send to us.

> Date:    Tue, 05 Oct 2004 09:55:59 EDT
> To:      Unidata Support <address@hidden>
> From:    Rich Signell <address@hidden>
> Subject: Re: 20041004: horizontal section in ocean-s-coordinate?
> 
> Fantastic.   Again!
> 
> I've played with new idv.jar supporting horizontal slicing, and I've 
> even added vertical sections, isosurfaces, and it doesn't break.   It 
> looks great from here!