20040915: Geopotential Isosurface Plot

Hi BIll-

Sorry for the delay in responding to this.

> I am viewing model data, and am trying to display "geopotential height of the tropopause" using a 
> 3-D "topography" like display. When I use the GFS grid "Geopotential height at surface of the 
> earth" I get the expected result. So, I tried a parallel approach.
>
> I have written a jython method to calculate the height of the trop. from the GFS grid 
> "pressure at tropopause".
>
> def heightTrop(p):
>    """ height of pressures from U.S. Standard Atmosphere """
>    # calculate height for a constant lapse rate (6.5 C/km) atmosphere
>    z=(288.15/.0065)*(1-(p/101325.)**(287.05*.0065/9.806))
>    # change height in stratosphere to using isothermal (216.65 K)
>    for i in range(len(z)):
>      if p[i] <  22634.: 
>        z[i]=11000+216.65*287.05*(22634.-p[i])/(p[i]*9.806)
>    return z
>
> I have created a formula:
> name - z_sfc
> formula - heightTrop(pressureTrop)
> displays include - plan views, probes, and topography.
>
> I have displayed contours, and they look fine.
>
> After I select the topography display type, create the display, and choose the 
> GFS grid, I get an error message:
> An error has occurred: ControlDescriptor.Creating display
> Unable to handle units of Generic_1_nullUnit_28
>
> This lools like a problem with units ??? 

Yes.  When you do all the multiplication, you are just using
numbers and ignoring units.  In the process, you end up with
a unit that is not covertible with a unit of height.  Thus,
the error message.

> Am I on the right track ???

What you really want to do is take the pressure field and
resample the geopotential height field at each x,y,p point
to end up with a field heights at each pressure point.
This can be done, but the solution is not straight forward.
Basically, you need to take the pressure field at each time
step, convert it to a sampling set, resample the corresponding
geopotential height field at that timestep and display it
as a topopgraphy field.  The following code can be used to
do that:

def makeSampleFrom2DPressureField(field):
   import ucar.unidata.data.grid.GridUtil as gu
   import ucar.visad.quantities.AirPressure as ap
   from visad import *
   if gu.isTimeSequence(field):
     field = field.getSample(0)
   domain = gu.getSpatialDomain(field)
   samples = domain.getSamples()
   lengths = domain.getLengths()
   dtype = domain.getType().getDomain()
   rtypes = dtype.getRealComponents()
   dunits = domain.getSetUnits()
   dcs = domain.getCoordinateSystem()
   rsamples = field.getFloats()
   rangetype = rangeType(field)
   if isinstance(rangetype, RealTupleType):
      rangetype = rangetype.getComponent(0)
   print "samples = ", len(samples)
   newD = [samples[0], samples[1], rsamples[0]]
   newU = [dunits[0], dunits[1], rangetype.getDefaultUnit()]
   newCS = CartesianProductCoordinateSystem(dcs, 
ap.getStandardAtmosphereCS())
   newDT = RealTupleType(rtypes[0], rtypes[1], rangetype, newCS, None)
   newSet = Gridded3DSet(newDT, newD, lengths[0], lengths[1],
                         None, newU, None)
   return newSet

def makeTropopauseSurface(pressures, heights):
   import ucar.unidata.data.grid.GridUtil as gu
   from visad import *
   timeset = gu.getTimeSet(pressures)
   newField = None
   for i in range(timeset.length):
       pressureSet = makeSampleFrom2DPressureField(pressures[i])
       resampledHeight = gu.resampleGrid(heights[i], pressureSet)
       resampledHeight = gu.make2DGridFromSlice(resampledHeight)
       if i == 0:
          ftype = FunctionType(domainType(heights), 
resampledHeight.getType())
          newField = FieldImpl(ftype, timeset)
       newField.setSample(i, resampledHeight, 0)
   return newField

Note this requires an understanding of the VisAD data model.
It also assumes that the pressure and height fields are
on the same projection.

Now, create a Formula that would be:

name: troposurface
formula: makeTropopauseSurface(pressures, heights)

For displays, you can select topography.

Execute the formula and select the pressure at the tropopause for the 
pressure field and the geopotential height at isobaric levels for the 
height field.  Display it as topography.  The gaps are where the 
tropopause is above 100 hPa.  If you want to display contours in
3D, you could comment out the line:

       resampledHeight = gu.make2DGridFromSlice(resampledHeight)

See if this gives you the desired results and let me know what
other twists you'd like on this.