Source code for epygram.geometries.ProjectedGeometry

#!/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 math
import copy
import sys

import footprints
from footprints import proxy as fpx

from epygram import epygramError, config
from epygram.config import rounding_decimal as _rd
from epygram.util import (degrees_nearest_mod, Angle,
                          write_formatted,
                          as_numpy_array,
                          is_scalar, Comparator)

from .VGeometry import VGeometry
from .AbstractGeometry import RectangularGridGeometry

epylog = footprints.loggers.getLogger(__name__)



[docs]class ProjectedGeometry(RectangularGridGeometry): """ Handles the geometry for a Projected 3-Dimensions Field. """ _ghost_attributes = RectangularGridGeometry._ghost_attributes + ['_proj'] 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 among ['lambert', 'mercator', 'polar_stereographic', 'space_view'] :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', ['lambert', 'mercator', 'polar_stereographic', 'space_view']) self.add_attr_dict('projection') self.name = name self.projection = projection super(ProjectedGeometry, self).__init__(grid, dimensions, vcoordinate, position_on_horizontal_grid, geoid) import pyproj def compute_center_proj(p, center): if center == self.grid['input_position']: self._center_lon = self.grid['input_lon'] self._center_lat = self.grid['input_lat'] else: # x1, y1: coordinates in non rotated proj of input point x1, y1 = p(float(self.grid['input_lon'].get('degrees')), float(self.grid['input_lat'].get('degrees'))) # offset between center and input points is known in rotated proj # dx, dy is the offset in non rotated proj (dx, dy) = self._rotate_axis( round((center[0] - self.grid['input_position'][0]) * self.grid['X_resolution'], _rd), round((center[1] - self.grid['input_position'][1]) * self.grid['Y_resolution'], _rd), 'xy2ll') # xc, yc: coordinates of center point in non rotated proj xc = x1 + dx yc = y1 + dy center_lon, center_lat = p(xc, yc, inverse=True) self._center_lon = Angle(center_lon, 'degrees') self._center_lat = Angle(center_lat, 'degrees') projdict = {'lambert':'lcc', 'mercator':'merc', 'polar_stereographic':'stere', 'space_view':'geos'} proj = projdict[self.name] # build proj if self.grid['LAMzone'] is not None: nx = self.dimensions['X_CIzone'] ny = self.dimensions['Y_CIzone'] io = self.dimensions.get('X_CIoffset', 0) jo = self.dimensions.get('Y_CIoffset', 0) centerPoint = (io + (float(nx) - 1) / 2., jo + (float(ny) - 1) / 2.) # Coordinates of center point else: nx = self.dimensions['X'] ny = self.dimensions['Y'] centerPoint = ((float(nx) - 1) / 2., (float(ny) - 1) / 2.) # Coordinates of center point # 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 = pyproj.Proj(proj=proj, lon_0=self.projection['reference_lon'].get('degrees'), lat_1=lat_1, lat_2=lat_2, **self.geoid) compute_center_proj(p, centerPoint) x0, y0 = p(self._center_lon.get('degrees'), self._center_lat.get('degrees')) self._proj = pyproj.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 = pyproj.Proj(proj=proj, lon_0=self.projection['reference_lon'].get('degrees'), lat_ts=lat_ts, **self.geoid) compute_center_proj(p, centerPoint) x0, y0 = p(self._center_lon.get('degrees'), self._center_lat.get('degrees')) self._proj = pyproj.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 = pyproj.Proj(proj=proj, lon_0=self.projection['reference_lon'].get('degrees'), lat_0=lat_0, lat_ts=lat_ts, **self.geoid) compute_center_proj(p, centerPoint) x0, y0 = p(self._center_lon.get('degrees'), self._center_lat.get('degrees')) self._proj = pyproj.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) elif self.name == 'space_view': latSat = self.projection['satellite_lat'].get('degrees') lonSat = self.projection['satellite_lon'].get('degrees') height = self.projection['satellite_height'] # Height above ellipsoid if latSat != 0: raise epygramError("Only space views with satellite_lat=0 are allowed") p = pyproj.Proj(proj=proj, h=height, lon_0=lonSat, **self.geoid) compute_center_proj(p, centerPoint) x0, y0 = p(self._center_lon.get('degrees'), self._center_lat.get('degrees')) self._proj = pyproj.Proj(proj=proj, h=height, lon_0=lonSat, x_0=-x0, y_0=-y0, **self.geoid) else: raise NotImplementedError("projection: " + self.name) @property 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)
[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['grid'].pop(k) almost_self['grid']['center'] = self.getcenter() almost_other['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. :param 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.""" if self.name == 'lambert': sets_of_keys = (['reference_lon', 'reference_lat', 'rotation'], ['reference_lon', 'secant_lat1', 'secant_lat2', 'rotation']) elif self.name in ('polar_stereographic', 'mercator'): sets_of_keys = (['reference_lon', 'reference_lat', 'rotation'], ['reference_lon', 'reference_lat', 'secant_lat', 'rotation']) elif self.name == 'space_view': sets_of_keys = ['satellite_lat', 'satellite_lon', 'satellite_height', 'rotation', 'reference_lon'] sets_of_keys = (sets_of_keys, sets_of_keys) else: raise NotImplementedError("projection: " + self.name) 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])) for a in ['reference_lon', 'reference_lat', 'rotation', 'secant_lat1', 'secant_lat2', 'secant_lat', 'satellite_lat', 'satellite_lon']: assert isinstance(self.projection.get(a, Angle(0., 'degrees')), Angle) grid_keys = ['LAMzone', 'X_resolution', 'Y_resolution', 'input_lon', 'input_lat', 'input_position'] if set(self.grid.keys()) != set(grid_keys): raise epygramError("grid attribute must consist in keys: " + str(grid_keys)) assert isinstance(self.grid['input_lon'], Angle) assert isinstance(self.grid['input_lat'], Angle) 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)) # if self.projection['rotation'].get('degrees') != 0.0: # epylog.warning('*rotation* != 0. may not have been thoroughly tested...') # TOBECHECKED: here and there, ...
[docs] def select_subzone(self, subzone): """ If a LAMzone defines the geometry, select only the *subzone* from it and return a new geometry object. :param subzone: among ('C', 'CI'). """ assert subzone in ('C', 'CI'), \ 'unknown subzone : ' + subzone if self.grid.get('LAMzone') is not None: newgeom = self.deepcopy() if subzone == 'CI': io = self.dimensions.get('X_CIoffset', 0) jo = self.dimensions.get('Y_CIoffset', 0) centerPoint = (io + (float(self.dimensions['X_CIzone']) - 1) / 2., jo + (float(self.dimensions['Y_CIzone']) - 1) / 2.) # Coordinates of center point newgeom.grid['LAMzone'] = subzone newgeom.dimensions = {'X':self.dimensions['X_CIzone'], 'Y':self.dimensions['Y_CIzone'], 'X_CIzone':self.dimensions['X_CIzone'], 'Y_CIzone':self.dimensions['Y_CIzone'], 'X_Iwidth':self.dimensions['X_Iwidth'], 'Y_Iwidth':self.dimensions['Y_Iwidth'], 'X_Czone':self.dimensions['X_Czone'], 'Y_Czone':self.dimensions['Y_Czone']} elif subzone == 'C': centerPoint = ((float(self.dimensions['X_Czone']) - 1) / 2., (float(self.dimensions['Y_Czone']) - 1) / 2.) # Coordinates of center point newgeom.grid['LAMzone'] = None newgeom.dimensions = {'X':self.dimensions['X_Czone'], 'Y':self.dimensions['Y_Czone']} newgeom.grid['input_lon'] = self._center_lon newgeom.grid['input_lat'] = self._center_lat newgeom.grid['input_position'] = centerPoint else: newgeom = self return newgeom
[docs] def getcenter(self): """ Returns the coordinate of the grid center as a tuple of Angles (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) Xresolution = self.grid['X_resolution'] Yresolution = self.grid['Y_resolution'] if self.grid.get('LAMzone', None) is 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.get('X_CIoffset', 0) + float(Xpoints - 1) / 2. j0 = self.dimensions.get('Y_CIoffset', 0) + 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): """ Return the (i, j) indexes of point *(x, y)*, in the 2D matrix of gridpoints. :param x: X coordinate of point in the projection :param y: Y coordinate of point in the 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) Xresolution = self.grid['X_resolution'] Yresolution = self.grid['Y_resolution'] if self.grid.get('LAMzone', None) is 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.get('X_CIoffset', 0) + float(Xpoints - 1) / 2. j0 = self.dimensions.get('Y_CIoffset', 0) + 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): """ 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) lon = degrees_nearest_mod(lon, self.projection['reference_lon'].get('degrees')) xy = self._proj(lon, lat) 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 projection :param y: Y coordinate of point in the 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') ll = self._proj(a, b, inverse=True) # mask invalid values lons, lats = ll lons_is_scalar, lats_is_scalar = is_scalar(lons), is_scalar(lats) lons, lats = as_numpy_array(lons), as_numpy_array(lats) mask = numpy.logical_or(lons == config.pyproj_invalid_values, lats == config.pyproj_invalid_values) lons = numpy.ma.masked_where(mask, lons) lats = numpy.ma.masked_where(mask, lats) if lons_is_scalar: lons = lons[0] if lats_is_scalar: lats = lats[0] ll = (lons, lats) return ll
[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. """ 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: (*i,j*) are float. """ return self.xy2ij(*self.ll2xy(lon, lat), position=position)
[docs] def map_factor(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. :param position: grid position with respect to the model cell. Defaults to self.position_on_horizontal_grid. """ kwargs_vcoord = {'typeoffirstfixedsurface': 255, 'position_on_grid': '__unknown__', 'levels': [0]} vcoordinate = VGeometry(**kwargs_vcoord) geometry = self.deepcopy() geometry.vcoordinate=vcoordinate f = fpx.field(structure='H2D', geometry=geometry, fid={'geometry':'Map Factor'}, units='-') lats = self.get_lonlat_grid(position=position)[1] data = self.map_factor(lats) f.setdata(data) return f
[docs] def mesh_area(self, lat): """ Compute the area of a mesh/gridpoint, given its latitude. """ return self.grid['X_resolution'] * self.grid['Y_resolution'] / (self.map_factor(lat)**2)
[docs] def mesh_area_field(self, position=None): """ Returns a new field whose data is the mesh area of gridpoints, i.e. X_resolution x Y_resolution / m^2, where m is the local map factor. :param position: grid position with respect to the model cell. Defaults to self.position_on_horizontal_grid. """ m = self.map_factor_field(position=position) dxdy = self.grid['X_resolution'] * self.grid['Y_resolution'] m.setdata(dxdy / m.data**2) m.fid['geometry'] = 'Mesh Area' m.units = 'm^2' return m
[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.") (x1, y1) = self.ll2xy(*end1) (x2, y2) = self.ll2xy(*end2) xy_linspace = list(zip(numpy.linspace(x1, x2, num=num), numpy.linspace(y1, y2, num=num))) return [tuple(numpy.around(self.xy2ll(*xy), 8)) for xy in xy_linspace]
[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. """ (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 100 points along the line # between point1 and point2 mean_map_factor = numpy.array([self.map_factor(lat) for (_, lat) in self.linspace(end1, end2, 10)]).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. :param lon: longitude of the point in degrees :param lat: latitude of the point in degrees """ 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. :param i: X index of point in the 2D matrix of gridpoints :param j: Y index of point 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 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) beta = self.projection['rotation'].get('degrees') return (numpy.degrees(numpy.arctan2(x2 - x1, y2 - y1)) - beta + 180.) % 360. - 180.
[docs] def make_section_geometry(self, end1, end2, points_number=None, resolution=None, position=None): """ Returns a projected 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 rotation = numpy.arctan2(y2 - y1, x2 - x1) + self.projection['rotation'].get('radians') vcoordinate = VGeometry(typeoffirstfixedsurface=255, levels=[]) grid = {'LAMzone':None, 'X_resolution':resolution, 'Y_resolution':resolution, 'input_lon':Angle(end1[0], 'degrees'), 'input_lat':Angle(end1[1], 'degrees'), 'input_position':(0, 0)} projection = dict(self.projection) projection['rotation'] = Angle(rotation, 'radians') dimensions = {'X':points_number, 'Y':1} kwargs_geom = dict(name=self.name, grid=grid, dimensions=dimensions, projection=projection, geoid=self.geoid, position_on_horizontal_grid='center' if position is None else position, vcoordinate=vcoordinate) if self.geoid: kwargs_geom['geoid'] = self.geoid return ProjectedGeometry(**kwargs_geom)
[docs] def compass_grid(self, subzone=None, position=None): """ Get the compass grid, i.e. the angle between Y-axis and North for each gridpoint. :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. """ (lons, _) = self.get_lonlat_grid(subzone=subzone, position=position) if self.name == 'lambert': if self.secant_projection: k = numpy.copysign(self._K, self.projection['secant_lat1'].get('degrees') + self.projection['secant_lat1'].get('degrees')) else: k = numpy.copysign(self._K, self.projection['reference_lat'].get('degrees')) elif self.name == 'mercator': k = 0. elif self.name == 'polar_stereographic': k = self.projection['reference_lat'].get('cos_sin')[1] theta = k * (lons - self.projection['reference_lon'].get('degrees')) - \ self.projection.get('rotation', 0.).get('degrees') return theta
[docs] def reproject_wind_on_lonlat(self, u, v, lon=None, lat=None, map_factor_correction=True, reverse=False): """ Reprojects a wind vector (u, v) on the grid onto real axes, i.e. with components on true zonal/meridian axes. :param u: the u == zonal-on-the-grid component of wind :param v: the v == meridian-on-the-grid component of wind :param lon: longitudes of points in degrees, if u/v are not vectors on the whole grid :param lat: latitudes of points in degrees, if u/v are not vectors on the whole grid :param map_factor_correction:, applies a correction of magnitude due to map factor. :param reverse: if True, apply the reverse reprojection. """ theta = (-1) ** int(reverse) * self.compass_grid() costheta = numpy.cos(-theta * numpy.pi / 180.) sintheta = numpy.sin(-theta * numpy.pi / 180.) if numpy.shape(u) == costheta.shape: # whole grid if map_factor_correction: m = self.map_factor_field().getdata() else: m = numpy.ones(costheta.shape) if reverse: m = 1. / m ru = m * (u * costheta - v * sintheta) rv = m * (u * sintheta + v * costheta) else: # some points ru = numpy.ndarray(numpy.shape(u)) rv = numpy.ndarray(numpy.shape(v)) for k in range(len(u)): if map_factor_correction: m = self.map_factor(lat[k]) else: m = 1. if reverse: m = 1. / m (i, j) = self.ll2ij(lon[k], lat[k]) i = int(i) j = int(j) c = costheta[j, i] s = sintheta[j, i] ru[k] = m * (u[k] * c - v[k] * s) rv[k] = m * (u[k] * s + v[k] * c) return (ru, rv)
def _what_projection(self, out=sys.stdout, arpifs_var_names=False): """ Writes in file a summary of the projection of the field's grid. :param out: the output open file-like object :param arpifs_var_names: if True, prints the equivalent 'arpifs' variable names. """ projection = self.projection grid = self.grid dimensions = self.dimensions varname = '' projmap = {'regular_lonlat':'Regular Lon/Lat', 'lambert':'Lambert (conformal conic)', 'mercator':'Mercator', 'polar_stereographic':'Polar Stereographic', 'space_view':'Space View'} write_formatted(out, "Kind of projection", projmap[self.name]) if self.name == 'space_view': write_formatted(out, "Satellite Longitude in deg", projection['satellite_lon'].get('degrees')) write_formatted(out, "Satellite Latitude in deg", projection['satellite_lat'].get('degrees')) write_formatted(out, "Reference Longitude in deg", projection['reference_lon'].get('degrees')) else: if self.secant_projection: if self.name == 'lambert': write_formatted(out, "Secant Latitude 1 in deg", projection['secant_lat1'].get('degrees')) write_formatted(out, "Secant Latitude 2 in deg", projection['secant_lat2'].get('degrees')) else: write_formatted(out, "Secant Latitude in deg", projection['secant_lat'].get('degrees')) else: write_formatted(out, "Sinus of Reference Latitude", projection['reference_lat'].get('cos_sin')[1]) if arpifs_var_names: varname = ' (ELAT0)' write_formatted(out, "Reference Latitude in deg" + varname, projection['reference_lat'].get('degrees')) if arpifs_var_names: varname = ' (ELON0)' write_formatted(out, "Reference Longitude in deg" + varname, projection['reference_lon'].get('degrees')) write_formatted(out, "Angle of rotation in deg", projection['rotation'].get('degrees')) if self.grid.get('LAMzone', False) in ('C', False): (lons, lats) = self.get_lonlat_grid() corners = self.gimme_corners_ll() else: (lons, lats) = self.get_lonlat_grid(subzone='CI') corners = self.gimme_corners_ll(subzone='CI') if arpifs_var_names: varname = ' (ELONC)' write_formatted(out, "Center Longitude (of C+I) in deg" + varname, self._center_lon.get('degrees')) if arpifs_var_names: varname = ' (ELATC)' write_formatted(out, "Center Latitude (of C+I) in deg" + varname, self._center_lat.get('degrees')) if arpifs_var_names: varname = ' (EDELX)' write_formatted(out, "Resolution in X, in metres" + varname, grid['X_resolution']) if arpifs_var_names: varname = ' (EDELY)' write_formatted(out, "Resolution in Y, in metres" + varname, grid['Y_resolution']) if self.grid['LAMzone'] is None: dimX = dimensions['X'] dimY = dimensions['Y'] else: dimX = dimensions['X_CIzone'] dimY = dimensions['Y_CIzone'] if arpifs_var_names: varname = ' (ELX)' write_formatted(out, "Domain width (of C+I) in X, in metres" + varname, grid['X_resolution'] * (dimX - 1)) if arpifs_var_names: varname = ' (ELY)' write_formatted(out, "Domain width (of C+I) in Y, in metres" + varname, grid['Y_resolution'] * (dimY - 1)) write_formatted(out, "Max Longitude (of C+I) in deg", lons.max()) write_formatted(out, "Min Longitude (of C+I) in deg", lons.min()) write_formatted(out, "Max Latitude (of C+I) in deg", lats.max()) write_formatted(out, "Min Latitude (of C+I) in deg", lats.min()) write_formatted(out, "Low-Left corner (of C+I) Longitude in deg", corners['ll'][0]) write_formatted(out, "Low-Left corner (of C+I) Latitude in deg", corners['ll'][1]) write_formatted(out, "Low-Right corner (of C+I) Longitude in deg", corners['lr'][0]) write_formatted(out, "Low-Right corner (of C+I) Latitude in deg", corners['lr'][1]) write_formatted(out, "Upper-Left corner (of C+I) Longitude in deg", corners['ul'][0]) write_formatted(out, "Upper-Left corner (of C+I) Latitude in deg", corners['ul'][1]) write_formatted(out, "Upper-Right corner (of C+I) Longitude in deg", corners['ur'][0]) write_formatted(out, "Upper-Right corner (of C+I) Latitude in deg", corners['ur'][1]) # FIXME: cleanme def __eq__(self, other): # """Test of equality by recursion on the object's attributes.""" """ if self.__class__ == other.__class__ and \ set(self._attributes.keys()) == set(other._attributes.keys()): for attr in self._attributes.keys(): if attr == 'grid': # the same grid could be defined from different inputs selfgrid = {k:v for k, v in self.grid.items() if k not in ['input_lon', 'input_lat', 'input_position']} othergrid = {k:v for k, v in self.grid.items() if k not in ['input_lon', 'input_lat', 'input_position']} ok = (selfgrid == othergrid and self.getcenter() == other.getcenter()) else: ok = self._attributes[attr] == other._attributes[attr] if not ok: break else: ok = False return ok def __hash__(self): # known issue __eq__/__hash__ must be defined both or none, else inheritance is broken return super(ProjectedGeometry, self).__hash__()"""