Fields ====== .. highlight:: python >>> 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 :ref:`fid `, that identifies the field nature - a :ref:`geometry `, that defines the position of the of the field's data; detailed in :doc:`../library/geometries_lib` - a :ref:`validity `, that defines the temporal validity of the field's data - a :ref:`spectral geometry `, if the field is spectral. >>> field.structure 'H2D' >>> type(field.geometry) >>> type(field.validity) In the following, we focus on H2DField, but all kind of fields behave similarly. ----------------------------------------------------------- .. _tuto-fid: 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... ----------------------------------------------------------- .. _tuto-validity: Validity -------- The ``validity`` attribute of a field is a list of :class:`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. ----------------------------------------------------------- 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.* .. _tuto-spectral: Spectralness ^^^^^^^^^^^^ Spectral transforms are done ("in place") through the two methods :meth:`Field.sp2gp() ` and :meth:`Field.gp2sp(a_sp_geom) `: >>> field = r.readfield('S090TEMPERATURE') >>> field.spectral True The :class:`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 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 :meth:`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) >>> 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 :meth:`Field.operation() ` method: any of the four basic operations (+,*,-,/) with scalars or any :mod:`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) >>> dx.spectral False >>> dx.max() 0.0051387105385038408 - of 2D fields can be computed spectra (:class:`epygram.spectra.Spectrum`): >>> t.sp2gp() >>> s = t.dctspectrum() >>> type(s) 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 :meth:`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. ----------------------------------------------------------- Building Vector Fields ---------------------- Wind fields (for instance) can be re-assembled from their U/V components into :doc:`H2DVectorField <../library/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 :meth:`Field.compute_xy_spderivatives` method enable to compute derivatives of spectral fields, the :meth:`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) - direction/module: to compute a wind direction or wind module field from vectors: >>> wind.sp2gp() >>> ff = wind.to_module() >>> type(ff) ----------------------------------------------------------- Plots ----- The :meth:`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, ax = t.plotfield() then either >>> fig.show() >>> fig.savefig('my_figure.png') .. note:: as the most time-consuming part of creating the plot is the creation of the underlying Basemap object (:mod:`mpl_toolkits.basemap`), when creating several plots one should save and reuse the basemap object (cf. :doc:`../library/geometries_lib`). 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, ax = t.plotfield(use_basemap=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, ax = geop.plotfield(use_basemap=bm, graphicmode='contourlines') >>> fig, ax = t.plotfield(use_basemap=bm, over=(fig, ax))