#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
H2DField:
Contains the class that handle a Horizontal 2D field.
"""
import copy
import numpy
from epygram import config, util, epygramError
from epygram.base import Field, FieldValidity
from epygram.geometries import H2DGeometry, SpectralGeometry
[docs]class H2DField(Field):
"""
Horizontal 2-Dimensions field class.
A field is defined by its identifier 'fid',
its data, its geometry (gridpoint and optionally spectral),
and its validity.
The natural being of a field is gridpoint, so that:
a field always has a gridpoint geometry, but it has a spectral geometry only
in case it is spectral.
"""
_collector = ('field',)
_footprint = dict(
attr=dict(
geometry=dict(type=H2DGeometry),
validity=dict(
type=FieldValidity,
optional=True,
default=FieldValidity()),
spectral_geometry=dict(
type=SpectralGeometry,
optional=True),
processtype=dict(
optional=True,
info="Generating process.")
)
)
@property
[docs] def spectral(self):
"""
Returns True if the field is spectral.
"""
return self.spectral_geometry != None
##############
# ABOUT DATA #
##############
[docs] def sp2gp(self):
"""
Transforms the spectral field into gridpoint, according to its spectral
geometry. Replaces data in place.
The spectral transform subroutine is actually included in the spectral
geometry's *sp2gp()* method.
"""
if self.spectral:
if self.geometry.rectangular_grid:
# LAM
gpdims = {}
for dim in ['X', 'Y', 'X_CIzone', 'Y_CIzone']:
gpdims[dim] = self.geometry.dimensions[dim]
else:
# global
gpdims = {}
for dim in ['lat_number', 'lon_number_by_lat']:
gpdims[dim] = self.geometry.dimensions[dim]
data2d = self.geometry.reshape_data(self.spectral_geometry.sp2gp(self.data, gpdims))
self._attributes['spectral_geometry'] = None
self.setdata(data2d)
[docs] def gp2sp(self, spectral_geometry):
"""
Transforms the gridpoint field into spectral space, according to the
*spectral geometry* mandatorily passed as argument. Replaces data in
place.
The spectral transform subroutine is actually included in the spectral
geometry's *gp2sp()* method.
"""
if not self.spectral:
if not isinstance(spectral_geometry, SpectralGeometry):
raise epygramError("a spectral geometry (SpectralGeometry\
instance) must be passed as argument to\
gp2sp()")
if self.geometry.rectangular_grid:
# LAM
if self.geometry.grid['LAMzone'] != 'CIE':
raise epygramError("this field is not bi-periodicized: it\
cannot be transformed into spectral\
space.")
gpdims = {}
for dim in ['X', 'Y', 'X_CIzone', 'Y_CIzone']:
gpdims[dim] = self.geometry.dimensions[dim]
data1d = self.data.flatten()
else:
# global
gpdims = {}
for dim in ['lat_number', 'lon_number_by_lat']:
gpdims[dim] = self.geometry.dimensions[dim]
data1d = self.data.compressed()
self._attributes['spectral_geometry'] = spectral_geometry
self.setdata(spectral_geometry.gp2sp(data1d, gpdims))
[docs] def getdata(self, subzone=None):
"""
Returns the field data, with 2D shape if the field is not spectral,
1D if spectral.
- *subzone*: optional, among ('C', 'CI'), for LAM fields only, returns
the data resp. on the C or C+I zone.
Default is no subzone, i.e. the whole field.
Shape of 2D data: \n
- Rectangular grids:\n
grid[0,0] is SW, grid[-1,-1] is NE \n
grid[0,-1] is SE, grid[-1,0] is NW
- Gauss grids:\n
grid[0,:Nj] is first (Northern) band of latitude, masked after
Nj = number of longitudes for latitude j \n
grid[-1,:Nj] is last (Southern) band of latitude (idem).
"""
if not self.spectral and subzone and \
'LAMzone' in self.geometry.grid.keys():
data = self.geometry.extract_subzone(self.data, subzone)
else:
data = self.data
return data
[docs] def setdata(self, data):
"""
Sets field data, checking *data* to be:
- 2D if not spectral,
- 1D if spectral.
"""
if not self.spectral and len(numpy.shape(data)) != 2:
raise epygramError("gridpoint data must be 2D array.")
elif self.spectral and len(numpy.shape(data)) != 1:
raise epygramError("spectral data must be 1D array.")
super(H2DField, self).setdata(data)
[docs] def getvalue_ij(self, i, j):
"""
Returns the value of the field on point of indices (*i, j*).
Take care (*i, j*) is python-indexing, ranging from 0 to dimension - 1.
"""
if not self.geometry.point_is_inside_domain_ij(i, j):
raise ValueError("point is out of field domain.")
value = numpy.copy(self.data[j, i])
if numpy.shape(value) in ((1,), ()):
value = float(value)
return value
[docs] def getvalue_ll(self, lon, lat,
interpolation='nearest',
neighborinfo=False):
"""
Returns the value of the field on point of coordinates (*lon, lat*): \n
- if *interpolation == 'nearest'* (default), returns the value of the
nearest neighboring gridpoint;
- if *interpolation == 'linear'*, computes and returns the field value
with linear spline interpolation;
- if *interpolation == 'cubic'*, computes and returns the field value
with cubic spline interpolation.
If *neighborinfo* is set to **True**, returns a tuple
*(value, (lon, lat))*, with *(lon, lat)* being the actual coordinates
of the neighboring gridpoint (only for *interpolation == 'nearest'*).
*lon* and *lat* may be longer than 1.
Warning: for interpolation on Gauss geometries, requires the
:mod:`pyproj` module.
"""
if self.spectral:
raise epygramError("field must be gridpoint to get value of a\
lon/lat point.")
try:
N = len(lon)
except Exception:
N = 1
lon = numpy.array([lon])
lat = numpy.array([lat])
if interpolation == 'nearest':
(ri, rj) = self.geometry.nearest_points(lon, lat, 'nearest')
value = self.getvalue_ij(ri, rj)
if neighborinfo:
(lon, lat) = self.geometry.ij2ll(ri, rj)
if numpy.shape(lon) in ((1,), ()):
lon = float(lon)
lat = float(lat)
value = (value, (lon, lat))
elif interpolation in ('linear', 'cubic'):
from scipy.interpolate import interp2d
nvalue = numpy.zeros(N)
for n in range(N):
points = self.geometry.nearest_points(lon[n], lat[n],
interpolation)
lons = [self.geometry.ij2ll(i, j)[0] for (i, j) in points]
lats = [self.geometry.ij2ll(i, j)[1] for (i, j) in points]
values = [self.data[j, i] for (i, j) in points]
f = interp2d(lons, lats, values, kind=interpolation)
value = f(lon[n], lat[n])
nvalue[n] = value
value = nvalue
try:
value = float(value)
except ValueError:
pass
return copy.copy(value)
###################
# PRE-APPLICATIVE #
###################
# (but useful and rather standard) !
# [so that, subject to continuation through updated versions,
# including suggestions/developments by users...]
[docs] def plotfield(self,
subzone=None,
label=None,
gisquality='i',
specificproj=None,
existingbasemap=None,
minmax=None,
graphicmode='colorshades',
levelsnumber=21,
colormap='jet',
center_cmap_on_0=False,
drawrivers=False,
llgridstep=None,
colorbar='right',
minmax_in_title=True,
departments=False,
existingfigure=None,
pointsize=20,
contourcolor='blue',
contourwidth=1,
contourlabel=True):
"""
Makes a simple plot of the field, with a number of options.
Requires :mod:`matplotlib`
Options: \n
- *subzone*: among ('C', 'CI'), for LAM fields only, plots the data
resp. on the C or C+I zone. \n
Default is no subzone, i.e. the whole field.
- *gisquality*: among ('c', 'l', 'i', 'h', 'f') -- by increasing
quality. Defines the quality for GIS elements (coastlines,
countries boundaries...).
Cf. 'basemap' doc for more details.
- *specificproj*: among ('kav7', 'cyl', 'ortho', 'nsperXXXX'),
to represent data on these kinds of projections. \n
Else, the default projection is the grid projection for LAM fields,
'moll' for global fields. \n
Cf. 'basemap' doc for more details. The XXXX of 'nsperXXXX' defines
the satellite height in km.
- *existingbasemap*: as making Basemap is the most time-consuming step,
it is possible to bring an already existing basemap object.
In that case, all the above options are ignored,
overwritten by those of the *existingbasemap*.
- *label*: title for the plot. Default is field identifier.
- *minmax*: defines the min and max values for the plot colorbar. \n
Syntax: [min, max]. [0.0, max] also works. Default is min/max of the
field.
- *graphicmode*: among ('colorshades', 'contourlines', 'points').
- *levelsnumber*: number of levels for contours and colorbar.
- *colormap*: name of the **matplotlib** colormap to use.
- *center_cmap_on_0*: aligns the colormap center on the value 0.
- *drawrivers*: to add rivers on map.
- *colorbar*: if *False*, hide colorbar the plot; else, befines the
colorbar position, among ('bottom', 'right'). Defaults to 'right'.
- *llgridstep*: specific step between two lon/lat gridlines, in
degrees. \n
Default depends on domain size.
- *minmax_in_title*: if True and minmax != None, adds min and max
values in title.
- *departments*: if True, adds the french departments on map (instead
of countries).
- *existingfigure*: to plot the field over an existing figure (e.g.
contourlines over colorshades).
Be aware that no check is done between the *existingfigure* basemap
and either the *existingbasemap* or the one that is created from the
field geometry: there might be inconsistency.
- *pointsize*: size of points for *graphicmode* == 'points'.
- *contourcolor*: color or colormap to be used for 'contourlines'
graphicmode. It can be either a legal html color name, or a colormap
name.
- *contourwidth*: width of contours for 'contourlines' graphicmode.
- *contourlabel*: displays labels on contours.
This method uses (hence requires) 'matplotlib' and 'basemap' libraries.
"""
import matplotlib.pyplot as plt
from matplotlib.colors import cnames
plt.rc('font', family='serif')
plt.rc('figure', figsize=config.plotsizes)
if colormap not in plt.colormaps():
if colormap in config.colormaps:
util.add_cmap(colormap)
else:
from mpl_toolkits.basemap import cm
if colormap in cm.datad.keys():
plt.register_cmap(name=colormap, cmap=cm.__dict__[colormap])
if graphicmode == 'contourlines':
if contourcolor in cnames:
colormap = None
else:
if contourcolor not in plt.colormaps():
util.add_cmap(contourcolor)
colormap = contourcolor
contourcolor = None
if self.spectral:
raise epygramError("please convert to gridpoint with sp2gp()\
method before plotting.")
if existingfigure == None:
f = plt.figure()
else:
f = existingfigure
academic = self.geometry.name == 'academic'
if existingbasemap == None:
bm = self.geometry.make_basemap(gisquality=gisquality,
subzone=subzone,
specificproj=specificproj)
else:
bm = existingbasemap
if academic:
(lons, lats) = self.geometry.get_lonlat_grid(subzone=subzone)
x, y = lons, lats
else:
bm.drawcoastlines()
if departments:
import json
with open(config.installdir + '/data/departments.json', 'r') as dp:
depts = json.load(dp)[1]
for d in range(len(depts)):
for part in range(len(depts[d])):
dlon = depts[d][part][0]
dlat = depts[d][part][1]
(x, y) = bm(dlon, dlat)
bm.plot(x, y, color='0.25')
else:
bm.drawcountries()
if drawrivers:
bm.drawrivers(color='blue')
meridians = numpy.arange(0, 360, 20)
parallels = numpy.arange(-90, 90, 20)
if self.geometry.rectangular_grid:
corners = self.geometry.gimme_corners_ll(subzone=subzone)
minlon = min([c[0] for c in (corners['ul'], corners['ll'])])
maxlon = max([c[0] for c in (corners['ur'], corners['lr'])])
maxlon = util.degrees_nearest_mod(maxlon, minlon)
minlat = min([c[1] for c in (corners['ll'], corners['lr'])])
maxlat = max([c[1] for c in (corners['ul'], corners['ur'])])
if llgridstep:
meridians = numpy.arange(int(minlon) - llgridstep,
maxlon + llgridstep,
llgridstep)
parallels = numpy.arange(int(minlat) - llgridstep,
maxlat + llgridstep,
llgridstep)
else:
if abs(maxlon - minlon) < 15.:
meridians = numpy.arange(0, 360, 2)
elif abs(maxlon - minlon) < 25.:
meridians = numpy.arange(0, 360, 5)
elif abs(maxlon - minlon) < 90.:
meridians = numpy.arange(0, 360, 10)
if abs(maxlat - minlat) < 15.:
parallels = numpy.arange(-90, 90, 2)
elif abs(maxlat - minlat) < 25.:
parallels = numpy.arange(-90, 90, 5)
elif abs(maxlat - minlat) < 90.:
parallels = numpy.arange(-90, 90, 10)
elif llgridstep:
meridians = numpy.arange(0, 360, llgridstep)
parallels = numpy.arange(-90, 90, llgridstep)
if bm.projection in ('ortho', 'nsper'):
bm.drawparallels(parallels, labels=[False, False, False, False])
else:
bm.drawparallels(parallels, labels=[True, False, False, False])
if bm.projection in ('spstere', 'npstere'):
bm.drawmeridians(meridians, labels=[True, False, False, True])
elif bm.projection in ('ortho', 'moll', 'nsper'):
bm.drawmeridians(meridians, labels=[False, False, False, False])
else:
bm.drawmeridians(meridians, labels=[False, False, False, True])
bm.drawmeridians([0], labels=[False] * 4, linewidth=1,
dashes=[10, 1])
bm.drawparallels([0], labels=[False] * 4, linewidth=1,
dashes=[10, 1])
(lons, lats) = self.geometry.get_lonlat_grid(subzone=subzone)
x, y = bm(lons, lats)
data = numpy.ma.masked_outside(self.getdata(subzone=subzone),
- config.mask_outside,
config.mask_outside)
m = data.min()
M = data.max()
if minmax != None:
if minmax_in_title:
minmax_in_title = '(min: ' + \
'{: .{precision}{type}}'.format(m, type='E', precision=3) + \
' // max: ' + \
'{: .{precision}{type}}'.format(M, type='E', precision=3) + ')'
try:
m = float(minmax[0])
except ValueError:
m = data.min()
try:
M = float(minmax[1])
except ValueError:
M = data.max()
else:
minmax_in_title = ''
if abs(m - M) > config.epsilon:
levels = numpy.linspace(m, M, levelsnumber)
vmin = vmax = None
if center_cmap_on_0:
vmax = max(abs(m), M)
vmin = -vmax
else:
raise epygramError("cannot plot uniform field.")
L = int((levelsnumber - 1) // 15) + 1
hlevels = [levels[l] for l in range(len(levels) - (L / 3 + 1)) if
l % L == 0] + [levels[-1]]
if graphicmode == 'colorshades':
if not self.geometry.rectangular_grid:
xf = numpy.ma.masked_where(data.mask, x).compressed()
yf = numpy.ma.masked_where(data.mask, y).compressed()
zf = data.compressed()
pf = bm.contourf(xf, yf, zf, levels,
tri=True)
else:
pf = bm.contourf(x, y, data, levels,
cmap=colormap, vmin=vmin, vmax=vmax)
if colorbar:
if academic:
cb = f.colorbar(pf, ticks=hlevels, ax=bm)
else:
cb = bm.colorbar(pf, location=colorbar, pad="5%",
ticks=hlevels)
if minmax_in_title != '':
cb.set_label(minmax_in_title)
elif graphicmode == 'contourlines':
if not self.geometry.rectangular_grid:
xf = x.compressed()
yf = y.compressed()
zf = data.compressed()
pf = bm.contour(xf, yf, zf, levels=levels, colors=contourcolor,
cmap=colormap, linewidths=contourwidth,
tri=True)
#FIXME: problem of duplicate points with arpege grid
else:
pf = bm.contour(x, y, data, levels=levels, colors=contourcolor,
cmap=colormap, linewidths=contourwidth)
if contourlabel:
f.axes[0].clabel(pf, colors=contourcolor, cmap=colormap)
elif graphicmode == 'points':
xf = x.flatten()
yf = y.flatten()
zf = numpy.ma.masked_outside(data.flatten(), m, M)
pf = bm.scatter(xf, yf, c=zf, s=pointsize, marker=',', linewidths=0,
cmap=colormap, vmin=m, vmax=M)
if colorbar:
if academic:
cb = f.colorbar(pf, ticks=hlevels, ax=bm)
else:
cb = bm.colorbar(pf, location=colorbar, pad="5%",
ticks=hlevels)
if minmax_in_title != '':
cb.set_label(minmax_in_title)
if label == None:
f.axes[0].set_title(str(self.fid) + "\n" + str(self.validity.get()))
else:
f.axes[0].set_title(label)
return f
[docs] def stats(self, subzone=None):
"""
Computes some basic statistics on the field, as a dict containing:
{'min', 'max', 'mean', 'std', 'quadmean', 'nonzero'}.
See each of these methods for details.
- *subzone*: optional, among ('C', 'CI'), for LAM fields only, plots
the data resp. on the C or C+I zone. \n
Default is no subzone, i.e. the whole field.
"""
return {'min':self.min(subzone=subzone),
'max':self.max(subzone=subzone),
'mean':self.mean(subzone=subzone),
'std':self.std(subzone=subzone),
'quadmean':self.quadmean(subzone=subzone),
'nonzero':self.nonzero(subzone=subzone)}
[docs] def min(self, subzone=None):
"""Returns the minimum value of data."""
data = self.getdata(subzone=subzone)
return numpy.ma.masked_outside(data,
- config.mask_outside,
config.mask_outside).min()
[docs] def max(self, subzone=None):
"""Returns the maximum value of data."""
data = self.getdata(subzone=subzone)
return numpy.ma.masked_outside(data,
- config.mask_outside,
config.mask_outside).max()
[docs] def mean(self, subzone=None):
"""Returns the mean value of data."""
data = self.getdata(subzone=subzone)
return numpy.ma.masked_outside(data,
- config.mask_outside,
config.mask_outside).mean()
[docs] def std(self, subzone=None):
"""Returns the standard deviation of data."""
data = self.getdata(subzone=subzone)
return numpy.ma.masked_outside(data,
- config.mask_outside,
config.mask_outside).std()
[docs] def quadmean(self, subzone=None):
"""Returns the quadratic mean of data."""
data = self.getdata(subzone=subzone)
return numpy.sqrt((numpy.ma.masked_outside(data,
- config.mask_outside,
config.mask_outside) ** 2).mean())
[docs] def nonzero(self, subzone=None):
"""
Returns the number of non-zero values (whose absolute
value > config.epsilon).
"""
data = self.getdata(subzone=subzone)
return numpy.count_nonzero(abs(numpy.ma.masked_outside(data,
- config.mask_outside,
config.mask_outside)) > config.epsilon)
[docs] def dctspectrum(self, subzone=None):
"""
Returns the DCT spectrum of the field, as a
:class:`epygram.spectra.Spectrum` instance.
"""
import epygram.spectra as esp
variances = esp.dctspectrum(self.getdata(subzone=subzone))
spectrum = esp.Spectrum(variances[1:],
name=str(self.fid),
resolution=self.geometry.grid['X_resolution'] / 1000.,
mean2=variances[0])
return spectrum
[docs] def global_shift_center(self, longitude_shift):
"""
For global RegLLGeometry grids only !
Shifts the center of the geometry (and the data accordingly) by
*longitude_shift* (in degrees). *longitude_shift* has to be a multiple
of the grid's resolution in longitude.
"""
if self.geometry.name != 'regular_lonlat':
raise epygramError("only for regular lonlat geometries.")
self.geometry.global_shift_center(longitude_shift)
n = longitude_shift / self.geometry.grid['X_resolution'].get('degrees')
self.setdata(numpy.concatenate((self.data[:, n:], self.data[:, 0:n]),
axis=1))
[docs] def what(self, out, vertical_geometry=True, cumulativeduration=True):
"""
Writes in file a summary of the field.
Args: \n
- *out*: the output open file-like object (duck-typing: *out*.write()
only is needed).
- *vertical_geometry*: if True, writes the vertical geometry of the
field.
"""
firstcolumn_width = 50
secondcolumn_width = 16
sepline = '{:-^{width}}'.format('', width=firstcolumn_width + \
secondcolumn_width + 1) + '\n'
def write_formatted(dest, label, value):
dest.write('{:<{width}}'.format(label, width=firstcolumn_width) \
+ ':' \
+ '{:>{width}}'.format(str(value), width=secondcolumn_width) \
+ '\n')
v_geometry = self.geometry.vcoordinate
geometryname = self.geometry.name
grid = self.geometry.grid
dimensions = self.geometry.dimensions
validity = self.validity
out.write("################\n")
out.write("### VALIDITY ###\n")
out.write("################\n")
write_formatted(out, "Validity", validity.get())
write_formatted(out, "Basis", validity.getbasis())
write_formatted(out, "Term", validity.term())
if cumulativeduration:
write_formatted(out, "Duration for cumulative quantities",
validity.cumulativeduration())
write_formatted(out, "Kind of producting process", self.processtype)
out.write(sepline)
out.write("\n")
out.write("###########################\n")
out.write("### HORIZONTAL GEOMETRY ###\n")
out.write("###########################\n")
write_formatted(out, "Rectangular grid ( = LAM or reg. Lon/Lat)",
self.geometry.rectangular_grid)
if self.geometry.rectangular_grid:
if 'LAMzone' in grid.keys():
write_formatted(out, "Zone", grid['LAMzone'])
write_formatted(out, "Total points in X (NDLON)",
dimensions['X'])
write_formatted(out, "Total points in Y (NDGL)",
dimensions['Y'])
if grid['LAMzone'] == 'CIE':
write_formatted(out, "Points of C+I in X",
dimensions['X_CIzone'])
write_formatted(out, "Points of C+I in Y",
dimensions['Y_CIzone'])
write_formatted(out,
"Low-left X offset for I zone (NDLUN-1)",
dimensions['X_CIoffset'])
write_formatted(out,
"Low-left Y offset for I zone (NDGUN-1)",
dimensions['Y_CIoffset'])
write_formatted(out, "Width of I strip in X",
dimensions['X_Iwidth'])
write_formatted(out, "Width of I strip in Y",
dimensions['Y_Iwidth'])
elif grid['LAMzone'] == 'CI':
write_formatted(out, "Width of I strip in X",
dimensions['X_Iwidth'])
write_formatted(out, "Width of I strip in Y",
dimensions['Y_Iwidth'])
if self.spectral:
write_formatted(out, "Truncation in X (NMSMAX)",
self.spectral_geometry.truncation['in_X'])
write_formatted(out, "Truncation in Y (NSMAX)",
self.spectral_geometry.truncation['in_Y'])
else:
write_formatted(out, "Total points in X (NDLON)",
dimensions['X'])
write_formatted(out, "Total points in Y (NDGL)",
dimensions['Y'])
out.write(sepline)
if self.geometry.projected_geometry:
projection = self.geometry.projection
# projection
projmap = {'regular_lonlat':'Regular Lon/Lat',
'lambert':'Lambert (conformal conic)',
'mercator':'Mercator',
'polar_stereographic':'Polar Stereographic'}
if geometryname == 'polar_stereographic':
write_formatted(out, "Kind of projection",
projection['pole'] + projmap[geometryname])
else:
write_formatted(out, "Kind of projection",
projmap[geometryname])
write_formatted(out, "Sinus of Reference Latitude",
projection['reference_lat'].get('cos_sin')[1])
write_formatted(out, "Reference Longitude in deg (ELON0)",
projection['reference_lon'].get('degrees'))
write_formatted(out, "Reference Latitude in deg (ELAT0)",
projection['reference_lat'].get('degrees'))
(lons, lats) = self.geometry.get_lonlat_grid(subzone='CI')
corners = self.geometry.gimme_corners_ll(subzone='CI')
write_formatted(out, "Center Longitude (of C+I) in deg (ELONC)",
projection['center_lon'].get('degrees'))
write_formatted(out, "Center Latitude (of C+I) in deg (ELATC)",
projection['center_lat'].get('degrees'))
write_formatted(out, "Resolution in X, in metres (EDELX)",
grid['X_resolution'])
write_formatted(out, "Resolution in Y, in metres (EDELY)",
grid['Y_resolution'])
write_formatted(out,
"Domain width (of C+I) in X, in metres (ELX)",
grid['X_resolution'] * dimensions['X_CIzone'])
write_formatted(out,
"Domain width (of C+I) in Y, in metres (ELY)",
grid['Y_resolution'] * dimensions['Y_CIzone'])
write_formatted(out, "Max Longitude (of C+I) in deg", lons.max())
write_formatted(out, "Min Longitude (of C+I) in deg", lons.min())
write_formatted(out, "Max Latitude (of C+I) in deg", lats.max())
write_formatted(out, "Min Latitude (of C+I) in deg", lats.min())
write_formatted(out,
"Low-Left corner (of C+I) Longitude in deg",
corners['ll'][0])
write_formatted(out,
"Low-Left corner (of C+I) Latitude in deg",
corners['ll'][1])
write_formatted(out,
"Low-Right corner (of C+I) Longitude in deg",
corners['lr'][0])
write_formatted(out,
"Low-Right corner (of C+I) Latitude in deg",
corners['lr'][1])
write_formatted(out,
"Upper-Left corner (of C+I) Longitude in deg",
corners['ul'][0])
write_formatted(out,
"Upper-Left corner (of C+I) Latitude in deg",
corners['ul'][1])
write_formatted(out,
"Upper-Right corner (of C+I) Longitude in deg",
corners['ur'][0])
write_formatted(out,
"Upper-Right corner (of C+I) Latitude in deg",
corners['ur'][1])
elif geometryname == 'regular_lonlat':
(lons, lats) = self.geometry.get_lonlat_grid()
corners = self.geometry.gimme_corners_ll()
write_formatted(out, "Kind of Geometry", 'Regular Lon/Lat')
write_formatted(out, "Center Longitude in deg (ELONC)",
grid['center_lon'].get('degrees'))
write_formatted(out, "Center Latitude in deg (ELATC)",
grid['center_lat'].get('degrees'))
write_formatted(out, "Resolution in X, in deg (EDELX)",
grid['X_resolution'].get('degrees'))
write_formatted(out, "Resolution in Y, in deg (EDELY)",
grid['Y_resolution'].get('degrees'))
write_formatted(out, "Domain width in X, in deg (ELX)",
grid['X_resolution'].get('degrees') * dimensions['X'])
write_formatted(out, "Domain width in Y, in deg (ELY)",
grid['Y_resolution'].get('degrees') * dimensions['Y'])
write_formatted(out, "Max Longitude in deg", lons.max())
write_formatted(out, "Min Longitude in deg", lons.min())
write_formatted(out, "Max Latitude in deg", lats.max())
write_formatted(out, "Min Latitude in deg", lats.min())
write_formatted(out, "Low-Left corner Longitude in deg",
corners['ll'][0])
write_formatted(out, "Low-Left corner Latitude in deg",
corners['ll'][1])
write_formatted(out, "Low-Right corner Longitude in deg",
corners['lr'][0])
write_formatted(out, "Low-Right corner Latitude in deg",
corners['lr'][1])
write_formatted(out, "Upper-Left corner Longitude in deg",
corners['ul'][0])
write_formatted(out, "Upper-Left corner Latitude in deg",
corners['ul'][1])
write_formatted(out, "Upper-Right corner Longitude in deg",
corners['ur'][0])
write_formatted(out, "Upper-Right corner Latitude in deg",
corners['ur'][1])
elif geometryname == 'academic':
write_formatted(out, "Kind of Geometry", 'Academic')
else:
# global
gridmap = {'reduced_gauss':'Reduced Gauss',
'rotated_reduced_gauss':'Rotated Reduced Gauss'}
write_formatted(out, "Grid", gridmap[geometryname])
if geometryname == 'rotated_reduced_gauss':
write_formatted(out, "Pole Longitude",
grid['pole_lon'].get('degrees'))
write_formatted(out, "Pole Latitude",
grid['pole_lat'].get('degrees'))
write_formatted(out, "Dilatation coefficient",
grid['dilatation_coef'])
write_formatted(out, "Number of latitudes",
dimensions['lat_number'])
write_formatted(out, "Maximum number of longitudes on a parallel",
dimensions['max_lon_number'])
if self.spectral:
write_formatted(out, "Truncation",
self.spectral_geometry.truncation['max'])
out.write(sepline)
out.write("\n")
if vertical_geometry:
out.write("#########################\n")
out.write("### VERTICAL GEOMETRY ###\n")
out.write("#########################\n")
for k, v in v_geometry.items():
write_formatted(out, k, v)
out.write(sepline)
out.write("\n")
#############
# OPERATORS #
#############
def _check_operands(self, other):
"""
Internal method to check compatibility of terms in operations on fields.
"""
if isinstance(other, self.__class__):
if self.spectral != other.spectral:
raise epygramError("cannot operate a spectral field with a\
non-spectral field.")
if self.geometry.dimensions != other.geometry.dimensions:
raise epygramError("operations on fields cannot be done if\
fields do not share their gridpoint\
dimensions.")
if self.spectral_geometry != other.spectral_geometry:
raise epygramError("operations on fields cannot be done if\
fields do not share their spectral\
geometry.")
else:
super(H2DField, self)._check_operands(other)
def __add__(self, other):
"""
Definition of addition, 'other' being:
- a scalar (integer/float)
- another Field of the same subclass.
Returns a new Field whose data is the resulting operation,
with 'fid' = {'op':'+'} and null validity.
"""
newfield = super(H2DField, self)._add(other,
geometry=self.geometry,
spectral_geometry=self.spectral_geometry)
return newfield
def __mul__(self, other):
"""
Definition of multiplication, 'other' being:
- a scalar (integer/float)
- another Field of the same subclass.
Returns a new Field whose data is the resulting operation,
with 'fid' = {'op':'*'} and null validity.
"""
newfield = super(H2DField, self)._mul(other,
geometry=self.geometry,
spectral_geometry=self.spectral_geometry)
return newfield
def __sub__(self, other):
"""
Definition of substraction, 'other' being:
- a scalar (integer/float)
- another Field of the same subclass.
Returns a new Field whose data is the resulting operation,
with 'fid' = {'op':'-'} and null validity.
"""
newfield = super(H2DField, self)._sub(other,
geometry=self.geometry,
spectral_geometry=self.spectral_geometry)
return newfield
def __div__(self, other):
"""
Definition of division, 'other' being:
- a scalar (integer/float)
- another Field of the same subclass.
Returns a new Field whose data is the resulting operation,
with 'fid' = {'op':'/'} and null validity.
"""
newfield = super(H2DField, self)._div(other,
geometry=self.geometry,
spectral_geometry=self.spectral_geometry)
return newfield