#!/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 sys
import footprints
from bronx.syntax.arrays import stretch_array
from epygram.util import (Angle,
positive_longitudes, longitudes_between_minus180_180,
write_formatted,
as_numpy_array, moveaxis)
from .AbstractGeometry import RectangularGridGeometry
epylog = footprints.loggers.getLogger(__name__)
[docs]class UnstructuredGeometry(RectangularGridGeometry):
"""Handles the geometry for an unstructured 3-Dimensions Field."""
def __init__(self, name, grid, dimensions, vcoordinate,
position_on_horizontal_grid='center', geoid=None):
"""
:param name: Name of geometrical type of representation of points on the Globe.
Name must be 'unstructured'
: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.
"""
self.add_attr_inlist('name', ['unstructured',
'DDH:point', 'DDH:ij_point', 'DDH:quadrilateral',
'DDH:rectangle', 'DDH:globe', 'DDH:zonal_bands'])
self.name = name
super(UnstructuredGeometry, self).__init__(grid, dimensions, vcoordinate,
position_on_horizontal_grid, geoid)
def _consistency_check(self):
"""Check that the geometry is consistent."""
grid_keys = ['latitudes', 'longitudes']
if set(self.grid.keys()) != set(grid_keys) and \
set(self.grid.keys()) != set(['DDH_domain']):
raise epygramError("grid attribute must consist in keys: " +
str(grid_keys) + " or " +
str(['DDH_domain']))
[docs] def get_lonlat_grid(self,
subzone=None,
position=None,
d4=False,
nb_validities=0,
force_longitudes=None):
"""
Returns a tuple of two tables containing one the longitude of each
point, the other the latitude, with 2D shape.
:param 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.
:param position: position of lonlat grid with respect to the model cell.
Defaults to self.position_on_horizontal_grid.
:param d4: - if True, returned values are shaped in a 4 dimensions array
- if False, shape of returned values is determined with respect to geometry.
d4=True requires nb_validities > 0
:param nb_validities: number of validities represented in data values
:param force_longitudes: if 'positive', the longitudes will be forced positive
if ']-180,180]', the longitudes will be in the ]-180, 180] interval
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
"""
if self._getoffset(position) != (0., 0.):
raise epygramError('We can only retrieve latitude and longitude of mass point on an unstructured grid')
lons = self.grid['longitudes']
if not isinstance(lons, numpy.ndarray):
lons = numpy.array(lons)
lats = self.grid['latitudes']
if not isinstance(self.grid['latitudes'], numpy.ndarray):
lats = numpy.array(lats)
if len(lons.shape) == 1:
lons = self.reshape_data(lons)
lats = self.reshape_data(lats)
if subzone and self.grid.get('LAMzone', None):
lons = self.extract_subzone(lons, subzone)
lats = self.extract_subzone(lats, subzone)
if d4:
lons, lats = self._reshape_lonlat_4d(lons, lats, nb_validities)
elif not d4 and nb_validities != 0:
raise ValueError("*nb_validities* must be 0 when d4==False")
else:
lons = lons.squeeze()
lats = lats.squeeze()
if force_longitudes == 'positive':
lons = positive_longitudes(lons)
elif force_longitudes == ']-180,180]':
lons = longitudes_between_minus180_180(lons)
return (lons, lats)
[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.
"""
(lons, lats) = self.get_lonlat_grid(position=position)
if lons.ndim == 1:
assert numpy.all(j == 0), "j must be 0 when geometry does not contain j dimension"
return (lons[i], lats[i])
else:
return (lons[j, i], lats[j, i])
[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.
"""
lon, lat = as_numpy_array(lon).flatten(), as_numpy_array(lat).flatten()
i = numpy.zeros(len(lon))
j = numpy.zeros(len(lon))
(lons, lats) = self.get_lonlat_grid(position=position)
for n in range(len(lon)):
where = (numpy.abs(lons - lon[n]) + numpy.abs(lats - lat[n]) == 0).nonzero()
if len(where[0]) == 0:
raise epygramError("No point found with these coordinates.")
elif len(where[0]) > 1:
raise epygramError("Several points have the same coordinates.")
if self.datashape['j']:
i[n], j[n] = where[::-1]
else:
i[n] = where[0]
return (i.squeeze(), j.squeeze())
[docs] def nearest_points(self, lon, lat, request,
position=None,
external_distance=None,
squeeze=True):
"""
Returns the (i, j) position of the points needed to perform an
interpolation. This is a list of (lon, lat) tuples.
:param lon: longitude of point in degrees.
:param lat: latitude of point in degrees.
:param request: criteria for selecting the points, among:
* {'n':'1'} - the nearest point
:param position: position in the model cell of the lat lon position.
Defaults to self.position_on_horizontal_grid.
:param external_distance: can be a dict containing the target point value
and an external field on the same grid as self, to which the distance
is computed within the 4 horizontally nearest points; e.g.
{'target_value':4810, 'external_field':a_3DField_with_same_geometry}.
If so, the nearest point is selected with
distance = |target_value - external_field.data|
:param squeeze: True to suppress useless dimensions
:rtype: general output form is [list, list, ..., list]
with as many list items as the length of lon/lat.
Each list item is of the form [tuple, tuple, ..., tuple]
with as many tuples as the request implies. A tuple
represents one of the nearest points associated with one
value taken from lon/lat. Each tuple as the form
(i, j).
Dimensions with a length of one are removed except if
squeeze is False. If squeeze is True and if request
implies only one nearest point, the list item of the general
output form is replaced by the tuple item; if length of
lon/lat is one, the output is directly the list item of
the general output form. Hence, if length of lon/lat is
one and the request implies only one point, the output is
a tuple.
In case of a simple square request, output is actually
an array. Otherwise, the output is as described (it cannot
be an array because the number of nearest points can vary
with the entry point).
"""
if request != {'n':'1'}:
raise NotImplementedError("*request* != '{'n':'1'}' for UnstructuredGeometry.")
if external_distance is not None:
raise NotImplementedError("*external_distance* is not None for UnstructuredGeometry.")
lon, lat = as_numpy_array(lon).flatten(), as_numpy_array(lat).flatten()
(lons, lats) = self.get_lonlat_grid(position=position)
result = []
for (one_lon, one_lat) in zip(lon, lat):
dist = (one_lon - lons) ** 2 + (one_lat - lats) ** 2
i = dist.argmin()
if lons.ndim == 2:
result.append(numpy.unravel_index(i, lons.shape)[::-1])
else:
result.append((i, 0))
result = numpy.array(result)
if squeeze:
result = result.squeeze()
return result
[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.
:param lon: longitude of point in degrees.
:param lat: latitude of point in degrees.
"""
near_point = moveaxis(self.nearest_points(lon, lat, request={'n':'1'}), 0, -1)
return self.resolution_ij(*near_point)
[docs] def resolution_ij(self, i, j, position=None):
"""
Returns the distance to the nearest point of (i,j) point.
:param i: X index of point in the 2D matrix of gridpoints
:param j: Y index of point in the 2D matrix of gridpoints
"""
(lons, lats) = self.get_lonlat_grid(position=position)
i, j = as_numpy_array(i), as_numpy_array(j)
assert len(i) == len(j), "Both coordinates must have the same length"
result = numpy.zeros(len(i))
for k in range(len(i)):
dist = self.distance((lons[j[k], i[k]], lats[j[k], i[k]]),
(stretch_array(lons), stretch_array(lats)))
result[k] = dist[dist != 0].min()
return result.squeeze()
[docs] def getcenter(self, position=None):
"""
Returns the coordinate of the grid center as a tuple of Angles
(center_lon, center_lat).
Caution: this is computed as the raw average of all grid points.
A barycentric computation would be more adequate.
"""
(lons, lats) = self.get_lonlat_grid(position=position)
return (Angle(lons.mean(), 'degrees'),
Angle(lats.mean(), 'degrees'))
def _what_grid(self, out=sys.stdout):
"""
Writes in file a summary of the grid of the field.
:param out: the output open file-like object
"""
(lons, lats) = self.get_lonlat_grid()
corners = self.gimme_corners_ll()
write_formatted(out, "Kind of Geometry", 'Unstructured')
write_formatted(out, "Max Longitude in deg", lons.max())
write_formatted(out, "Min Longitude in deg", lons.min())
write_formatted(out, "Max Latitude in deg", lats.max())
write_formatted(out, "Min Latitude in deg", lats.min())
write_formatted(out, "Low-Left corner Longitude in deg",
corners['ll'][0])
write_formatted(out, "Low-Left corner Latitude in deg",
corners['ll'][1])
write_formatted(out, "Low-Right corner Longitude in deg",
corners['lr'][0])
write_formatted(out, "Low-Right corner Latitude in deg",
corners['lr'][1])
write_formatted(out, "Upper-Left corner Longitude in deg",
corners['ul'][0])
write_formatted(out, "Upper-Left corner Latitude in deg",
corners['ul'][1])
write_formatted(out, "Upper-Right corner Longitude in deg",
corners['ur'][0])
write_formatted(out, "Upper-Right corner Latitude in deg",
corners['ur'][1])