Source code for epygram.geometries.H2DGeometry

#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
H2DGeometry:

Contains the classes for Horizontal 2D geometries of fields.
"""

import numpy
import math

from footprints import FootprintBase, FPDict, proxy as fpx

from epygram import epygramError, config, epylog
from epygram.util import RecursiveObject, degrees_nearest_mod, Angle



[docs]class H2DGeometry(RecursiveObject, FootprintBase): """ Handles the geometry for a Horizontal 2-Dimensions Field. Abstract mother class. """ _abstract = True _collector = ('geometry',) _footprint = dict( attr=dict( structure=dict( values=set(['H2D'])), name=dict( values=set(['lambert', 'mercator', 'polar_stereographic', 'regular_lonlat', 'rotated_reduced_gauss', 'reduced_gauss'])), grid=dict( type=FPDict, info="Handles description of the horizontal grid."), dimensions=dict( type=FPDict, info="Handles grid dimensions."), vcoordinate=dict( type=FPDict, optional=True, info="Handles vertical geometry parameters."), position_on_grid=dict( type=FPDict, optional=True, info="Position of points w/r to the (horizontal/vertical) grid:\ mass or flux points.", default=FPDict({'vertical':'__unknown__', 'horizontal':'__unknown__'}), values=set(FPDict({'vertical':v, 'horizontal':h}) for v in ['mass', 'flux', '__unknown__'] for h in ['upper-right', 'upper-left', 'lower-left', 'lower-right', 'center-left', 'center-right', 'lower-center', 'upper-center', 'center', '__unknown__'])) ) ) @property
[docs] def rectangular_grid(self): """ Is the grid rectangular ? """ return isinstance(self, H2DRectangularGridGeometry)
@property
[docs] def projected_geometry(self): """ Is the geometry a projection ? """ return 'projection' in self._attributes.keys()
def __init__(self, *args, **kwargs): """Constructor. See its footprint for arguments.""" super(H2DGeometry, self).__init__(*args, **kwargs) # Checks ! self._consistency_check() def _getoffset(self, position=None): """Returns the offset to use for this position.""" if position is not None: pos = position else: pos = self.position_on_grid['horizontal'] if pos == '__unknown__': raise epygramError("position_on_grid['horizontal'] must be\ defined.") return {'upper-right' : (.5, .5), 'upper-left' : (-.5, .5), 'lower-left' : (-.5, -.5), 'lower-right' : (.5, -.5), 'center-left' : (-.5, 0.), 'center-right': (.5, 0.), 'lower-center': (0., -.5), 'upper-center': (0., .5), 'center' : (0., 0.)}[pos]
[docs]class H2DRectangularGridGeometry(H2DGeometry): """ Handles the geometry for a rectangular Horizontal 2-Dimensions Field. Abstract. """ _abstract = True _collector = ('geometry',) _footprint = dict( attr=dict( name=dict( values=set(['lambert', 'mercator', 'polar_stereographic', 'regular_lonlat', 'academic'])) ) )
[docs] def get_lonlat_grid(self, subzone=None, position=None): """ Returns a tuple of two tables containing one the longitude of each point, the other the latitude, with 2D shape. - *subzone*: optional, among ('C', 'CI'), for LAM grids only, returns the grid resp. for the C or C+I zone off the C+I+E zone. \n Default is no subzone, i.e. the whole field. - *position*: position of lonlat grid with respect to the model cell. Defaults to position_on_grid['horizontal']. Shape of 2D data on Rectangular grids: \n - grid[0,0] is SW, grid[-1,-1] is NE \n - grid[0,-1] is SE, grid[-1,0] is NW """ Jmax = self.dimensions['Y'] Imax = self.dimensions['X'] igrid = numpy.zeros(Jmax * Imax) jgrid = numpy.zeros(Jmax * Imax) for j in range(Jmax): for i in range(Imax): igrid[j * Imax + i] = i jgrid[j * Imax + i] = j (lons, lats) = self.ij2ll(igrid, jgrid, position) lons = self.reshape_data(lons) lats = self.reshape_data(lats) if subzone and 'LAMzone' in self.grid.keys(): lons = self.extract_subzone(lons, subzone) lats = self.extract_subzone(lats, subzone) return (lons, lats)
[docs] def extract_subzone(self, data, subzone): """ Extracts the subzone C or CI from a LAM field. Args: \n - *data*: the 2D data, of dimension concording with geometry. - *subzone*: optional, among ('C', 'CI'), for LAM grids only, extracts the data resp. from the C or C+I zone off the C+I(+E) zone. """ if subzone not in ('C', 'CI'): raise epygramError("only possible values for 'subzone' are 'C' or\ 'CI'.") if not 'LAMzone' in self.grid.keys(): raise epygramError("method for LAM grids only.") if self.grid['LAMzone'] == 'CIE': # remove E-zone x1 = self.dimensions['X_CIoffset'] x2 = self.dimensions['X_CIoffset'] + self.dimensions['X_CIzone'] y1 = self.dimensions['Y_CIoffset'] y2 = self.dimensions['Y_CIoffset'] + self.dimensions['Y_CIzone'] edata = data[y1:y2, x1:x2] else: edata = data if subzone == 'C': # remove I-zone x1 = self.dimensions['X_Iwidth'] x2 = -self.dimensions['X_Iwidth'] y1 = self.dimensions['Y_Iwidth'] y2 = -self.dimensions['Y_Iwidth'] edata = edata[y1:y2, x1:x2] return edata
[docs] def reshape_data(self, data): """ Returns a 2D data reshaped from 1D, according to geometry. - *data*: the 1D data, of dimension concording with geometry. """ return data.reshape((self.dimensions['Y'], self.dimensions['X']))
[docs] def gimme_corners_ij(self, subzone=None): """ Returns the indices (i, j) of the four corners of a rectangular grid, as a dict(corner=(i, j)) with corner in: \n ll = lower-left / lr = lower-right / ur = upper-right / ul = upper-left. (0, 0) is the lower-left corner of the grid. - *subzone*: for LAM fields, returns the corners of the subzone. """ if not 'LAMzone' in self.grid.keys() or self.grid['LAMzone'] == None: ll = (0, 0) lr = (self.dimensions['X'] - 1, 0) ul = (0, self.dimensions['Y'] - 1) ur = (self.dimensions['X'] - 1, self.dimensions['Y'] - 1) elif self.grid['LAMzone'] == 'CIE': if subzone in (None, 'CIE'): ll = (0, 0) lr = (self.dimensions['X'] - 1, 0) ul = (0, self.dimensions['Y'] - 1) ur = (self.dimensions['X'] - 1, self.dimensions['Y'] - 1) elif subzone == 'CI': ll = (self.dimensions['X_CIoffset'], self.dimensions['Y_CIoffset']) lr = (self.dimensions['X_CIoffset'] + self.dimensions['X_CIzone'] - 1, self.dimensions['Y_CIoffset']) ul = (self.dimensions['X_CIoffset'], self.dimensions['Y_CIoffset'] + self.dimensions['Y_CIzone'] - 1) ur = (self.dimensions['X_CIoffset'] + self.dimensions['X_CIzone'] - 1, self.dimensions['Y_CIoffset'] + self.dimensions['Y_CIzone'] - 1) elif subzone == 'C': ll = (self.dimensions['X_CIoffset'] + self.dimensions['X_Iwidth'], self.dimensions['Y_CIoffset'] + self.dimensions['Y_Iwidth']) lr = (self.dimensions['X_CIoffset'] + self.dimensions['X_Iwidth'] + self.dimensions['X_Czone'] - 1, self.dimensions['Y_CIoffset'] + self.dimensions['Y_Iwidth']) ul = (self.dimensions['X_CIoffset'] + self.dimensions['X_Iwidth'], self.dimensions['Y_CIoffset'] + self.dimensions['Y_Iwidth'] + self.dimensions['Y_Czone'] - 1) ur = (self.dimensions['X_CIoffset'] + self.dimensions['X_Iwidth'] + self.dimensions['X_Czone'] - 1, self.dimensions['Y_CIoffset'] + self.dimensions['Y_Iwidth'] + self.dimensions['Y_Czone'] - 1) elif self.grid['LAMzone'] == 'CI': if subzone in (None, 'CI'): ll = (0, 0) lr = (self.dimensions['X_CIzone'] - 1, 0) ul = (0, self.dimensions['Y_CIzone'] - 1) ur = (self.dimensions['X_CIzone'] - 1, self.dimensions['Y_CIzone'] - 1) elif subzone == 'C': ll = (self.dimensions['X_Iwidth'], self.dimensions['Y_Iwidth']) lr = (self.dimensions['X_Iwidth'] + self.dimensions['X_Czone'] - 1, self.dimensions['Y_Iwidth']) ul = (self.dimensions['X_Iwidth'], self.dimensions['Y_Iwidth'] + self.dimensions['Y_Czone'] - 1) ur = (self.dimensions['X_Iwidth'] + self.dimensions['X_Czone'] - 1, self.dimensions['Y_Iwidth'] + self.dimensions['Y_Czone'] - 1) return {'ll':ll, 'lr':lr, 'ul':ul, 'ur':ur}
[docs] def gimme_corners_ll(self, subzone=None, position=None): """ Returns the lon/lat of the four corners of a rectangular grid, as a dict(corner=(lon, lat)) with corner in: \n ll = lower-left / lr = lower-right / ur = upper-right / ul = upper-left. - *subzone*: for LAM grids, returns the corners of the subzone. - *position*: position of corners with respect to the model cell. Defaults to position_on_grid['horizontal']. """ corners = self.gimme_corners_ij(subzone=subzone) for c in corners.keys(): i = corners[c][0] j = corners[c][1] corners[c] = self.ij2ll(i, j, position) return corners
[docs] def point_is_inside_domain_ll(self, lon, lat, margin=-0.1, subzone=None, position=None): """ Returns True if the point(s) of lon/lat coordinates is(are) inside the field. Args: \n - *lon*: longitude of point(s) in degrees. - *lat*: latitude of point(s) in degrees. - *margin*: considers the point inside if at least 'margin' points far from the border. The -0.1 default is a safety for precision errors. - *subzone*: considers only a subzone among ('C', 'CI') of the domain. - *position*: position of the grid with respect to the model cell. Defaults to position_on_grid['horizontal']. """ (Xmin, Ymin) = self.gimme_corners_ij(subzone)['ll'] (Xmax, Ymax) = self.gimme_corners_ij(subzone)['ur'] try: N = len(lon) except Exception: N = 1 if N == 1: inside = Xmin + margin <= self.ll2ij(lon, lat, position)[0] <= Xmax - margin \ and Ymin + margin <= self.ll2ij(lon, lat, position)[1] <= Ymax - margin else: inside = [] for i in range(N): inside.append(Xmin + margin <= self.ll2ij(lon[i], lat[i], position)[0] <= Xmax - margin \ and Ymin + margin <= self.ll2ij(lon[i], lat[i], position)[1] <= Ymax - margin) return inside
[docs] def point_is_inside_domain_ij(self, i, j, margin=-0.1, subzone=None): """ Returns True if the point(s) of lon/lat coordinates is(are) inside the field. Args: \n - *lon*: longitude of point(s) in degrees. - *lat*: latitude of point(s) in degrees. - *margin*: considers the point inside if at least 'margin' points far from the border. The -0.1 default is a safety for precision errors. - *subzone*: considers only a subzone among ('C', 'CI') of the domain. """ (Xmin, Ymin) = self.gimme_corners_ij(subzone)['ll'] (Xmax, Ymax) = self.gimme_corners_ij(subzone)['ur'] try: N = len(i) except Exception: N = 1 if N == 1: inside = Xmin + margin <= i <= Xmax - margin \ and Ymin + margin <= j <= Ymax - margin else: inside = [] for n in range(N): inside.append(Xmin + margin <= i[n] <= Xmax - margin \ and Ymin + margin <= j[n] <= Ymax - margin) return inside
[docs] def nearest_points(self, lon, lat, interpolation, position=None): """ Returns the (i, j) position of the points needed to perform an interpolation. This is a list of (lon, lat) tuples. Args: \n - *lon*: longitude of point in degrees. - *lat*: latitude of point in degrees. - *interpolation* can be: - 'nearest' to get the nearest point only - 'linear' to get the 2*2 points bordering the (*lon*, *lat*) position - 'cubic' to get the 4*4 points bordering the (*lon*, *lat*) position - *position*: position in the model cell of the lat lon position. Defaults to position_on_grid['horizontal']. """ if not self.point_is_inside_domain_ll(lon, lat, position=position): raise ValueError("point (" + str(lon) + ", " + str(lat) + \ ") is out of field domain.") (ic, jc) = self.ll2ij(lon, lat, position) if interpolation == 'nearest': points = [(numpy.rint(ic).astype('int'), numpy.rint(jc).astype('int'))] elif interpolation == 'linear': points = [(i, j) for i in [numpy.floor(ic).astype('int') + o for o in [0, 1]] for j in [numpy.floor(jc).astype('int') + o for o in [0, 1]]] elif interpolation == 'cubic': points = [(i, j) for i in [numpy.floor(ic).astype('int') + o for o in [-1, 0, 1, 2]] for j in [numpy.floor(jc).astype('int') + o for o in [-1, 0, 1, 2]]] else: raise epygramError("'interpolation' can only be 'nearest',\ 'linear' or 'cubic'.") for point in points: if not self.point_is_inside_domain_ij(*point): raise epygramError("point (" + str(lon) + ", " + str(lat) + \ ") too close to field domain borders.") return points[0] if interpolation == 'nearest' else points
[docs]class H2DAcademicGeometry(H2DRectangularGridGeometry): """Handles the geometry for an academic Horizontal 2-Dimensions Field.""" _collector = ('geometry',) _footprint = dict( attr=dict( name=dict( values=set(['academic'])) ) ) def _consistency_check(self): """Check that the geometry is consistent.""" grid_keys = ['LAMzone', 'X_resolution', 'Y_resolution'] if set(self.grid.keys()) != set(grid_keys): raise epygramError("grid attribute must consist in keys: " + \ str(grid_keys)) LAMzone_values = [None, 'CI', 'CIE'] if self.grid['LAMzone'] not in LAMzone_values: raise epygramError("grid['LAMzone'] must be among " + \ str(LAMzone_values)) dimensions_keys = ['X', 'Y'] if self.grid['LAMzone'] in ('CI', 'CIE'): dimensions_keys.extend(['X_Iwidth', 'Y_Iwidth', 'X_Czone', 'Y_Czone', 'X_CIzone', 'Y_CIzone']) if self.grid['LAMzone'] == 'CIE': dimensions_keys.extend(['X_CIoffset', 'Y_CIoffset']) if set(self.dimensions.keys()) != set(dimensions_keys): raise epygramError("dimensions attribute must consist in keys: " + \ str(dimensions_keys)) def _getoffset(self, position=None): """ Returns the offset to use for this position. Replaces the method defined in H2DRectangularGridGeometry to deal with 1D or 2D simulations. """ offset = super(H2DAcademicGeometry, self)._getoffset(position) if self.dimensions['X'] == 1: offset = (0, offset[1]) if self.dimensions['Y'] == 1: offset = (offset[1], 0) return offset def ij2xy(self, i, j, position=None): """ (*i, j*) being the coordinates in the 2D matrix of gridpoints, \n (*x, y*) being the coordinates in the projection. - *position*: position to return with respect to model cell. Defaults to position_on_grid['horizontal']. Note that origin of coordinates in projection is the center of the C+I domain. """ if isinstance(i, list) or isinstance(i, tuple): i = numpy.array(i) if isinstance(j, list) or isinstance(j, tuple): j = numpy.array(j) (oi, oj) = self._getoffset(position) return ((i + oi) * self.grid['X_resolution'], (j + oj) * self.grid['Y_resolution']) def xy2ij(self, x, y, position=None): """ (*x, y*) being the coordinates in the projection, \n (*i, j*) being the coordinates in the 2D matrix of gridpoints. - *position*: position represented by (x,y) with respect to model cell. Defaults to position_on_grid['horizontal']. Caution: (*i,j*) are float (the nearest grid point is the nearest integer). Note that origin of coordinates in projection is the center of the C+I domain. """ if isinstance(x, list) or isinstance(x, tuple): x = numpy.array(x) if isinstance(y, list) or isinstance(y, tuple): y = numpy.array(y) (oi, oj) = self._getoffset(position) return (x / self.grid['X_resolution'] - oi, y / self.grid['Y_resolution'] - oj) def ij2ll(self, i, j, position=None): """ (*i, j*) being the coordinates in the 2D matrix of gridpoints, \n (*lon, lat*) being the lon/lat coordinates in degrees. - *position*: lat lon position to return with respect to the model cell. Defaults to position_on_grid['horizontal']. !!! This routine has no sense for this geometry, it is identity. """ if isinstance(i, list) or isinstance(i, tuple): i = numpy.array(i) if isinstance(j, list) or isinstance(j, tuple): j = numpy.array(j) (oi, oj) = self._getoffset(position) return (i + 1 + oi, j + 1 + oj) def ll2ij(self, lon, lat, position=None): """ (*lon, lat*) being the lon/lat coordinates in degrees, \n (*i, j*) being the coordinates in the 2D matrix of gridpoints. - *position*: lat lon position with respect to the model cell. Defaults to position_on_grid['horizontal']. !!! This routine has no sense for this geometry, it is identity. """ if isinstance(lon, list) or isinstance(lon, tuple): lon = numpy.array(lon) if isinstance(lat, list) or isinstance(lat, tuple): lat = numpy.array(lat) (oi, oj) = self._getoffset(position) return (lon - 1 - oi, lat - 1 - oj) def ll2xy(self, lon, lat): """ (*lon, lat*) being the lon/lat coordinates in degrees, \n (*x, y*) being the coordinates in the projection. Note that origin of coordinates in projection is the center of the C+I domain. """ if isinstance(lon, list) or isinstance(lon, tuple): lon = numpy.array(lon) if isinstance(lat, list) or isinstance(lat, tuple): lat = numpy.array(lat) return self.ij2xy(*self.ll2ij(lon, lat)) def xy2ll(self, x, y): """ (*x, y*) being the coordinates in the projection, \n (*lon, lat*) being the lon/lat coordinates in degrees. - *position*: lat lon position with respect to the model cell. Defaults to position_on_grid['horizontal']. Note that origin of coordinates in projection is the center of the C+I domain. """ if isinstance(x, list) or isinstance(x, tuple): x = numpy.array(x) if isinstance(y, list) or isinstance(y, tuple): y = numpy.array(y) return self.ij2ll(*self.xy2ij(x, y)) def distance(self, end1, end2): """ Computes the distance between two points along a straight line in the geometry. *end1* must be a tuple (lon, lat). *end2* must be a tuple (lon, lat). """ return numpy.sqrt(((end1[0] - end2[0]) * self.grid['X_resolution']) ** 2 + \ ((end1[1] - end2[1]) * self.grid['Y_resolution']) ** 2) def linspace(self, end1, end2, num): """ Returns evenly spaced points over the specified interval. Points are lined up in the geometry. *end1* must be a tuple (lon, lat). *end2* must be a tuple (lon, lat). *num* is the number of points, including point1 and point2. """ if num < 2: raise epygramError("'num' must be at least 2.") return zip(numpy.linspace(end1[0], end2[0], num=num), \ numpy.linspace(end1[1], end2[1], num=num)) def resolution_ll(self, lon, lat): """ Returns the local resolution at the nearest point of lon/lat. It's the distance between this point and its closest neighbour. *point* must be a tuple (lon, lat). """ return min(self.grid['X_resolution'], self.grid['Y_resolution']) def resolution_ij(self, i, j): """ Returns the distance to the nearest point of (i,j) point. (*i, j*) being the coordinates in the 2D matrix of gridpoints. """ return min(self.grid['X_resolution'], self.grid['Y_resolution']) def azimuth(self, end1, end2): """ Initial bearing from *end1* to *end2* points in geometry. *end1* must be a tuple (lon, lat). *end2* must be a tuple (lon, lat). """ (x1, y1) = self.ll2xy(*end1) (x2, y2) = self.ll2xy(*end2) return (numpy.degrees(numpy.arctan2(x2 - x1, y2 - y1)) + 180.) % 360. - 180.
[docs]class H2DRegLLGeometry(H2DRectangularGridGeometry): """ Handles the geometry for a Regular Lon/Lat Horizontal 2-Dimensions Field. """ _collector = ('geometry',) _footprint = dict( attr=dict( name=dict( values=set(['regular_lonlat'])) ) ) def _consistency_check(self): """ Check that the geometry is consistent. """ if self.name == 'regular_lonlat': grid_keys = ['center_lon', 'center_lat', 'X_resolution', 'Y_resolution'] if set(self.grid.keys()) != set(grid_keys): raise epygramError("grid attribute must consist in keys: " + \ str(grid_keys)) dimensions_keys = ['X', 'Y'] if set(self.dimensions.keys()) != set(dimensions_keys): raise epygramError("dimensions attribute must consist in keys: " + str(dimensions_keys))
[docs] def ij2xy(self, i, j, position=None): """ (*i, j*) being the coordinates in the 2D matrix of gridpoints, \n (*x, y*) being the coordinates in the projection. - *position*: position to return with respect to model cell. Defaults to position_on_grid['horizontal']. Note that origin of coordinates in projection is the center of the C+I domain. """ if isinstance(i, list) or isinstance(i, tuple): i = numpy.array(i) if isinstance(j, list) or isinstance(j, tuple): j = numpy.array(j) (oi, oj) = self._getoffset(position) Xpoints = self.dimensions['X'] Ypoints = self.dimensions['Y'] Xresolution = self.grid['X_resolution'].get('degrees') Yresolution = self.grid['Y_resolution'].get('degrees') Xorigin = self.grid['center_lon'].get('degrees') Yorigin = self.grid['center_lat'].get('degrees') # origin of coordinates is the center of domain i0 = float(Xpoints - 1) / 2. j0 = float(Ypoints - 1) / 2. x = Xorigin + (i - i0 + oi) * Xresolution y = Yorigin + (j - j0 + oj) * Yresolution return (x, y)
[docs] def xy2ij(self, x, y, position=None): """ (*x, y*) being the coordinates in the projection, \n (*i, j*) being the coordinates in the 2D matrix of gridpoints. - *position*: position represented by (x,y) with respect to model cell. Defaults to position_on_grid['horizontal']. Caution: (*i,j*) are float (the nearest grid point is the nearest integer). Note that origin of coordinates in projection is the center of the C+I domain. """ if isinstance(x, list) or isinstance(x, tuple): x = numpy.array(x) if isinstance(y, list) or isinstance(y, tuple): y = numpy.array(y) (oi, oj) = self._getoffset(position) Xpoints = self.dimensions['X'] Ypoints = self.dimensions['Y'] Xresolution = self.grid['X_resolution'].get('degrees') Yresolution = self.grid['Y_resolution'].get('degrees') Xorigin = self.grid['center_lon'].get('degrees') Yorigin = self.grid['center_lat'].get('degrees') # origin of coordinates is the center of domain i0 = float(Xpoints - 1) / 2. j0 = float(Ypoints - 1) / 2. i = i0 + (x - Xorigin) / Xresolution - oi j = j0 + (y - Yorigin) / Yresolution - oj return (i, j)
[docs] def ij2ll(self, i, j, position=None): """ (*i, j*) being the coordinates in the 2D matrix of gridpoints, \n (*lon, lat*) being the lon/lat coordinates in degrees. - *position*: lat lon position to return with respect to the model cell. Defaults to position_on_grid['horizontal']. """ return self.ij2xy(i, j, position)
[docs] def ll2ij(self, lon, lat, position=None): """ (*lon, lat*) being the lon/lat coordinates in degrees, \n (*i, j*) being the coordinates in the 2D matrix of gridpoints. - *position*: lat lon position with respect to the model cell. Defaults to position_on_grid['horizontal']. Caution: (*i,j*) are float. """ lon = degrees_nearest_mod(lon, self.grid['center_lon'].get('degrees')) return self.xy2ij(lon, lat, position)
[docs] def ll2xy(self, lon, lat): """ (*lon, lat*) being the lon/lat coordinates in degrees, \n (*x, y*) being the coordinates in the projection. Note that origin of coordinates in projection is the center of the C+I domain. """ if isinstance(lon, list) or isinstance(lon, tuple): lon = numpy.array(lon) if isinstance(lat, list) or isinstance(lat, tuple): lat = numpy.array(lat) return self.ij2xy(*self.ll2ij(lon, lat))
[docs] def xy2ll(self, x, y): """ (*x, y*) being the coordinates in the projection, \n (*lon, lat*) being the lon/lat coordinates in degrees. Note that origin of coordinates in projection is the center of the C+I domain. """ if isinstance(x, list) or isinstance(x, tuple): x = numpy.array(x) if isinstance(y, list) or isinstance(y, tuple): y = numpy.array(y) return self.ij2ll(*self.xy2ij(x, y))
[docs] def make_basemap(self, gisquality='i', subzone=None, specificproj=None, zoom=None): """ Returns a :class:`matplotlib.basemap.Basemap` object of the 'ad hoc' projection (if available). This is designed to avoid explicit handling of deep horizontal geometry attributes. Args: \n - *gisquality*: defines the quality of GIS contours, cf. Basemap doc. \n Possible values (by increasing quality): 'c', 'l', 'i', 'h', 'f'. - *subzone*: defines the LAM subzone to be included, in LAM case, among: 'C', 'CI'. - *specificproj*: enables to make basemap on the specified projection, among: 'kav7', 'cyl', 'ortho', 'nsperXXXX' (cf. Basemap doc). \n The XXXX of 'nsperXXXX' defines the satellite height in km. \n Overwritten by *zoom*. - *zoom*: specifies the lon/lat borders of the map, implying hereby a 'cyl' projection. Must be a dict(lonmin=, lonmax=, latmin=, latmax=).\n Overwrites *specificproj*. """ from mpl_toolkits.basemap import Basemap # corners (llcrnrlon, llcrnrlat) = self.ij2ll(*self.gimme_corners_ij(subzone)['ll']) (urcrnrlon, urcrnrlat) = self.ij2ll(*self.gimme_corners_ij(subzone)['ur']) # make basemap if zoom != None: # zoom case specificproj = 'cyl' llcrnrlon = zoom['lonmin'] llcrnrlat = zoom['latmin'] urcrnrlon = zoom['lonmax'] urcrnrlat = zoom['latmax'] if specificproj == None: # defaults if llcrnrlat <= -90.0 + config.epsilon or \ urcrnrlat >= 90.0 - config.epsilon: proj = 'cyl' else: proj = 'merc' b = Basemap(resolution=gisquality, projection=proj, llcrnrlon=llcrnrlon, llcrnrlat=llcrnrlat, urcrnrlon=urcrnrlon, urcrnrlat=urcrnrlat) else: # specificproj lon0 = self.grid['center_lon'].get('degrees') lat0 = self.grid['center_lat'].get('degrees') if specificproj == 'kav7': b = Basemap(resolution=gisquality, projection=specificproj, lon_0=lon0, llcrnrlon=llcrnrlon, llcrnrlat=llcrnrlat, urcrnrlon=urcrnrlon, urcrnrlat=urcrnrlat) elif specificproj == 'ortho': b = Basemap(resolution=gisquality, projection=specificproj, lon_0=lon0, lat_0=lat0) elif specificproj == 'cyl': b = Basemap(resolution=gisquality, projection=specificproj, llcrnrlon=llcrnrlon, llcrnrlat=llcrnrlat, urcrnrlon=urcrnrlon, urcrnrlat=urcrnrlat) elif specificproj == 'moll': b = Basemap(resolution=gisquality, projection=specificproj, lon_0=lon0, llcrnrlon=llcrnrlon, llcrnrlat=llcrnrlat, urcrnrlon=urcrnrlon, urcrnrlat=urcrnrlat) elif specificproj[0:5] == 'nsper': try: sat_height = float(specificproj[5:]) * 1000. except Exception: epylog.warning("syntax error with specificproj:" + \ specificproj + \ " => replaced by 'nsper3000'.") sat_height = 3000.*1000. b = Basemap(resolution=gisquality, projection=specificproj[0:5], lon_0=lon0, lat_0=lat0, satellite_height=sat_height) return b
[docs] def global_shift_center(self, longitude_shift): """ Shifts the center of the geometry by *longitude_shift* (in degrees). *longitude_shift* has to be a multiple of the grid's resolution in longitude. """ corners = self.gimme_corners_ll() if abs(abs(degrees_nearest_mod(corners['ur'][0] - corners['ul'][0], 0.)) - self.grid['X_resolution'].get('degrees')) \ < config.epsilon: if abs(longitude_shift % self.grid['X_resolution'].get('degrees')) > config.epsilon: raise epygramError("*longitude_shift* has to be a multiple of \ the grid's resolution in longitude.") self.grid['center_lon'] = Angle(self.grid['center_lon'].get('degrees') + longitude_shift, 'degrees') else: raise epygramError("unable to shift center if \ longitude_max - longitude_min != X_resolution.")
[docs] def linspace(self, end1, end2, num): """ Returns evenly spaced points over the specified interval. Points are lined up in the geometry. *end1* must be a tuple (lon, lat). *end2* must be a tuple (lon, lat). *num* is the number of points, including point1 and point2. """ if num < 2: raise epygramError("'num' must be at least 2.") (x1, y1) = self.ll2xy(*end1) (x2, y2) = self.ll2xy(*end2) xy_linspace = zip(numpy.linspace(x1, x2, num=num), numpy.linspace(y1, y2, num=num)) return [self.xy2ll(*xy) for xy in xy_linspace]
[docs] def distance(self, end1, end2): """ Computes the distance between two points along a straight line in the geometry. *end1* must be a tuple (lon, lat). *end2* must be a tuple (lon, lat). Warning: requires the :mod:`pyproj` module. """ import pyproj plast = end1 distance = 0 g = pyproj.Geod(ellps='sphere') for p in self.linspace(end1, end2, 1000)[1:]: distance += g.inv(plast[0], plast[1], *p)[2] plast = p return distance
[docs] def resolution_ll(self, lon, lat): """ Returns the local resolution at the nearest point of lon/lat. It's the distance between this point and its closest neighbour. *point* must be a tuple (lon, lat). """ return self.resolution_ij(*self.ll2ij(lon, lat))
[docs] def resolution_ij(self, i, j): """ Returns the distance to the nearest point of (i,j) point. (*i, j*) being the coordinates in the 2D matrix of gridpoints. """ (iint, jint) = (numpy.rint(i).astype('int'), numpy.rint(j).astype('int')) points_list = [(iint + oi, jint + oj) for oi in [-1, 0, 1] for oj in [-1, 0, 1] if (oi, oj) != (0, 0)] return numpy.array([self.distance(self.ij2ll(iint, jint), self.ij2ll(*p)) for p in points_list]).min()
[docs] def azimuth(self, end1, end2): """ Initial bearing from *end1* to *end2* points in geometry. *end1* must be a tuple (lon, lat). *end2* must be a tuple (lon, lat). """ (x1, y1) = self.ll2xy(*end1) (x2, y2) = self.ll2xy(*end2) return (numpy.degrees(numpy.arctan2(x2 - x1, y2 - y1)) + 180.) % 360. - 180.
[docs]class H2DProjectedGeometry(H2DRectangularGridGeometry): """ Handles the geometry for a Projected Horizontal 2-Dimensions Field. """ _collector = ('geometry',) _footprint = dict( attr=dict( name=dict( values=set(['lambert', 'mercator', 'polar_stereographic'])), projection=dict( type=FPDict, info="Handles projection information."), projtool=dict( optional=True, default=config.default_projtool, info="For projected geometries, to use pyproj or\ epygram.myproj."), geoid=dict( type=FPDict, optional=True, default=FPDict(config.default_geoid), info="For projected geometries only.") ) ) @property
[docs] def secant_projection(self): """ Is the projection secant to the sphere ? (or tangent)""" return 'secant_lat' in self.projection \ or 'secant_lat1' in self.projection
def __init__(self, *args, **kwargs): """Constructor. See its footprint for arguments.""" super(H2DProjectedGeometry, self).__init__(*args, **kwargs) if self.projtool == 'pyproj': import pyproj projtool = pyproj projdict = {'lambert':'lcc', 'mercator':'merc', 'polar_stereographic':'stere'} elif self.projtool == 'myproj': from epygram import myproj projtool = myproj projdict = {'lambert':'lambert', 'mercator':'mercator', 'polar_stereographic':'polar_stereographic'} proj = projdict[self.name] # build proj # CAUTION: x0, y0 are deeply linked with ij2xy and xy2ij methods: # the origin of the grid is the center of the domain ! if self.name == 'lambert': if self.secant_projection: lat_1 = self.projection['secant_lat1'].get('degrees') lat_2 = self.projection['secant_lat2'].get('degrees') m1 = math.cos(math.radians(lat_1)) m2 = math.cos(math.radians(lat_2)) t1 = math.tan(math.pi / 4. - math.radians(lat_1) / 2.) t2 = math.tan(math.pi / 4. - math.radians(lat_2) / 2.) self._K = (math.log(m1) - math.log(m2)) / \ (math.log(t1) - math.log(t2)) else: lat_1 = self.projection['reference_lat'].get('degrees') lat_2 = self.projection['reference_lat'].get('degrees') self._K = abs(self.projection['reference_lat'].get('cos_sin')[1]) p = projtool.Proj(proj=proj, lon_0=self.projection['reference_lon'].get('degrees'), lat_1=lat_1, lat_2=lat_2, **self.geoid) x0, y0 = p(self.projection['center_lon'].get('degrees'), self.projection['center_lat'].get('degrees')) self._proj = projtool.Proj(proj=proj, lon_0=self.projection['reference_lon'].get('degrees'), lat_1=lat_1, lat_2=lat_2, x_0=-x0, y_0=-y0, **self.geoid) elif self.name == 'mercator': if self.secant_projection: lat_ts = self.projection['secant_lat'].get('degrees') else: lat_ts = 0. p = projtool.Proj(proj=proj, lon_0=self.projection['reference_lon'].get('degrees'), lat_ts=lat_ts, **self.geoid) x0, y0 = p(self.projection['center_lon'].get('degrees'), self.projection['center_lat'].get('degrees')) self._proj = projtool.Proj(proj=proj, lon_0=self.projection['reference_lon'].get('degrees'), lat_ts=lat_ts, x_0=-x0, y_0=-y0, **self.geoid) elif self.name == 'polar_stereographic': lat_0 = self.projection['reference_lat'].get('degrees') if self.secant_projection: lat_ts = self.projection['secant_lat'].get('degrees') else: lat_ts = self.projection['reference_lat'].get('degrees') p = projtool.Proj(proj=proj, lon_0=self.projection['reference_lon'].get('degrees'), lat_0=lat_0, lat_ts=lat_ts, **self.geoid) x0, y0 = p(self.projection['center_lon'].get('degrees'), self.projection['center_lat'].get('degrees')) self._proj = projtool.Proj(proj=proj, lon_0=self.projection['reference_lon'].get('degrees'), lat_0=lat_0, lat_ts=lat_ts, x_0=-x0, y_0=-y0, **self.geoid) else: raise NotImplementedError("projection: " + self.name) def _consistency_check(self): """ Check that the geometry is consistent. """ if self.name == 'lambert': sets_of_keys = (['reference_lon', 'reference_lat', 'center_lon', 'center_lat'], ['reference_lon', 'secant_lat1', 'secant_lat2', 'center_lon', 'center_lat']) elif self.name in ('polar_stereographic', 'mercator'): sets_of_keys = (['reference_lon', 'reference_lat', 'center_lon', 'center_lat'], ['reference_lon', 'reference_lat', 'secant_lat', 'center_lon', 'center_lat']) if set(self.projection.keys()) != set(sets_of_keys[0]) \ and set(self.projection.keys()) != set(sets_of_keys[1]): print self.projection.keys() raise epygramError("attributes for projection " + self.name + \ " must consist in keys: " + \ str(sets_of_keys[0]) + " or " + \ str(sets_of_keys[1])) grid_keys = ['LAMzone', 'X_resolution', 'Y_resolution'] if set(self.grid.keys()) != set(grid_keys): raise epygramError("grid attribute must consist in keys: " + \ str(grid_keys)) LAMzone_values = [None, 'CI', 'CIE'] if self.grid['LAMzone'] not in LAMzone_values: raise epygramError("grid['LAMzone'] must be among " + \ str(LAMzone_values)) dimensions_keys = ['X', 'Y'] if self.grid['LAMzone'] in ('CI', 'CIE'): dimensions_keys.extend(['X_Iwidth', 'Y_Iwidth', 'X_Czone', 'Y_Czone', 'X_CIzone', 'Y_CIzone']) if self.grid['LAMzone'] == 'CIE': dimensions_keys.extend(['X_CIoffset', 'Y_CIoffset']) if set(self.dimensions.keys()) != set(dimensions_keys): raise epygramError("dimensions attribute must consist in keys: " + \ str(dimensions_keys))
[docs] def ij2xy(self, i, j, position=None): """ (*i, j*) being the coordinates in the 2D matrix of CIE gridpoints, \n (*x, y*) being the coordinates in the projection. - *position*: position to return with respect to model cell. Defaults to position_on_grid['horizontal']. Note that origin of coordinates in projection is the center of the C+I domain. """ if isinstance(i, list) or isinstance(i, tuple): i = numpy.array(i) if isinstance(j, list) or isinstance(j, tuple): j = numpy.array(j) Xresolution = self.grid['X_resolution'] Yresolution = self.grid['Y_resolution'] if self.grid['LAMzone'] == None: Xpoints = self.dimensions['X'] Ypoints = self.dimensions['Y'] else: Xpoints = self.dimensions['X_CIzone'] Ypoints = self.dimensions['Y_CIzone'] Xorigin = 0. Yorigin = 0. # origin of coordinates is the center of CI domain i0 = self.dimensions['X_CIoffset'] + float(Xpoints - 1) / 2. j0 = self.dimensions['Y_CIoffset'] + float(Ypoints - 1) / 2. (oi, oj) = self._getoffset(position) x = Xorigin + (i - i0 + oi) * Xresolution y = Yorigin + (j - j0 + oj) * Yresolution return (x, y)
[docs] def xy2ij(self, x, y, position=None): """ (*x, y*) being the coordinates in the projection, \n (*i, j*) being the coordinates in the 2D matrix of CIE gridpoints. - *position*: position represented by (x,y) with respect to model cell. Defaults to position_on_grid['horizontal']. Caution: (*i,j*) are float (the nearest grid point is the nearest integer). Note that origin of coordinates in projection is the center of the C+I domain. """ if isinstance(x, list) or isinstance(x, tuple): x = numpy.array(x) if isinstance(y, list) or isinstance(y, tuple): y = numpy.array(y) Xresolution = self.grid['X_resolution'] Yresolution = self.grid['Y_resolution'] if self.grid['LAMzone'] == None: Xpoints = self.dimensions['X'] Ypoints = self.dimensions['Y'] else: Xpoints = self.dimensions['X_CIzone'] Ypoints = self.dimensions['Y_CIzone'] Xorigin = 0. Yorigin = 0. # origin of coordinates is the center of CI domain i0 = self.dimensions['X_CIoffset'] + float(Xpoints - 1) / 2. j0 = self.dimensions['Y_CIoffset'] + float(Ypoints - 1) / 2. (oi, oj) = self._getoffset(position) i = i0 + (x - Xorigin) / Xresolution - oi j = j0 + (y - Yorigin) / Yresolution - oj return (i, j)
[docs] def ll2xy(self, lon, lat): """ (*lon, lat*) being the lon/lat coordinates in degrees, \n (*x, y*) being the coordinates in the projection. Note that origin of coordinates in projection is the center of the C+I domain. """ if isinstance(lon, list) or isinstance(lon, tuple): lon = numpy.array(lon) if isinstance(lat, list) or isinstance(lat, tuple): lat = numpy.array(lat) lon = degrees_nearest_mod(lon, self.projection['reference_lon'].get('degrees')) xy = self._proj(lon, lat) return xy
[docs] def xy2ll(self, x, y): """ (*x, y*) being the coordinates in the projection, \n (*lon, lat*) being the lon/lat coordinates in degrees. Note that origin of coordinates in projection is the center of the C+I domain. """ if isinstance(x, list) or isinstance(x, tuple): x = numpy.array(x) if isinstance(y, list) or isinstance(y, tuple): y = numpy.array(y) ll = self._proj(x, y, inverse=True) return ll
[docs] def ij2ll(self, i, j, position=None): """ (*i, j*) being the coordinates in the 2D matrix of CIE gridpoints, \n (*lon, lat*) being the lon/lat coordinates in degrees. - *position*: lat lon position to return with respect to the model cell. Defaults to position_on_grid['horizontal']. """ return self.xy2ll(*self.ij2xy(i, j, position))
[docs] def ll2ij(self, lon, lat, position=None): """ (*lon, lat*) being the lon/lat coordinates in degrees, \n (*i, j*) being the coordinates in the 2D matrix of CIE gridpoints. - *position*: lat lon position with respect to the model cell. Defaults to position_on_grid['horizontal']. Caution: (*i,j*) are float. """ return self.xy2ij(*self.ll2xy(lon, lat), position=position)
[docs] def map_factor_lat(self, lat): """Returns the map factor at the given latitude(s) *lat* in degrees.""" lat = numpy.radians(lat) if self.name == 'mercator': if self.secant_projection: lat_0 = self.projection['secant_lat'].get('radians') else: lat_0 = 0. m = numpy.cos(lat_0) / numpy.cos(lat) elif self.name == 'polar_stereographic': if self.secant_projection: lat_0 = self.projection['secant_lat'].get('radians') else: lat_0 = self.projection['reference_lat'].get('radians') m = (1. + numpy.copysign(1., lat_0) * numpy.sin(lat_0)) / \ (1. + numpy.copysign(1., lat_0) * numpy.sin(lat)) elif self.name == 'lambert': if self.secant_projection: lat_0 = self.projection['secant_lat2'].get('radians') else: lat_0 = self.projection['reference_lat'].get('radians') m = (numpy.cos(lat_0) / numpy.cos(lat)) ** (1. - self._K) \ * ((1. + numpy.copysign(1., lat_0) * numpy.sin(lat_0)) \ / (1. + numpy.copysign(1., lat_0) * numpy.sin(lat))) ** self._K else: raise epygramError("projection " + self.name + " not implemented.") return m
[docs] def map_factor_field(self, position=None): """ Returns a new field whose data is the map factor over the field. - *position*: grid position with respect to the model cell. Defaults to position_on_grid['horizontal']. """ f = fpx.field(geometry=self, fid=FPDict({'geometry':'Map Factor'})) lats = self.get_lonlat_grid(position=position)[1] data = self.map_factor_lat(lats) f.setdata(data) return f
[docs] def make_basemap(self, gisquality='i', subzone=None, specificproj=None, zoom=None): """ Returns a :class:`matplotlib.basemap.Basemap` object of the 'ad hoc' projection (if available). This is designed to avoid explicit handling of deep horizontal geometry attributes. Args: \n - *gisquality*: defines the quality of GIS contours, cf. Basemap doc. \n Possible values (by increasing quality): 'c', 'l', 'i', 'h', 'f'. - *subzone*: defines the LAM subzone to be included, in LAM case, among: 'C', 'CI'. - *specificproj*: enables to make basemap on the specified projection, among: 'kav7', 'cyl', 'ortho', 'nsperXXXX' (cf. Basemap doc). \n The XXXX of 'nsperXXXX' defines the satellite height in km. \n Overwritten by *zoom*. - *zoom*: specifies the lon/lat borders of the map, implying hereby a 'cyl' projection. Must be a dict(lonmin=, lonmax=, latmin=, latmax=). \n Overwrites *specificproj*. """ from mpl_toolkits.basemap import Basemap if self.secant_projection: lat = 'secant_lat' else: lat = 'reference_lat' # corners (llcrnrlon, llcrnrlat) = self.ij2ll(*self.gimme_corners_ij(subzone)['ll']) (urcrnrlon, urcrnrlat) = self.ij2ll(*self.gimme_corners_ij(subzone)['ur']) # make basemap if zoom != None: # zoom case specificproj = 'cyl' llcrnrlon = zoom['lonmin'] llcrnrlat = zoom['latmin'] urcrnrlon = zoom['lonmax'] urcrnrlat = zoom['latmax'] if specificproj == None: # defaults if self.name == 'lambert': if self.secant_projection: lat_0 = (self.projection['secant_lat1'].get('degrees') \ + self.projection['secant_lat2'].get('degrees')) / 2. b = Basemap(resolution=gisquality, projection='lcc', lat_1=self.projection['secant_lat1'].get('degrees'), lat_2=self.projection['secant_lat2'].get('degrees'), lat_0=lat_0, lon_0=self.projection['reference_lon'].get('degrees'), llcrnrlon=llcrnrlon, llcrnrlat=llcrnrlat, urcrnrlon=urcrnrlon, urcrnrlat=urcrnrlat) else: b = Basemap(resolution=gisquality, projection='lcc', lat_0=self.projection['reference_lat'].get('degrees'), lon_0=self.projection['reference_lon'].get('degrees'), llcrnrlon=llcrnrlon, llcrnrlat=llcrnrlat, urcrnrlon=urcrnrlon, urcrnrlat=urcrnrlat) elif self.name == 'mercator': b = Basemap(resolution=gisquality, projection='merc', lat_ts=self.projection[lat].get('degrees'), lon_0=self.projection['center_lon'].get('degrees'), llcrnrlon=llcrnrlon, llcrnrlat=llcrnrlat, urcrnrlon=urcrnrlon, urcrnrlat=urcrnrlat) elif self.name == 'polar_stereographic': b = Basemap(resolution=gisquality, projection='stere', lat_ts=self.projection[lat].get('degrees'), lat_0=numpy.copysign(90., self.projection[lat].get('degrees')), lon_0=self.projection['reference_lon'].get('degrees'), llcrnrlon=llcrnrlon, llcrnrlat=llcrnrlat, urcrnrlon=urcrnrlon, urcrnrlat=urcrnrlat) else: # specificproj lon0 = self.projection['center_lon'].get('degrees') lat0 = self.projection['center_lat'].get('degrees') if specificproj == 'kav7': b = Basemap(resolution=gisquality, projection=specificproj, lon_0=lon0, llcrnrlon=llcrnrlon, llcrnrlat=llcrnrlat, urcrnrlon=urcrnrlon, urcrnrlat=urcrnrlat) elif specificproj == 'ortho': b = Basemap(resolution=gisquality, projection=specificproj, lon_0=lon0, lat_0=lat0) elif specificproj == 'cyl': b = Basemap(resolution=gisquality, projection=specificproj, llcrnrlon=llcrnrlon, llcrnrlat=llcrnrlat, urcrnrlon=urcrnrlon, urcrnrlat=urcrnrlat) elif specificproj == 'moll': b = Basemap(resolution=gisquality, projection=specificproj, lon_0=lon0, llcrnrlon=llcrnrlon, llcrnrlat=llcrnrlat, urcrnrlon=urcrnrlon, urcrnrlat=urcrnrlat) elif specificproj[0:5] == 'nsper': try: sat_height = float(specificproj[5:]) * 1000. except ValueError: epylog.warning("syntax error with specificproj:" + \ specificproj + \ " => replaced by 'nsper3000'.") sat_height = 3000.*1000. b = Basemap(resolution=gisquality, projection=specificproj[0:5], lon_0=lon0, lat_0=lat0, satellite_height=sat_height) return b
[docs] def linspace(self, end1, end2, num): """ Returns evenly spaced points over the specified interval. Points are lined up in the geometry. *end1* must be a tuple (lon, lat). *end2* must be a tuple (lon, lat). *num* is the number of points, including point1 and point2. """ if num < 2: raise epygramError("'num' must be at least 2.") (x1, y1) = self.ll2xy(*end1) (x2, y2) = self.ll2xy(*end2) xy_linspace = zip(numpy.linspace(x1, x2, num=num), \ numpy.linspace(y1, y2, num=num)) return [self.xy2ll(*xy) for xy in xy_linspace]
[docs] def distance(self, end1, end2): """ Computes the distance between two points along a straight line in the geometry. *end1* must be a tuple (lon, lat). *end2* must be a tuple (lon, lat). """ (x1, y1) = self.ll2xy(*end1) (x2, y2) = self.ll2xy(*end2) #distance on the projected surface distance = numpy.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2) # map factor computed as the mean over 1000 points along the line # between point1 and point2 mean_map_factor = numpy.array([self.map_factor_lat(lat) for (lon, lat) in self.linspace(end1, end2, 1000)]).mean() return distance / mean_map_factor
[docs] def resolution_ll(self, lon, lat): """ Returns the local resolution at the nearest point of lon/lat. It's the distance between this point and its closest neighbour. *point* must be a tuple (lon, lat). """ return self.resolution_ij(*self.ll2ij(lon, lat))
[docs] def resolution_ij(self, i, j): """ Returns the distance to the nearest point of (i,j) point. (*i, j*) being the coordinates in the 2D matrix of gridpoints. """ (iint, jint) = (numpy.rint(i).astype('int'), numpy.rint(j).astype('int')) points_list = [(iint + oi, jint + oj) for oi in [-1, 0, 1] for oj in [-1, 0, 1] if (oi, oj) != (0, 0)] return numpy.array([self.distance(self.ij2ll(iint, jint), self.ij2ll(*p)) for p in points_list]).min()
[docs] def azimuth(self, end1, end2): """ Initial bearing from *end1* to *end2* points in geometry. *end1* must be a tuple (lon, lat). *end2* must be a tuple (lon, lat). """ (x1, y1) = self.ll2xy(*end1) (x2, y2) = self.ll2xy(*end2) return (numpy.degrees(numpy.arctan2(x2 - x1, y2 - y1)) + 180.) % 360. - 180.
[docs]class H2DGaussGeometry(H2DGeometry): """ Handles the geometry for a Global Gauss grid Horizontal 2-Dimensions Field. """ _collector = ('geometry',) _footprint = dict( attr=dict( name=dict( values=set(['rotated_reduced_gauss', 'reduced_gauss'])), ) ) def _consistency_check(self): """Check that the geometry is consistent.""" grid_keys = ['dilatation_coef', 'latitudes'] if self.name == 'rotated_reduced_gauss': grid_keys.extend(['pole_lon', 'pole_lat']) if set(self.grid.keys()) != set(grid_keys): raise epygramError("grid attribute must consist in keys: " + \ str(grid_keys)) dimensions_keys = ['max_lon_number', 'lat_number', 'lon_number_by_lat'] if set(self.dimensions.keys()) != set(dimensions_keys): raise epygramError("dimensions attribute must consist in keys: " + \ str(dimensions_keys))
[docs] def ij2ll(self, i, j, position=None): """ (*i, j*) being the coordinates in the 2D matrix of gridpoints, \n (*lon, lat*) being the lon/lat coordinates in degrees. - *position*: lat lon position to return with respect to the model cell. Defaults to position_on_grid['horizontal']. Computation adapted from arpifs/transform/tracare.F90 """ if self._getoffset(position) != (0., 0.): raise NotImplementedError("horizontal staggered grids for reduced\ gauss grid are not implemented.") if self.name == 'rotated_reduced_gauss': KTYP = 2 elif self.name == 'reduced_gauss': KTYP = 1 if isinstance(i, list) or isinstance(i, tuple) or\ isinstance(i, numpy.ndarray): sinlats = [self.grid['latitudes'][n].get('cos_sin')[1] for n in j] lons = [numpy.pi * 2 * i[n] / self.dimensions['lon_number_by_lat'][j[n]] for n in range(len(j))] PSLAC = numpy.array(sinlats) PCLOC = numpy.cos(lons) PSLOC = numpy.sin(lons) else: sinlat = self.grid['latitudes'][j].get('cos_sin')[1] lon = numpy.pi * 2 * i / self.dimensions['lon_number_by_lat'][j] PSLAC = sinlat PCLOC = numpy.cos(lon) PSLOC = numpy.sin(lon) PFACDI = self.grid['dilatation_coef'] ZDBLC = 2.0 * PFACDI ZC2P1 = PFACDI * PFACDI + 1.0 ZC2M1 = PFACDI * PFACDI - 1.0 ZEPS = config.epsilon if KTYP == 2: # Complete ARPEGE (PCLAP, PSLAP) = self.grid['pole_lat'].get('cos_sin') (PCLOP, PSLOP) = self.grid['pole_lon'].get('cos_sin') ZCLAC = numpy.sqrt(1. - PSLAC * PSLAC) ZA = 1. / (ZC2P1 + ZC2M1 * PSLAC) ZB = ZC2P1 * PSLAC + ZC2M1 ZC = ZDBLC * PCLAP * ZCLAC * PCLOC + ZB * PSLAP PSLAR = ZC * ZA ZD = ZA / numpy.sqrt(numpy.maximum(ZEPS, 1. - PSLAR * PSLAR)) ZE = ZB * PCLAP * PCLOP - ZDBLC * ZCLAC * \ (PSLAP * PCLOC * PCLOP - PSLOP * PSLOC) ZF = ZB * PCLAP * PSLOP - ZDBLC * ZCLAC * \ (PSLAP * PCLOC * PSLOP + PCLOP * PSLOC) PCLOR = ZE * ZD PSLOR = ZF * ZD PSLAR = numpy.maximum(-1., numpy.minimum(1., PSLAR)) PCLOR = numpy.maximum(-1., numpy.minimum(1., PCLOR)) PSLOR = numpy.maximum(-1., numpy.minimum(1., PSLOR)) else: # Schmidt PSLOR = PSLOC PCLOR = PCLOC PSLAR = (ZC2P1 * PSLAC + ZC2M1) / (ZC2P1 + ZC2M1 * PSLAC) PSLAR = numpy.maximum(-1., numpy.minimum(1., PSLAR)) lon = numpy.degrees(numpy.arctan2(PSLOR, PCLOR)) lat = numpy.degrees(numpy.arcsin(PSLAR)) return (lon, lat)
[docs] def ll2ij(self, lon, lat, position=None): """ (*lon, lat*) being the lon/lat coordinates in degrees, \n (*i, j*) being the coordinates in the 2D matrix of gridpoints. - *position*: lat lon position with respect to the model cell. Defaults to position_on_grid['horizontal']. """ if isinstance(lon, list) or isinstance(lon, tuple): lon = numpy.array(lon) if isinstance(lat, list) or isinstance(lat, tuple): lat = numpy.array(lat) return self.nearest_points(lon, lat, 'nearest', position)
[docs] def get_lonlat_grid(self, position=None, **kwargs): """ Returns a tuple of two tables containing one the longitude of each point, the other the latitude, with 2D shape. Shape of 2D data in Gauss grids: \n - grid[0, 0:Nj] is first (Northern) band of latitude, masked after Nj = number of longitudes for latitude j \n - grid[-1, 0:Nj] is last (Southern) band of latitude (idem). Args:\n - *position*: position of lonlat grid with respect to the model cell. Defaults to position_on_grid['horizontal']. """ # !!! **kwargs enables the method to receive arguments specific to # other geometries, useless here ! Do not remove. if hasattr(self, '_buffered_gauss_grid'): lons = self._buffered_gauss_grid[0] lats = self._buffered_gauss_grid[1] elif hasattr(self, '_bound_gauss_grid'): (lons, lats) = self._bound_gauss_grid() else: Jmax = self.dimensions['lat_number'] Imax = self.dimensions['lon_number_by_lat'] igrid = [] jgrid = [] for j in range(Jmax): for i in range(Imax[j]): igrid.append(i) jgrid.append(j) igrid = numpy.array(igrid) jgrid = numpy.array(jgrid) (lons, lats) = self.ij2ll(igrid, jgrid, position) lons = self.reshape_data(lons) lats = self.reshape_data(lats) if config.FA_buffered_gauss_grid: self._buffered_gauss_grid = (lons, lats) return (lons, lats)
[docs] def reshape_data(self, data): """ Returns a 2D data reshaped from 1D, according to geometry. - *data*: the 1D data, of dimension concording with geometry. """ rdata = numpy.ma.masked_all((self.dimensions['lat_number'], self.dimensions['max_lon_number'])) ind_end = 0 for j in range(self.dimensions['lat_number']): ind_begin = ind_end ind_end = ind_begin + self.dimensions['lon_number_by_lat'][j] rdata[j, 0:self.dimensions['lon_number_by_lat'][j]] = data[ind_begin:ind_end] return rdata
[docs] def make_basemap(self, gisquality='i', specificproj=None, zoom=None, **kwargs): """ Returns a :class:`matplotlib.basemap.Basemap` object of the 'ad hoc' projection (if available). This is designed to avoid explicit handling of deep horizontal geometry attributes. Args: \n - *gisquality*: defines the quality of GIS contours, cf. Basemap doc. \n Possible values (by increasing quality): 'c', 'l', 'i', 'h', 'f'. - *specificproj*: enables to make basemap on the specified projection, among: 'kav7', 'cyl', 'ortho', 'nsperXXXX' (cf. Basemap doc). \n The XXXX of 'nsperXXXX' defines the satellite height in km. \n Overwritten by *zoom*. - *zoom*: specifies the lon/lat borders of the map, implying hereby a 'cyl' projection. Must be a dict(lonmin=, lonmax=, latmin=, latmax=). \n Overwrites *specificproj*. """ # !!! **kwargs enables the method to receive arguments specific to # other geometries, useless here ! Do not remove. from mpl_toolkits.basemap import Basemap gisquality = 'l' # forced for Gauss, for time-consumption reasons... # corners llcrnrlon = -180 llcrnrlat = -90 urcrnrlon = 180 urcrnrlat = 90 # make basemap if zoom != None: # zoom case specificproj = 'cyl' llcrnrlon = zoom['lonmin'] llcrnrlat = zoom['latmin'] urcrnrlon = zoom['lonmax'] urcrnrlat = zoom['latmax'] if specificproj == None: # defaults if self.name == 'rotated_reduced_gauss': lon0 = self.grid['pole_lon'].get('degrees') elif self.name == 'reduced_gauss': lon0 = 0.0 b = Basemap(resolution=gisquality, projection='moll', lon_0=lon0) else: # specificproj if self.name == 'rotated_reduced_gauss': lon0 = self.grid['pole_lon'].get('degrees') lat0 = self.grid['pole_lat'].get('degrees') elif self.name == 'reduced_gauss': lon0 = 0.0 lat0 = 0.0 if specificproj == 'kav7': b = Basemap(resolution=gisquality, projection=specificproj, lon_0=lon0, llcrnrlon=llcrnrlon, llcrnrlat=llcrnrlat, urcrnrlon=urcrnrlon, urcrnrlat=urcrnrlat) elif specificproj == 'ortho': b = Basemap(resolution=gisquality, projection=specificproj, lon_0=lon0, lat_0=lat0) elif specificproj == 'cyl': b = Basemap(resolution=gisquality, projection=specificproj, llcrnrlon=llcrnrlon, llcrnrlat=llcrnrlat, urcrnrlon=urcrnrlon, urcrnrlat=urcrnrlat) elif specificproj == 'moll': b = Basemap(resolution=gisquality, projection=specificproj, lon_0=lon0, llcrnrlon=llcrnrlon, llcrnrlat=llcrnrlat, urcrnrlon=urcrnrlon, urcrnrlat=urcrnrlat) elif specificproj[0:5] == 'nsper': try: sat_height = float(specificproj[5:]) * 1000. except ValueError: epylog.warning("syntax error with specificproj:" + \ specificproj + \ " => replaced by 'nsper3000'.") sat_height = 3000.*1000. b = Basemap(resolution=gisquality, projection=specificproj[0:5], lon_0=lon0, lat_0=lat0, satellite_height=sat_height) return b
[docs] def distance(self, end1, end2): """ Computes the distance between two points along a Great Circle. *end1* must be a tuple (lon, lat). *end2* must be a tuple (lon, lat). Warning: requires the :mod:`pyproj` module. """ import pyproj g = pyproj.Geod(ellps='sphere') distance = g.inv(end1[0], end1[1], end2[0], end2[1])[2] return distance
[docs] def resolution_ll(self, lon, lat): """ Returns the local resolution at the nearest point of lon/lat. It's the distance between this point and its closest neighbour. *point* must be a tuple (lon, lat). """ return self.resolution_ij(*self.ll2ij(lon, lat))
[docs] def resolution_ij(self, i, j): """ Returns the distance to the nearest point of (i,j) point. (*i, j*) being the coordinates in the 2D matrix of gridpoints. """ (iint, jint) = (numpy.rint(i).astype('int'), numpy.rint(j).astype('int')) points_list = [] for oi in [-1, 0, 1]: for oj in [-1, 0, 1]: oj = 0 if (oi, oj) != (0, 0): pi = iint + oi pj = jint + oj if pj < self.dimensions['lat_number'] and pj >= 0: pi = pi % self.dimensions['lon_number_by_lat'][pj] points_list.append((pi, pj)) return numpy.array([self.distance(self.ij2ll(iint, jint), self.ij2ll(*p)) for p in points_list]).min()
[docs] def linspace(self, end1, end2, num): """ Returns evenly spaced points over the specified interval. Points are lined up along a Great Circle. *end1* must be a tuple (lon, lat). *end2* must be a tuple (lon, lat). *num* is the number of points, including point1 and point2. Warning: requires the :mod:`pyproj` module. """ import pyproj if num < 2: raise epygramError("'num' must be at least 2.") g = pyproj.Geod(ellps='sphere') transect = g.npts(end1[0], end1[1], end2[0], end2[1], num - 2) transect.insert(0, end1) transect.append(end2) return transect
[docs] def azimuth(self, end1, end2): """ Initial bearing from *end1* to *end2* points following a Great Circle. *end1* must be a tuple (lon, lat). *end2* must be a tuple (lon, lat). Warning: requires the :mod:`pyproj` module. """ import pyproj g = pyproj.Geod(ellps='sphere') return g.inv(end1[0], end1[1], end2[0], end2[1])[0]
[docs] def point_is_inside_domain_ll(self, lon, lat, margin=-0.1, position=None): """ Returns True if the point(s) of lon/lat coordinates is(are) inside the field. Args: \n - *lon*: longitude of point(s) in degrees. - *lat*: latitude of point(s) in degrees. - *margin*: considers the point inside if at least 'margin' points far from the border. The -0.1 default is a safety for precision errors. - *position*: position of the grid with respect to the model cell. Defaults to position_on_grid['horizontal']. """ try: N = len(lon) except Exception: N = 1 if N == 1: inside = True else: inside = [True for n in range(N)] return inside
[docs] def point_is_inside_domain_ij(self, i, j, margin=-0.1): """ Returns True if the point(s) of lon/lat coordinates is(are) inside the field. Args: \n - *lon*: longitude of point(s) in degrees. - *lat*: latitude of point(s) in degrees. - *margin*: considers the point inside if at least 'margin' points far from the border. The -0.1 default is a safety for precision errors. """ try: N = len(i) except Exception: N = 1 if N == 1: inside = True if j >= self.dimensions['lat_number'] or j < 0: inside = False if i >= self.dimensions['lon_number_by_lat'][j] or i < 0: inside = False else: inside = [self.point_is_inside_domain_ij(i[n], j[n], margin=margin) for n in range(N)] return inside
[docs] def nearest_points(self, lon, lat, interpolation, position=None): """ Returns a list of the (i, j) position of the points needed to perform an interpolation. Args: \n - *lon*: longitude of point in degrees. - *lat*: latitude of point in degrees. - *interpolation* can be: - nearest to get the nearest point only - linear to get the 2*2 points bordering the (*lon*, *lat*) position - cubic to get the 4*4 points bordering the (*lon*, *lat*) position - *position*: position in the model cell of the lat lon position. Defaults to position_on_grid['horizontal']. """ def rotateStretch(lon, lat): """ Internal method used to transfor true lon/lat into rotated and stretched lon/lat (*lon, lat*) being the lon/lat coordinates in degrees Computation adapted from arpifs/transform/trareca.F90. """ if self.name == 'rotated_reduced_gauss': KTYP = 2 elif self.name == 'reduced_gauss': KTYP = 1 if isinstance(lon, list) or isinstance(lon, tuple): lon = numpy.array(lon) if isinstance(lat, list) or isinstance(lat, tuple): lat = numpy.array(lat) lon = degrees_nearest_mod(lon, self.grid['pole_lon'].get('degrees')) PSLAR = numpy.sin(numpy.radians(lat)) PSLOR = numpy.sin(numpy.radians(lon)) PCLOR = numpy.cos(numpy.radians(lon)) PFACDI = self.grid['dilatation_coef'] ZC2P1 = PFACDI * PFACDI + 1. ZC2M1 = PFACDI * PFACDI - 1. ZCRAP = -ZC2M1 / ZC2P1 ZEPS = config.epsilon if KTYP == 2: # Complete ARPEGE (PCLAP, PSLAP) = self.grid['pole_lat'].get('cos_sin') (PCLOP, PSLOP) = self.grid['pole_lon'].get('cos_sin') ZCLAR = numpy.sqrt(1. - PSLAR * PSLAR) ZA = PSLAP * PSLAR + PCLAP * ZCLAR * \ (PCLOP * PCLOR + PSLOP * PSLOR) PSLAC = (ZCRAP + ZA) / (1. + ZA * ZCRAP) ZB = 1. / numpy.sqrt(numpy.maximum(ZEPS, 1. - ZA * ZA)) PCLOC = (PCLAP * PSLAR - PSLAP * ZCLAR * \ (PCLOP * PCLOR + PSLOP * PSLOR)) * ZB PSLOC = ZCLAR * (PSLOP * PCLOR - PCLOP * PSLOR) * ZB PSLAC = numpy.maximum(-1., numpy.minimum(1., PSLAC)) PCLOC = numpy.maximum(-1., numpy.minimum(1., PCLOC)) PSLOC = numpy.maximum(-1., numpy.minimum(1., PSLOC)) if isinstance(PSLOC, numpy.ndarray): for k in range(0, len(PSLOC)): if PCLOC[k] * PCLOC[k] + PSLOC[k] * PSLOC[k] < 0.5: PSLOC[k] = 1. PCLOC[k] = 0. else: if PCLOC * PCLOC + PSLOC * PSLOC < 0.5: PSLOC = 1. PCLOC = 0. else: # Schmidt PSLOC = PSLOR PCLOC = PCLOR PSLAC = (ZC2P1 * PSLAR - ZC2M1) / (ZC2P1 - ZC2M1 * PSLAR) PSLAC = numpy.maximum(-1., numpy.minimum(1., PSLAC)) lat = numpy.arcsin(PSLAC) lon = numpy.arctan2(PSLOC, PCLOC) % (numpy.pi * 2) return (lon, lat) def nearest_lats(latrs, num): """ Internal method used to find the nearest latitude circles. *latrs* is the rotated and streched latitude. *num* is: - 2 to find the two surrounding latitude cicles, - 3 to find the nearest latitude circles plus the one above and the one under. - 4 to find the two surrounding latitude circles plus the preceding one and the following one. """ if not numpy.all(numpy.array(latitudes[1:]) <= numpy.array(latitudes[:-1])): raise ValueError('latitudes must be in descending order') distmin = None nearest = None for n in range(0, len(latitudes)): dist = latrs - latitudes[n] if distmin == None or abs(dist) < abs(distmin): distmin = dist nearest = n if num == 2: lats = [nearest, (nearest - numpy.copysign(1, distmin)).astype('int')] elif num == 3: lats = [nearest - 1, nearest, nearest + 1] elif num == 4: lats = [(nearest + numpy.copysign(1, distmin)).astype('int'), nearest, (nearest - numpy.copysign(1, distmin)).astype('int'), (nearest - numpy.copysign(2, distmin)).astype('int')] return lats def nearest_lons(lonrs, latnum, num): """ Internal method used to find the nearest points on a latitude circle. Args:\n - *lonrs* is the rotated and streched longitude. - *num* is: - 1 to find the nearest point, - 2 to find the two surrounding points, - 4 to find the two surrounding points plus the preceding one and the following one. - *latnum*: if -1 (resp. -2), we search for the opposite longitude on the first (resp. second) latitude circle. The same is true for the other pole. """ if latnum == -2: j = 1 lonnummax = self.dimensions['lon_number_by_lat'][j] i = ((lonrs - numpy.pi) % (numpy.pi * 2)) * \ lonnummax / (numpy.pi * 2) elif latnum == -1: j = 0 lonnummax = self.dimensions['lon_number_by_lat'][j] i = ((lonrs - numpy.pi) % (numpy.pi * 2)) * \ lonnummax / (numpy.pi * 2) elif latnum == len(latitudes): j = latnum - 1 lonnummax = self.dimensions['lon_number_by_lat'][j] i = ((lonrs - numpy.pi) % (numpy.pi * 2)) * \ lonnummax / (numpy.pi * 2) elif latnum == len(latitudes) + 1: j = latnum - 2 lonnummax = self.dimensions['lon_number_by_lat'][j] i = ((lonrs - numpy.pi) % (numpy.pi * 2)) * \ lonnummax / (numpy.pi * 2) else: j = latnum lonnummax = self.dimensions['lon_number_by_lat'][latnum] i = lonrs * lonnummax / (numpy.pi * 2) if num == 1: lons = (numpy.rint(i).astype('int') % lonnummax, j) elif num == 2: lons = [(numpy.floor(i).astype('int') % lonnummax, j), ((numpy.floor(i).astype('int') + 1) % lonnummax, j)] elif num == 4: lons = [((numpy.floor(i).astype('int') - 1) % lonnummax, j), (numpy.floor(i).astype('int') % lonnummax, j), ((numpy.floor(i).astype('int') + 1) % lonnummax, j), ((numpy.floor(i).astype('int') + 2) % lonnummax, j)] return lons if self._getoffset(position) != (0., 0.): raise NotImplementedError("horizontal staggered grids for reduced\ gauss grid are not implemented.") if isinstance(lon, list) or isinstance(lon, tuple): lon = numpy.array(lon) if isinstance(lat, list) or isinstance(lat, tuple): lat = numpy.array(lat) latitudes = [self.grid['latitudes'][n].get('radians') for n in range(0, self.dimensions['lat_number'])] if interpolation == 'nearest': def nearest(lon, lat, lonrs, latrs): """ Internal method used to find the nearest point. lon/lat are the true coordinate, lonrs/latrs are rotated and streched coordinates. """ distance = None nearpoint = None for latnum in nearest_lats(latrs, 3): point = nearest_lons(lonrs, latnum, 1) dist = self.distance((lon, lat), self.ij2ll(*point)) if distance == None or dist < distance: nearpoint = point distance = dist return nearpoint (lonrs, latrs) = rotateStretch(lon, lat) if isinstance(lat, numpy.ndarray): j = numpy.zeros(len(lat), dtype=numpy.int) i = numpy.zeros(len(lat), dtype=numpy.int) for k in range(0, len(lat)): (i[k], j[k]) = nearest(lon[k], lat[k], lonrs[k], latrs[k]) else: (i, j) = nearest(lon, lat, lonrs, latrs) points = [(i, j)] elif interpolation in ['linear', 'cubic']: def nearests(lonrs, latrs, n): """ Internal methods used to find the n*n points surrunding the point lonrs/latrs, lonrs/latrs are rotated and stretched coordinates. """ points = [] for j in nearest_lats(latrs, n): points.extend(nearest_lons(lonrs, j, n)) return points (lonrs, latrs) = rotateStretch(lon, lat) n = dict(linear=2, cubic=4)[interpolation] if isinstance(lat, numpy.ndarray): i = numpy.zeros((len(lat), n * n), dtype=numpy.int) j = numpy.zeros((len(lat), n * n), dtype=numpy.int) for k in range(len(lat)): points = nearests(lonrs[k], latrs[k], n) for p in range(n * n): i[k, p], j[k, p] = points[p] points = [(i[:, p], j[:, p]) for p in range(n * n)] else: points = nearests(lonrs, latrs, n) else: raise epygramError("'interpolation' can only be 'nearest', 'linear'\ or 'cubic'.") return points[0] if interpolation == 'nearest' else points