NOTICE: This version of the NSF Unidata web site (archive.unidata.ucar.edu) is no longer being updated.
Current content can be found at unidata.ucar.edu.
To learn about what's going on, see About the Archive Site.
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_28This 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.