Re: [idvusers] assign vertical coordinate (value) to flat ADDE image

• To: Don Murray <Don.Murray@xxxxxxxx>, "idvusers@xxxxxxxxxxxxxxxx" <idvusers@xxxxxxxxxxxxxxxx>
• Subject: Re: [idvusers] assign vertical coordinate (value) to flat ADDE image
• From: Valentijn Venus <venus@xxxxxx>
• Date: Wed, 16 Feb 2011 16:36:03 +0100

Hi Don,

> It's like the reverse of GridUtil.makeDomain2D, but harder. ;-)

Yes, I have seen a method in Unidata's GridUtil package that slices a 3D grid 
to obtain a 2D fieldimpl.... And yes, the reverse is indeed harder ;-(

The reverse: construct a 3D grid from a manifold of 2D singleBanded ADDE 
images, all sharing the same coordinate systems, yet with a missing z-level.

The 20 discrete (MODIS) vertical levels to the following pressure-levels (in 
mb):
 5,  10,  20,  30,  50,  70,  100 , 150 ,  200 , 250 ,  300 ,   400  ,
500 , 620 ,  700 ,   780 ,  850 ,  920 ,  950, 1000

Where the lowest index 0=5mb, ..., 20=1000mb.

For each level there are 3 fields: the retrieved temperature (K)¸ retrieved 
moisture (dewpoint) temperature, and a geopotential height (msl) field.

Below this message a _very_ ugly attempt but I clearly goes beyond my 
capacity.... Attached to this message I provide an IDV-bundle containing the 
input data.

Cheers, Tyn


def modisSwathTo3DSet(a,b):
  from visad import FunctionType, FieldImpl, Gridded3DDoubleSet, Set, Unit
  from visad import RealType
  from visad import RealTupleType
  from visad import Gridded3DSet
  from visad import CartesianProductCoordinateSystem

  from ucar.unidata.data.grid import GridUtil
  import ucar.visad.quantities.AirPressure as ap

  timedom = a.getDomainSet()
  numTimes = timedom.getLength()

  dom1 = GridUtil.getSpatialDomain(a)
  xyCs = dom1.getCoordinateSystem()

  samples = dom1.getSamples(0)
  zCs = ap.getStandardAtmosphereCS()
  #print zCs
  print zCs.getLengths()
  newCS = CartesianProductCoordinateSystem(xyCs, zCs)

  sdomain = Gridded3DSet(RealTupleType.SpatialEarth3DTuple, samples,
                        dom1.getX().getLength(), dom1.getY().getLength(),
                        None,
                        None,
                        dom1.getSetErrors(), 0, 0)

  image = GridUtil.setSpatialDomain(a,sdomain)
  newTime = Gridded3DDoubleSet(RealType.Time, timedom.getSamples(0),
                               numTimes,
                               timedom.getCoordinateSystem(),
                               timedom.getSetUnits(),
                               timedom.getSetErrors(), 0)
  newType = FunctionType(RealTupleType.Time1DTuple, image[0].getType())


Faculty of Geo-Information Science and Earth Observation (ITC)
University of Twente
Chamber of Commerce: 501305360000

E-mail disclaimer
The information in this e-mail, including any attachments, is intended for the 
addressee only. If you are not the intended recipient, you are hereby notified 
that any disclosure, copying, distribution or action in relation to the content 
of this information is strictly prohibited. If you have received this e-mail by 
mistake, please delete the message and any attachment and inform the sender by 
return e-mail. ITC accepts no liability for any error or omission in the 
message content or for damage of any kind that may arise as a result of e-mail 
transmission.

Attachment: MODIS2Dto3D.xidv
Description: MODIS2Dto3D.xidv