3. Fields

>>> import epygram
>>> epygram.init_env()
>>> r = epygram.open('ICMSHAROM+0042', 'r')
>>> field = r.readfield('S058WIND.U.PHYS')

Components of a meteorological field

A meteorological field has (mainly):

>>> field.structure
'H2D'
>>> type(field.geometry)
<class 'epygram.geometries.H2DGeometry.H2DProjectedGeometry'>
>>> type(field.validity)
<class 'epygram.base.FieldValidityList'>

In the following, we focus on H2DField, but all kind of fields behave similarly.

3.1. Field identifier (fid)

The fid of a field handle its nature, i.e. the physical parameter and the material (surface, soil, atmosphere) it represents (often with a detailed vertical location). It is handled as a dict() which keys are formats names, because fields in files are stored under different fid depending on the format, e.g.

>>> field.fid
{'FA':'S058WIND.U.PHYS',
 'GRIB1':{'indicatorOfParameter':33,
          'indicatorOfTypeOfLevel':109,
          'level':58,
          'table2Version':1,}
 'GRIB2':{'discipline':0,
          'parameterCategory':3,
          'parameterNumber':2,
          'typeofFirstFixedSurface':119,
          'level':58,}
}

Take care that for a field to be written in a resource, the only-but-mandatory required field characteristics is that its fid in the resource format must be present!

In other words, a field whose fid is {'FA':'SURFTEMPERATURE', 'GRIB':{...}} will be writeable either in a GRIB or a FA file, but not in a LFI or any other format…

3.2. Validity

The validity attribute of a field is a list of epygram.base.FieldValidity, because a field can have a temporal dimension.

>>> field.validity[0].getbasis()
datetime.datetime(2014, 12, 1, 0, 0)
>>> field.validity[0].getbasis(fmt='IntStr')
20141201000000
>>> field.validity[0].term()
datetime.timedelta(1, 64800)

Some fields may have a cumulative duration, e.g. precipitation fields are accumulated:

>>> field = r.readfield('SURFACCGRAUPEL')
>>> field.validity[0].statistical_process_on_duration()
'accumulation'
>>> field.validity[0].cumulativeduration()
datetime.timedelta(0, 10800)

The statistical_process_on_duration, among which one can find min, max, average and so on, is coded as GRIB2 norm (http://apps.ecmwf.int/codes/grib/format/grib2/ctables/4/10), if available.

3.3. Fields useful methods

Note

Autocompletion, in interactive (i)Python session or smart editors, may be an even better (than doc/tuto) way to explore the available methods of objects.

3.3.1. Spectralness

Spectral transforms are done (“in place”) through the two methods Field.sp2gp() and Field.gp2sp(a_sp_geom):

>>> field = r.readfield('S090TEMPERATURE')
>>> field.spectral
True

The epygram.geometries.SpectralGeometry contains the kind of spectral space (bi-Fourier//LAM or Legendre//global), the truncation(s), and the actual spectral transforms routines.

>>> spgeom = field.spectral_geometry
>>> field.sp2gp()
>>> field.spectral
False
>>> field.spectral_geometry
None
>>> field.gp2sp(spgeom)
>>> field.spectral
True

3.3.2. Data

Some methods have been implemented to ease a comprehensive access to the data:

  • basic statistics:

    >>> field.sp2gp()
>>> field.stats()
{'std': 6.1556479204416092, 'nonzero': 2211840, 'quadmean': 280.99889135008499, 'min': 259.33480158698774, 'max': 293.21439026360537, 'mean': 280.93145950331785}

  • field value at some lon/lat point:

    >>> field.getvalue_ll(1.5, 45.6) # default interpolation = 'nearest'
274.278588481978
>>> field.getvalue_ll(1.5, 45.6, neighborinfo=True) # get info about the nearest neighbor gridpoint used
(274.278588481978, (1.4988145157970028, 45.60001936720281))
>>> field.getvalue_ll(1.5, 45.6, interpolation='linear')
274.25775674717687

Also, although the field’s data is accessible through its attribute data, it is strongly advised to access the data through the method Field.getdata(), because the internal storage of the data may differ from expected by the user.

  • modifying the field data should resemble:

    >>> data = field.getdata()
>>> type(data)
<type 'numpy.ndarray'>
>>> data.shape
(1536, 1440)
>>> data[100:800,500:600] += 10*numpy.random.rand(700,100)
>>> field.setdata(data)

    after what the field can of course be re-written in a resource.

  • some patterned operations on fields are facilitated through the Field.operation() method: any of the four basic operations (+,*,-,/) with scalars or any numpy function (exp, sin, log…):

    >>> field.operation('-', 273.15)  # e.g. go from K to °C
>>> field.operation('sin')  # does field.data = numpy.sin(field.data)

  • of spectral fields can also be computed horizontal derivatives:

    >>> t = r.readfield('S045TEMPERATURE')
>>> t.spectral
True
>>> (dx, dy) = t.compute_xy_spderivatives()
>>> type(dx)
<class 'epygram.fields.H2DField.H2DField'>
>>> dx.spectral
False
>>> dx.max()
0.0051387105385038408

  • of 2D fields can be computed spectra (epygram.spectra.Spectrum):

    >>> t.sp2gp()
>>> s = t.dctspectrum()
>>> type(s)
<class 'epygram.spectra.Spectrum'>

3.3.3. Operations between fields

Operations between fields can be done in two ways:

  • standard Python syntax; in case a new Field object is created, with uninitialized validity (what is the validity of an operation between two fields of potential different validity ?) and fid:

    >>> field90 = r.readfield('S090TEMPERATURE')
>>> field89 = r.readfield('S089TEMPERATURE')
>>> field_diff = field90 - field89

  • the Field.operation() method; in case the field values are modified “in place”:

    >>> field90 = r.readfield('S090WIND.U.PHYS')
>>> field89 = r.readfield('S089WIND.U.PHYS')
>>> field90.operation('+', field89)

In any case, a simple consistency check is done on the fields’ geometry, basically on their dimensions.

3.4. Building Vector Fields

Wind fields (for instance) can be re-assembled from their U/V components into H2DVectorField or D3VectorField for more integrated functionalities (re-projection, computation of derivatives or direction/module, plotting and so on…).

>>> u = r.readfield('S090WIND.U.PHYS')
>>> v = r.readfield('S090WIND.V.PHYS')
>>> wind = epygram.fields.make_vector_field(u,v)
>>> wind.sp2gp()

  • reprojection: FA wind fields are projected on the grid axes (here, a Lambert projection); let’s get the wind components on true zonal/meridian axes:

    >>> wind.getvalue_ij(0,0)
[0.5525041298918116, -2.8212975453933336]
>>> wind.reproject_wind_on_lonlat()
>>> wind.getvalue_ij(0,0)
[0.9307448483516318, -2.759376908801778]

  • derivatives: just as the Field.compute_xy_spderivatives() method enable to compute derivatives of spectral fields, the H2DVectorField.compute_vordiv() method enable to compute vorticity and divergence of a spectral wind field:

    >>> wind.gp2sp(r.spectral_geometry)
>>> (vor, div) = wind.compute_vordiv()
>>> type(vor)
<class 'epygram.fields.H2DField.H2DField'>

  • direction/module: to compute a wind direction or wind module field from vectors:

    >>> wind.sp2gp()
>>> ff = wind.to_module()
>>> type(ff)
<class 'epygram.fields.H2DField.H2DField'>

3.5. Plots

Cf. Plots sections of the Gallery.

3.6. 3D Plots

Fields can be plotted in a 3D view in three ways:

  • contour with Field.plot3DContour() method (which becomes a classical contour plot if field is 2D, vertical or horizontal)

  • volume with Field.plot3DVolume()

  • color with Field.plot3DColor() (corresponding to the matplotlib contourf method)

A vector field can be plotted using Field.plot3DVector() (to plot arrows) or Field.plot3DStream() to plot (stream lines or tubes).

>>> import vtk #We need to import vtk before epygram even if do not use it directly in the script
>>> import epygram
>>> epygram.init_env() #initialisation of environment, for FA/LFI and spectrals transforms sub-libraries
>>> r = epygram.open(filename, 'r', true3d=True)
>>>
>>> CF = r.readfield('S---CLOUD_FRACTI')
>>>
>>> #Set-up of the view
... offset = CF.geometry.gimme_corners_ll()['ll'] #We translate the domain
>>> hCoord = 'll' #We use lat/lon on the horizontal
>>> z_factor = 0.1 #0.1 horizontal degree of lat/lon is represented by the same length as one model level on the vertical
>>> ren = epygram.util.vtk_set_window((0.5, 0.5, 0.5), (800, 800))
>>>
>>> CF.plot3DContour(ren, [1.], color='White', hCoord=hCoord, offset=offset, z_factor=z_factor)
((vtkRenderingOpenGL2Python.vtkOpenGLActor)0x7f899c230390, (vtkRenderingOpenGL2Python.vtkOpenGLPolyDataMapper)0x7f89b8750a78)
>>>
>>> ren['interactor'].Start()