#!/usr/bin/env python
# -*- coding: utf-8 -*-
from footprints import proxy as fpx, FPDict
from epygram import config, SpectralGeometry, epygramError, H2DField, util
import numpy
[docs]def make_vector(fX, fY):
"""
Creates a new :class:`epygram.H2DVectorField` from two
:class:`epygram.H2DField` *fX, fY* representing resp.
the X and Y components of the vector in the field geometry.
"""
if not isinstance(fX, H2DField) or not isinstance(fY, H2DField):
raise epygramError("'fX', 'fY' must be H2DField.")
if fX.geometry.dimensions != fY.geometry.dimensions:
raise epygramError("'fX', 'fY' must be share their gridpoint dimensions.")
if fX.spectral_geometry != fY.spectral_geometry:
raise epygramError("'fX', 'fY' must be share their spectral geometry.")
f = fpx.field(fid=FPDict({'op':'make_vector()'}),
geometry=fX.geometry,
validity=fX.validity,
spectral_geometry=fX.spectral_geometry,
processtype=fX.processtype,
vector=True)
if fX.spectral:
data = numpy.zeros((len(fX.data), 2))
data[:,0] = fX.data[:]
data[:,1] = fY.data[:]
else:
if fX.geometry.rectangular_grid:
data = numpy.zeros((fX.geometry.dimensions['Y'],
fX.geometry.dimensions['X'],
2))
else:
data = numpy.ma.masked_all((fX.geometry.dimensions['lat_number'],
fX.geometry.dimensions['max_lon_number'],
2))
data[:,:,0] = fX.data[:,:]
data[:,:,1] = fY.data[:,:]
f.setdata(data)
return f
[docs]class H2DVectorField(H2DField):
"""
Horizontal 2-Dimensions Vector field class.
This is a H2DField whose data has 2 components, representing the X and Y components of a
vector projected on its geometry (the grid axes).
The attribute H2DVectorField.data is hence an array dimensioned as (nY, nX, 2).
"""
_collector = ('field',)
_footprint = dict(
attr = dict(
vector = dict(
type = bool,
values = set([True]))
)
)
##############
# ABOUT DATA #
##############
[docs] def sp2gp(self):
"""
Transforms the spectral field into gridpoint, according to its spectral geometry.
Replaces data on site.
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]
dataspX = self.data[:,0].flatten()
dataspY = self.data[:,1].flatten()
data2dX = self.geometry.reshape_data(self.spectral_geometry.sp2gp(dataspX, gpdims))
data2dY = self.geometry.reshape_data(self.spectral_geometry.sp2gp(dataspY, gpdims))
if self.geometry.rectangular_grid:
data2d = numpy.zeros((gpdims['Y'], gpdims['X'], 2))
else:
data2d = numpy.ma.masked_all((self.geometry.dimensions['lat_number'],
self.geometry.dimensions['max_lon_number'],
2))
data2d[:,:,0] = data2dX[:,:]
data2d[:,:,1] = data2dY[:,:]
self._attributes['spectral_geometry'] = None
self.setdata(data2d)
[docs] def gp2sp(self, spectral_geometry=None):
"""
Transforms the gridpoint field into spectral space, according to the spectral geometry
mandatorily passed as argument.
Replaces data on site.
The spectral transform subroutine is actually included in the spectral geometry's
*gp2sp()* method.
"""
if not self.spectral:
if spectral_geometry == None or 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]
data1dX = self.data[:,:,0].flatten()
data1dY = self.data[:,:,1].flatten()
else:
# global
gpdims = {}
for dim in ['lat_number', 'lon_number_by_lat']:
gpdims[dim] = self.geometry.dimensions[dim]
data1dX = self.data[:,:,0].compressed()
data1dY = self.data[:,:,1].compressed()
dataspX = spectral_geometry.gp2sp(data1dX, gpdims)
dataspY = spectral_geometry.gp2sp(data1dY, gpdims)
datasp = numpy.zeros((len(dataspX), 2))
datasp[:,0] = dataspX[:]
datasp[:,1] = dataspY[:]
self._attributes['spectral_geometry'] = spectral_geometry
self.setdata(spectral_geometry.gp2sp(datasp, gpdims))
[docs] def getdata(self, subzone=None):
"""
Returns the field data, with 2D shape if the field is not spectral,
1D if spectral. The third (resp. second) dimension specifies the X/Y component.
- 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: (x (0) being the X component, y (1) the Y one) \n
- Rectangular grids:\n
grid[0,0,x] is SW, grid[-1,-1,x] is NE \n
grid[0,-1,x] is SE, grid[-1,0,x] is NW
- Gauss grids:\n
grid[0,:Nj,x] is first (Northern) band of latitude, masked after Nj = number of longitudes for latitude j \n
grid[-1,:Nj,x] is last (Southern) band of latitude (idem).
"""
if not self.spectral and subzone and 'LAMzone' in self.geometry.grid.keys():
dataX = self.geometry.extract_subzone(self.data[:,:,0], subzone)
dataY = self.geometry.extract_subzone(self.data[:,:,1], subzone)
if subzone == 'CI':
dimX = self.geometry.dimensions['X_CIzone']
dimY = self.geometry.dimensions['Y_CIzone']
elif subzone == 'C':
dimX = self.geometry.dimensions['X_Czone']
dimY = self.geometry.dimensions['Y_Czone']
data = numpy.zeros((dimY, dimX, 2))
data[:,:,0] = dataX[:,:]
data[:,:,1] = dataY[:,:]
else:
data = self.data
return data
[docs] def setdata(self, data):
"""
Sets data, checking it to be:
- 2D (+1) if not spectral,
- 1D (+1) if spectral.
"""
if not self.spectral and len(numpy.shape(data)) != 3:
raise epygramError("gridpoint data must be 2D (+1) array.")
elif self.spectral and len(numpy.shape(data)) != 2:
raise epygramError("spectral data must be 1D (+1) array.")
super(H2DField, self).setdata(data)
[docs] def to_module(self):
"""
Returns a :class:`epygram.H2DField` whose data is the module of the
Vector field.
"""
if self.spectral:
fieldcopy = self.deepcopy()
fieldcopy.sp2gp()
datagp = fieldcopy.data
else:
datagp = self.data
module = numpy.sqrt(datagp[:,:,0]**2 + datagp[:,:,1]**2)
f = fpx.field(geometry=self.geometry,
fid=FPDict({'op':'H2DVectorField.to_module()'}),
validity=self.validity,
processtype=self.processtype)
f.setdata(module)
if self.spectral:
f.gp2sp(self.spectral_geometry)
return f
###################
# PRE-APPLICATIVE #
###################
# (but useful and rather standard) !
# [so that, subject to continuation through updated versions,
# including suggestions/developments by users...]
[docs] def plotfield(self, existingfigure=None, subzone=None, label=None,
gisquality='i', specificproj=None, existingbasemap=None,
drawrivers=False, llgridstep=None, subsampling=1):
"""
Makes a simple plot of the field, with a number of options.
Options: \n
- *existingfigure*: to plot the vectors over an existing figure (e.g. 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.
- *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...).
Default is 'i'. 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.
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*.
- *drawrivers*: to add rivers on map.
- *label*: title for the plot. Default is field identifier.
- *llgridstep*: specific step between two lon/lat gridlines, in degrees. \n
Default depends on domain size.
- *subsampling*: to subsample the number of gridpoints to plot.
Ex: *subsampling* = 10 will only plot one gridpoint upon 10.
This method uses (hence requires) 'matplotlib' and 'basemap' libraries.
"""
import matplotlib.pyplot as plt
plt.rc('font', family='serif')
plt.rc('figure', figsize=config.plotsizes)
if self.spectral:
raise epygramError("please convert to gridpoint with sp2gp() method before plotting.")
if self.geometry.projected_geometry and specificproj != None:
raise epygramError("cannot plot vector from a projected geometry on another geometry: would need vector rotation.")
if existingfigure == None:
f = plt.figure()
else:
f = existingfigure
if existingbasemap == None:
bm = self.geometry.make_basemap(gisquality=gisquality, subzone=subzone, specificproj=specificproj)
else:
bm = existingbasemap
bm.drawcoastlines()
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)
if bm.projection in ('ortho', 'nsper'):
bm.drawparallels(parallels, labels=[False,False,False,False])
else:
bm.drawparallels(parallels, labels=[False,True,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)
lons = lons[::subsampling, ::subsampling]
lats = lats[::subsampling, ::subsampling]
x, y = bm(lons, lats)
data = numpy.ma.masked_outside(self.getdata(subzone=subzone),
-config.mask_outside,
config.mask_outside)
if not self.geometry.rectangular_grid:
xf = x.compressed()
yf = y.compressed()
u = data[::subsampling, ::subsampling, 0].compressed()
v = data[::subsampling, ::subsampling, 1].compressed()
else:
xf = x.flatten()
yf = y.flatten()
u = data[::subsampling, ::subsampling, 0].flatten()
v = data[::subsampling, ::subsampling, 1].flatten()
bm.barbs(xf, yf, u, v)
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 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 self.geometry.rectangular_grid:
if not self.geometry.point_is_inside_domain(self.geometry.ij2ll(i, j)):
raise ValueError("point is out of field domain.")
elif 'reduced_gauss' in self.geometry.name:
if j >= self.geometry.dimensions['lat_number'] or j < 0:
raise ValueError("j-index out of grid.")
if i >= self.geometry.dimensions['lon_number_by_lat'][j] or i < 0:
raise ValueError("i-index out of grid: " + \
str(i) + " out of [0," + \
str(self.geometry.dimensions['lon_number_by_lat'][j])+"[.")
return self.data[j, i,:]
[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'*).
Warning: for interpolation on Gauss geometries, requires the
:mod:`pyproj` module.
"""
(i, j) = self.geometry.ll2ij(lon, lat)
ri = numpy.rint(i).astype('int')
rj = numpy.rint(j).astype('int')
if interpolation == 'nearest':
value = self.getvalue_ij(ri, rj)
if neighborinfo:
(lon, lat) = self.geometry.ij2ll(ri, rj)
value = (value, (lon, lat))
elif interpolation in ('linear', 'cubic'):
from scipy.interpolate import interp2d
(lons, lats) = self.geometry.get_lonlat_grid()
if self.geometry.rectangular_grid:
if not self.geometry.point_is_inside_domain((lon, lat)):
raise ValueError("point ("+str(lon)+", "+str(lat)+") is out of field domain.")
if interpolation == 'linear':
m = 1
elif interpolation == 'cubic':
m = 2
try:
lons = lons[rj-m:rj+m+1, ri-m:ri+m+1]
lats = lats[rj-m:rj+m+1, ri-m:ri+m+1]
valuesX = self.data[rj-m:rj+m+1, ri-m:ri+m+1, 0]
valuesY = self.data[rj-m:rj+m+1, ri-m:ri+m+1, 1]
except IndexError:
raise epygramError("point ("+str(lon)+", "+str(lat)+") too \
close to field domain borders (less \
than "+str(m)+" gridpoints) for interpolation.")
else: #Gauss
from pyproj import Geod
g = Geod(ellps='sphere') # fixed earth radius...
j_n = max(0, rj-2)
j_x = min(self.geometry.dimensions['lat_number'], rj+3)
lons_inter = lons[j_n:j_x, :].flatten()
lats_inter = lats[j_n:j_x, :].flatten()
values_interX = self.data[j_n:j_x, :, 0].compressed()
values_interY = self.data[j_n:j_x, :, 1].compressed()
points_inter = len(values_interX)
dist_inter = g.inv(lons_inter, lats_inter,
lon*numpy.ones(points_inter),
lat*numpy.ones(points_inter))[2]
if interpolation == 'linear':
points_num = 9
elif interpolation == 'cubic':
points_num = 16
sorted_dists = sorted(dist_inter)[0:points_num]
lons = numpy.zeros(points_num)
lats = numpy.zeros(points_num)
valuesX = numpy.zeros(points_num)
valuesY = numpy.zeros(points_num)
for i in range(points_num):
for k in range(len(values_interX)):
if dist_inter[k] == sorted_dists[i]:
lons[i] = lons_inter[k]
lats[i] = lats_inter[k]
valuesX[i] = values_interX[k]
valuesY[i] = values_interY[k]
break
fX = interp2d(lons, lats, valuesX, kind=interpolation)
fY = interp2d(lons, lats, valuesY, kind=interpolation)
valueX = fX(lon, lat)
valueY = fY(lon, lat)
if numpy.shape(valueX) in ((1,), ()):
valueX = float(valueX)
valueY = float(valueY)
value = (valueX, valueY)
return value
[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)
if self.spectral:
dataX = data[:,0]
dataY = data[:,1]
else:
dataX = data[:,:,0]
dataY = data[:,:,1]
return (numpy.ma.masked_outside(dataX, -config.mask_outside, config.mask_outside).min(),
numpy.ma.masked_outside(dataY, -config.mask_outside, config.mask_outside).min())
[docs] def max(self, subzone=None):
"""
Returns the maximum value of data.
"""
data = self.getdata(subzone=subzone)
if self.spectral:
dataX = data[:,0]
dataY = data[:,1]
else:
dataX = data[:,:,0]
dataY = data[:,:,1]
return (numpy.ma.masked_outside(dataX, -config.mask_outside, config.mask_outside).max(),
numpy.ma.masked_outside(dataY, -config.mask_outside, config.mask_outside).max())
[docs] def mean(self, subzone=None):
"""
Returns the mean value of data.
"""
data = self.getdata(subzone=subzone)
if self.spectral:
dataX = data[:,0]
dataY = data[:,1]
else:
dataX = data[:,:,0]
dataY = data[:,:,1]
return (numpy.ma.masked_outside(dataX, -config.mask_outside, config.mask_outside).mean(),
numpy.ma.masked_outside(dataY, -config.mask_outside, config.mask_outside).mean())
[docs] def std(self, subzone=None):
"""
Returns the standard deviation of data.
"""
data = self.getdata(subzone=subzone)
if self.spectral:
dataX = data[:,0]
dataY = data[:,1]
else:
dataX = data[:,:,0]
dataY = data[:,:,1]
return (numpy.ma.masked_outside(dataX, -config.mask_outside, config.mask_outside).std(),
numpy.ma.masked_outside(dataY, -config.mask_outside, config.mask_outside).std())
[docs] def quadmean(self, subzone=None):
"""
Returns the quadratic mean of data.
"""
data = self.getdata(subzone=subzone)
if self.spectral:
dataX = data[:,0]
dataY = data[:,1]
else:
dataX = data[:,:,0]
dataY = data[:,:,1]
return (numpy.sqrt((numpy.ma.masked_outside(dataX, -config.mask_outside, config.mask_outside)**2).mean()),
numpy.sqrt((numpy.ma.masked_outside(dataY, -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)
if self.spectral:
dataX = data[:,0]
dataY = data[:,1]
else:
dataX = data[:,:,0]
dataY = data[:,:,1]
return (numpy.count_nonzero(abs(numpy.ma.masked_outside(dataX, -config.mask_outside, config.mask_outside)) > config.epsilon),
numpy.count_nonzero(abs(numpy.ma.masked_outside(dataY, -config.mask_outside, config.mask_outside)) > config.epsilon))
#############
# OPERATORS #
#############
def _check_operands(self, other):
"""
Internal method to check compatibility of terms in operations on fields.
"""
if not 'vector' in other._attributes:
raise epygramError("cannot operate a Vector field with a non-Vector one.")
else:
super(H2DVectorField, self)._check_operands(other)