[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

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.