3. Fields

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

Components of a meteorological field

A meteorological field has (mainly):

  • a fid, that identifies the field nature
  • a geometry, that defines the position of the of the field’s data; detailed in Geometries.
  • a validity, that defines the temporal validity of the field’s data
  • a spectral geometry, if the field is spectral.
>>> 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 is a dict() which keys are formats names (with a distinction between GRIB1 & GRIB2), e.g.

>>> field.fid
{'FA':'S058WIND.U.PHYS',
 'GRIB':{'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 must have a fid of the resource format !

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 (~ Hovmöller fields).

>>> 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

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
>>> field.sp2gp()

3.3.2. Data

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

  • basic statistics:

    >>> 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 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

The Field.plotfield() method may be one of the most useful methods, with a number of options continually growing... Although referring to each kind of field actual documentation is highly recommended, here is a short introduction:

>>> t = r.readfield('CLSTEMPERATURE')  # 2m temperature
>>> fig = t.plotfield()

then either

>>> fig.show()
>>> fig.savefig('my_figure.png')

Tip: as the most time-consuming part of creating the plot is the creation of the underlying Basemap object (mpl_toolkits.basemap), when creating several plots one should save and reuse the basemap object (cf. Geometries). This will save much time.

>>> temp = r.readfields('S00*TEMPERATURE')  # read highest levels temperature
>>> bm = temp[0].geometry.make_basemap()  # specific basemap can be addressed through optional arguments
>>> for t in temp:
...     t.sp2gp()
...     fig = t.plotfield(existingbasemap=bm, ...)
...     fig.savefig(t.fid['FA'] + '.png')

Similarly, superposition of plots can be done:

>>> t = r.readfield('CLSTEMPERATURE')
>>> geop = r.readfield('SPECSURFGEOPOTEN')
>>> geop.sp2gp()
>>> bm = t.geometry.make_basemap()
>>> fig = geop.plotfield(existingbasemap=bm, graphicmode='contourlines')
>>> fig = t.plotfield(existingbasemap=bm, existingfigure=fig)