#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright (c) Météo France (2014-)
# This software is governed by the CeCILL-C license under French law.
# http://www.cecill.info
"""
Contains the classes for 3D geometries of fields.
"""
from __future__ import print_function, absolute_import, unicode_literals, division
import numpy
import copy
import sys
import footprints
from epygram import epygramError, config
from epygram.util import Angle, write_formatted, Comparator
from .VGeometry import VGeometry
from .AbstractGeometry import RectangularGridGeometry
epylog = footprints.loggers.getLogger(__name__)
[docs]class AcademicGeometry(RectangularGridGeometry):
"""Handles the geometry for an academic 3-Dimensions Field."""
def __init__(self, name, grid, dimensions, vcoordinate, projection,
position_on_horizontal_grid='__unknown__', geoid=None):
"""
:param name: Name of geometrical type of representation of points on the Globe.
Name must be 'academic'
:param grid: Handles description of the horizontal grid.
:param dimensions: Handles grid dimensions.
:param vcoordinate: Handles vertical geometry parameters.
:param position_on_horizontal_grid: Position of points w/r to the horizontal.
among: ['upper-right', 'upper-left',
'lower-left', 'lower-right',
'center-left', 'center-right',
'lower-center', 'upper-center',
'center', '__unknown__']
:param geoid: To specify geoid shape (of no meaning in this geometry).
"""
self.add_attr_inlist('name', ['academic'])
self.add_attr_dict('projection')
self.name = name
self.projection = projection
super(AcademicGeometry, self).__init__(grid, dimensions, vcoordinate,
position_on_horizontal_grid, geoid)
if self.grid['input_position'] != (0, 0):
raise NotImplementedError("For now, only input_position = (0, 0) is allowed for academic geometries.")
self._center_lon = (self.dimensions['X'] - 1) / 2.
self._center_lat = (self.dimensions['Y'] - 1) / 2.
[docs] def tolerant_equal(self, other, tolerance=config.epsilon):
if self.__class__ == other.__class__:
# create copies of inner objects to filter some ghost attributes
almost_self = {k:copy.deepcopy(self.__dict__[k])
for k in self.__dict__.keys()
if k not in self._ghost_attributes}
almost_other = {k: copy.deepcopy(other.__dict__[k])
for k in other.__dict__.keys()
if k not in other._ghost_attributes}
# (the same grid could be defined from different inputs)
for almost in (almost_self, almost_other):
for k in ('input_lon', 'input_lat', 'input_position'):
almost['_attributes']['grid'].pop(k)
almost_self['_attributes']['grid']['center'] = self.getcenter()
almost_other['_attributes']['grid']['center'] = other.getcenter()
return Comparator.are_equal(almost_self, almost_other, tolerance)
else:
return False
def _rotate_axis(self, x, y, direction):
"""
Internal method used to rotate x/y coordinates to handle rotated geometries.
*direction*, if evals to 'xy2ll', direction used when converting (x, y) to (lat, lon).
if evals to 'll2xy', direction used when converting (lat, lon) to (x, y).
"""
if self.projection['rotation'].get('degrees') == 0:
return (x, y)
else:
beta = self.projection['rotation'].get('radians')
if direction == 'xy2ll':
return (numpy.array(x) * numpy.cos(-beta) + numpy.array(y) * numpy.sin(-beta),
- numpy.array(x) * numpy.sin(-beta) + numpy.array(y) * numpy.cos(-beta))
elif direction == 'll2xy':
return (numpy.array(x) * numpy.cos(beta) + numpy.array(y) * numpy.sin(beta),
- numpy.array(x) * numpy.sin(beta) + numpy.array(y) * numpy.cos(beta))
else:
raise epygramError('Wrong direction of rotation.')
def _consistency_check(self):
"""Check that the geometry is consistent."""
grid_keys = ['LAMzone', 'X_resolution', 'Y_resolution',
'input_lat', 'input_lon', 'input_position']
if set(self.grid.keys()) != set(grid_keys) and \
set(self.grid.keys()) != set(grid_keys + ['longitude', 'latitude']):
raise epygramError("grid attribute must consist in keys: " +
str(grid_keys) + " or " +
str(grid_keys + ['longitude', 'latitude']))
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))
projection_keys = ['reference_dX', 'reference_dY', 'rotation']
if set(self.projection.keys()) != set(projection_keys):
raise epygramError("projection attribute must consist in keys: " +
str(projection_keys))
def _getoffset(self, position=None):
"""
Returns the offset to use for this position.
Replaces the method defined in RectangularGridGeometry to deal with
1D or 2D simulations.
"""
offset = super(AcademicGeometry, self)._getoffset(position)
if self.dimensions['X'] == 1:
offset = (0, offset[1])
if self.dimensions['Y'] == 1:
offset = (offset[0], 0)
return offset
[docs] def getcenter(self):
"""
Returns the coordinate of the grid center as a tuple
(center_lon, center_lat).
"""
return (self._center_lon, self._center_lat)
[docs] def ij2xy(self, i, j, position=None):
"""
Return the (x, y) coordinates of point *(i,j)*, in the projection.
:param i: X index of point in the 2D matrix of gridpoints
:param j: Y index of point in the 2D matrix of gridpoints
:param position: position to return with respect to the model cell.
Defaults to self.position_on_horizontal_grid.
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'])
[docs] def xy2ij(self, x, y, position=None):
"""
Return the (i, j) indexes of point *(x, y)*,
in the 2D matrix of gridpoints.
:param x: X coordinate of point in the academic projection
:param y: Y coordinate of point in the academic projection
:param position: position to return with respect to the model cell.
Defaults to self.position_on_horizontal_grid.
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)
[docs] def ij2ll(self, i, j, position=None):
"""
Return the (lon, lat) coordinates of point *(i,j)*, in degrees.
:param i: X index of point in the 2D matrix of gridpoints
:param j: Y index of point in the 2D matrix of gridpoints
:param position: lat lon position to return with respect to the model cell.
Defaults to self.position_on_horizontal_grid.
This routine has a special meaning for this geometry. It returns the (i, j)
position in the original grid (especially after a section extraction) with
an offset of one (for old tools compatibility).
"""
if isinstance(i, list) or isinstance(i, tuple):
i = numpy.array(i)
if isinstance(j, list) or isinstance(j, tuple):
j = numpy.array(j)
return self.xy2ll(*self.ij2xy(i, j, position))
[docs] def ll2ij(self, lon, lat, position=None):
"""
Return the (i, j) indexes of point *(lon, lat)* in degrees,
in the 2D matrix of gridpoints.
:param lon: longitude of point in degrees
:param lat: latitude of point in degrees
:param position: lat lon position to return with respect to the model cell.
Defaults to self.position_on_horizontal_grid.
Caution: the returned (*i,j*) are float.
This routine has a special meaning for this geometry. It returns the (lon, lat)
position in the original grid (especially after a section extraction) with
an offset of one (for old tools compatibility).
"""
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.xy2ij(*self.ll2xy(lon, lat), position=position)
[docs] def ll2xy(self, lon, lat):
"""
Return the (x, y) coordinates of point *(lon, lat)* in degrees.
:param lon: longitude of point in degrees
:param lat: latitude of point in degrees
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)
# TOBECHECKED: position of self.ij2xy(0, 0) inside the formula
xy = ((lon - self.grid['input_lon']) * self.projection['reference_dX'] + self.ij2xy(0, 0, 'center')[0],
(lat - self.grid['input_lat']) * self.projection['reference_dY'] + self.ij2xy(0, 0, 'center')[1])
return self._rotate_axis(*xy, direction='ll2xy')
[docs] def xy2ll(self, x, y):
"""
Return the (lon, lat) coordinates of point *(x, y)* in the 2D matrix
of gridpoints*(i,j)*, in degrees.
:param x: X coordinate of point in the academic projection
:param y: Y coordinate of point in the academic projection
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)
(a, b) = self._rotate_axis(x, y, direction='xy2ll')
# TOBECHECKED: position of self.ij2xy(0, 0) inside the formula
return ((a - self.ij2xy(0, 0, 'center')[0]) / self.projection['reference_dX'] + self.grid['input_lon'],
(b - self.ij2xy(0, 0, 'center')[1]) / self.projection['reference_dY'] + self.grid['input_lat'])
[docs] def distance(self, end1, end2):
"""
Computes the distance between two points along a straight line in the
geometry.
:param end1: must be a tuple (lon, lat) in degrees.
:param end2: must be a tuple (lon, lat) in degrees.
"""
return numpy.sqrt(((end1[0] - end2[0]) * self.grid['X_resolution']) ** 2 +
((end1[1] - end2[1]) * self.grid['Y_resolution']) ** 2)
[docs] def linspace(self, end1, end2, num):
"""
Returns evenly spaced points over the specified interval.
Points are lined up in the geometry.
:param end1: must be a tuple (lon, lat) in degrees.
:param end2: must be a tuple (lon, lat) in degrees.
:param num: number of points, including point1 and point2.
"""
if num < 2:
raise epygramError("'num' must be at least 2.")
return list(zip(numpy.linspace(end1[0], end2[0], num=num),
numpy.linspace(end1[1], end2[1], num=num)))
[docs] def resolution_ll(self, *_, **__):
"""Returns the minimum of X and Y resolution."""
return min(self.grid['X_resolution'], self.grid['Y_resolution'])
[docs] def resolution_ij(self, *_, **__):
"""Returns the minimum of X and Y resolution."""
return min(self.grid['X_resolution'], self.grid['Y_resolution'])
[docs] def plane_azimuth(self, end1, end2):
"""
Initial bearing from *end1* to *end2* points in plane local referential
geometry.
:param end1: must be a tuple (lon, lat) in degrees.
:param end2: must be a tuple (lon, lat) in degrees.
"""
(x1, y1) = self.ll2xy(*end1)
(x2, y2) = self.ll2xy(*end2)
return (numpy.degrees(numpy.arctan2(x2 - x1, y2 - y1)) + 180.) % 360. - 180.
[docs] def azimuth(self, end1, end2):
"""Same as plane_azimuth in this geometry."""
return self.plane_azimuth(end1, end2)
def _what_position(self, out=sys.stdout):
"""
Writes in file a summary of the position of the field.
:param out: the output open file-like object
"""
write_formatted(out, "Kind of Geometry", 'Academic')
if 'latitude' in self.grid:
write_formatted(out, "Latitude", self.grid['latitude'])
if 'longitude' in self.grid:
write_formatted(out, "Longitude", self.grid['longitude'])
def _what_projection(self, out=sys.stdout, **_):
"""
Writes in file a summary of the projection of the field's grid.
:param out: the output open file-like object
"""
projection = self.projection
write_formatted(out, "X resolution of the reference grid",
projection['reference_dX'])
write_formatted(out, "Y resolution of the reference grid",
projection['reference_dY'])
[docs] def make_section_geometry(self, end1, end2,
points_number=None,
resolution=None,
position=None):
"""
Returns a academic Geometry.
:param end1: must be a tuple (lon, lat) in degrees.
:param end2: must be a tuple (lon, lat) in degrees.
:param points_number: defines the total number of horizontal points of the
section (including ends). If None, defaults to a number computed from
the *ends* and the *resolution*.
:param resolution: defines the horizontal resolution to be given to the
field. If None, defaults to the horizontal resolution of the field.
:param position: defines the position of data in the grid (defaults to 'center')
"""
if resolution is not None and points_number is not None:
raise epygramError("only one of resolution and " +
" points_number can be given.")
(x1, y1) = self.ll2xy(*end1)
(x2, y2) = self.ll2xy(*end2)
distance = numpy.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2)
if resolution is None and points_number is None:
resolution = 0.5 * (self.resolution_ll(*end1) +
self.resolution_ll(*end2))
points_number = int(numpy.rint(distance / resolution)) + 1
if resolution > distance:
raise epygramError("'ends' are too near: pure" +
" interpolation between two gridpoints.")
if points_number is not None:
if points_number < 2:
raise epygramError("'points_number' must be at least 2.")
resolution = distance / (points_number - 1)
else:
points_number = int(numpy.rint(distance / resolution)) + 1
vcoordinate = VGeometry(typeoffirstfixedsurface=255,
levels=[])
grid = {'LAMzone':'CIE',
'X_resolution':resolution,
'Y_resolution':resolution,
'input_lat':end1[1],
'input_lon':end1[0],
'input_position':(0, 0)}
if 'longitude' in self.grid:
grid['longitude'] = self.grid['longitude']
if 'latitude' in self.grid:
grid['latitude'] = self.grid['latitude']
dimensions = {'X':points_number, 'Y':1,
'X_Iwidth':0, 'Y_Iwidth':0,
'X_Czone':0, 'Y_Czone':0,
'X_CIzone':points_number, 'Y_CIzone':1,
'X_CIoffset':0, 'Y_CIoffset':0}
rotation = numpy.arctan2(y2 - y1, x2 - x1)
projection = {'rotation':Angle(rotation, 'radians'),
'reference_dX':self.projection['reference_dX'],
'reference_dY':self.projection['reference_dY']}
kwargs_geom = dict(name=self.name,
grid=grid,
projection=projection,
dimensions=dimensions,
position_on_horizontal_grid='center' if position is None else position,
vcoordinate=vcoordinate
)
if self.geoid:
kwargs_geom['geoid'] = self.geoid
return AcademicGeometry(**kwargs_geom)