[visad] construct a 3D grid from multiple 2D grids

• To: "visad@xxxxxxxxxxxxxxxx" <visad@xxxxxxxxxxxxxxxx>
• Subject: [visad] construct a 3D grid from multiple 2D grids
• From: Valentijn Venus <venus@xxxxxx>
• Date: Tue, 30 Nov 2010 21:35:14 +0100

Hi,

In Unidata's GridUtil package there is a method to slice a 3D grid to
obtain a 2D fieldimpl.

Now, I basically want the reverse: construct a 3D grid from many 2D grids
with similar coordinate systems.

Attached 20 discrete (MODIS) vertical levels (nlevs=20) corresponding 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

For each level there is the retrieved temperature field and a geopotential
height field.

So the idea is to take each temperature field and resample the
geopotential height field at each x,y,z point to end up with a field
temperature at pressure point (assuming a standard atmosphere).
This can be done, but the solution is not straight forward. As I have
limited understanding of the VisAD data model, could you help me in the
right direction?

Attached an IDV-bundle containing the input 2D grid, and below this
message a _very_ ugly attempt to get this going.

Cheers, Tyn


def makeMODISSwathTo3DSet(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())
  timeImage = FieldImpl(newType, newTime)
  for i in range(numTimes):
    timeImage.setSample(i, image[i],0)
  return timeImage




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