#!/usr/bin/env python
# -*- coding: utf-8 -*-
import numpy
import math
import footprints
from footprints import FootprintBase, FPDict, proxy as fpx
from epygram.util import RecursiveObject, degrees_nearest_mod
from epygram import epygramError, config, epylog
[docs]class H2DGeometry(RecursiveObject, FootprintBase):
"""
Handles the geometry for a Horizontal 2-Dimensions Field.
"""
_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__'})),
)
)
@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()
[docs]class H2DRectangularGridGeometry(H2DGeometry):
"""
Handles the geometry for a Horizontal 2-Dimensions Field.
"""
_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):
"""
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.
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)
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 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 == 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 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):
"""
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.
"""
corners = self.gimme_corners_ij(subzone=subzone)
for c in corners.keys():
corners[c] = self.ij2ll(*corners[c])
return corners
[docs] def point_is_inside_domain(self, lon, lat, 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(lon)
except Exception:
N = 1
if N == 1:
inside = Xmin + margin <= self.ll2ij(lon, lat)[0] <= Xmax - margin \
and Ymin + margin <= self.ll2ij(lon, lat)[1] <= Ymax - margin
else:
inside = []
for i in range(N):
inside.append(Xmin + margin <= self.ll2ij(lon[i], lat[i])[0] <= Xmax - margin \
and Ymin + margin <= self.ll2ij(lon[i], lat[i])[1] <= Ymax - margin)
return inside
class H2DAcademicGeometry(H2DRectangularGridGeometry):
"""
Handles the geometry for a 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 ij2ll(self, i, j):
"""
(*i, j*) being the coordinates in the 2D matrix of gridpoints, \n
(*lon, lat*) being the lon/lat coordinates in degrees.
!!! This routine has no sense for this geometry, it is identity.
"""
return (i+1, j+1)
def ll2ij(self, lon, lat):
"""
(*lon, lat*) being the lon/lat coordinates in degrees, \n
(*i, j*) being the coordinates in the 2D matrix of gridpoints.
!!! This routine has no sense for this geometry, it is identity.
"""
return (lon-1, lat-1)
[docs]class H2DRegLLGeometry(H2DRectangularGridGeometry):
"""
Handles the geometry for a 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):
"""
(*i, j*) being the coordinates in the 2D matrix of gridpoints, \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(i, list) or isinstance(i, tuple):
i = numpy.array(i)
if isinstance(j, list) or isinstance(j, tuple):
j = numpy.array(j)
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')
i0 = float(Xpoints-1)/2. # origin of coordinates is the center of domain
j0 = float(Ypoints-1)/2. # origin of coordinates is the center of domain
x = Xorigin + (i-i0) * Xresolution
y = Yorigin + (j-j0) * Yresolution
return (x,y)
[docs] def xy2ij(self, x, y):
"""
(*x, y*) being the coordinates in the projection, \n
(*i, j*) being the coordinates in the 2D matrix of gridpoints.
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)
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')
i0 = float(Xpoints-1)/2. # origin of coordinates is the center of domain
j0 = float(Ypoints-1)/2. # origin of coordinates is the center of domain
i = i0 + (x-Xorigin)/Xresolution
j = j0 + (y-Yorigin)/Yresolution
return (i,j)
[docs] def ij2ll(self, i, j):
"""
(*i, j*) being the coordinates in the 2D matrix of gridpoints, \n
(*lon, lat*) being the lon/lat coordinates in degrees.
"""
return self.ij2xy(i, j)
[docs] def ll2ij(self, lon, lat):
"""
(*lon, lat*) being the lon/lat coordinates in degrees, \n
(*i, j*) being the coordinates in the 2D matrix of gridpoints.
Caution: (*i,j*) are float.
"""
lon = degrees_nearest_mod(lon, self.grid['center_lon'].get('degrees'))
return self.xy2ij(lon, lat)
[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]class H2DProjectedGeometry(H2DRectangularGridGeometry):
"""
Handles the geometry for a 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):
"""
(*i, j*) being the coordinates in the 2D matrix of gridpoints, \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(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.
i0 = float(Xpoints-1)/2. # origin of coordinates is the center of domain
j0 = float(Ypoints-1)/2. # origin of coordinates is the center of domain
x = Xorigin + (i-i0) * Xresolution
y = Yorigin + (j-j0) * Yresolution
return (x,y)
[docs] def xy2ij(self, x, y):
"""
(*x, y*) being the coordinates in the projection, \n
(*i, j*) being the coordinates in the 2D matrix of gridpoints.
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.
i0 = float(Xpoints-1)/2. # origin of coordinates is the center of domain
j0 = float(Ypoints-1)/2. # origin of coordinates is the center of domain
i = i0 + (x-Xorigin)/Xresolution
j = j0 + (y-Yorigin)/Yresolution
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):
"""
(*i, j*) being the coordinates in the 2D matrix of gridpoints, \n
(*lon, lat*) being the lon/lat coordinates in degrees.
"""
return self.xy2ll(*self.ij2xy(i, j))
[docs] def ll2ij(self, lon, lat):
"""
(*lon, lat*) being the lon/lat coordinates in degrees, \n
(*i, j*) being the coordinates in the 2D matrix of gridpoints.
Caution: (*i,j*) are float.
"""
return self.xy2ij(*self.ll2xy(lon, lat))
[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')
K = self._K
m = (numpy.cos(lat_0)/numpy.cos(lat))**(1.-K) \
* ((1.+numpy.copysign(1., lat_0)*numpy.sin(lat_0)) \
/ (1.+numpy.copysign(1., lat_0)*numpy.sin(lat)))**K
else:
raise epygramError("projection " + self.name + " not implemented.")
return m
[docs] def map_factor_field(self):
"""
Returns a new field whose data is the map factor over the field.
"""
f = fpx.field(geometry=self, fid=FPDict({'geometry':'Map Factor'}))
lats = self.get_lonlat_grid()[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
# corners
if self.name in ('lambert', 'mercator'):
(llcrnrlon, llcrnrlat) = self.ij2ll(*self.gimme_corners_ij(subzone)['ll'])
(urcrnrlon, urcrnrlat) = self.ij2ll(*self.gimme_corners_ij(subzone)['ur'])
elif self.name == 'polar_stereographic':
if self.secant_projection:
lat = 'secant_lat'
else:
lat = 'reference_lat'
if self.projection[lat].get('degrees') > 0.: psproj = 'npstere'
elif self.projection[lat].get('degrees') < 0.: psproj = 'spstere'
else: raise NotImplementedError("universal polar stereographic.")
latgrid = self.get_lonlat_grid(subzone=subzone)[1]
boundinglat = numpy.concatenate((latgrid[0,:], latgrid[-1,:],
latgrid[:,0], latgrid[:,-1]))
(llcrnrlon, llcrnrlat) = self.ij2ll(*self.gimme_corners_ij(subzone)['ll'])
(urcrnrlon, urcrnrlat) = self.ij2ll(*self.gimme_corners_ij(subzone)['ur'])
if psproj == 'npstere':
boundinglat = boundinglat.min()
elif psproj == 'spstere':
boundinglat = boundinglat.max()
# 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':
if self.secant_projection:
lat = 'secant_lat'
else:
lat = 'reference_lat'
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':
if self.secant_projection:
lat = 'secant_lat'
else:
lat = 'reference_lat'
b = Basemap(resolution=gisquality, projection=psproj,
lat_ts=self.projection[lat].get('degrees'),
lat_0=self.projection['center_lat'].get('degrees'),
lon_0=self.projection['reference_lon'].get('degrees'),
boundinglat=boundinglat)
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 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]class H2DGaussGeometry(H2DGeometry):
"""
Handles the geometry for a 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):
"""
(*i, j*) being the coordinates in the 2D matrix of gridpoints, \n
(*lon, lat*) being the lon/lat coordinates in degrees.
Computation adapted from arpifs/transform/tracare.F90
"""
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):
"""
(*lon, lat*) being the lon/lat coordinates in degrees, \n
(*i, j*) being the coordinates in the 2D matrix of gridpoints.
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)
latitudes = [self.grid['latitudes'][n].get('radians')
for n in range(0, self.dimensions['lat_number'])]
if isinstance(lat, numpy.ndarray):
j = numpy.zeros(len(lat), dtype=numpy.int)
for k in range(0, len(lat)):
distmin = 999
nearest = None
for n in range(0, len(latitudes)):
dist = abs(lat[k]-latitudes[n])
if dist < distmin:
distmin = dist
nearest = n
j[k] = nearest
i = numpy.rint([lon[k] * self.dimensions['lon_number_by_lat'][j[k]] / (numpy.pi * 2) for k in range(len(lon))])
else:
distmin = 999
nearest = None
for n in range(0, len(latitudes)):
dist = abs(lat-latitudes[n])
if dist < distmin:
distmin = dist
nearest = n
j = nearest
i = numpy.rint(lon * self.dimensions['lon_number_by_lat'][j] / (numpy.pi * 2))
return (i, j)
[docs] def get_lonlat_grid(self, **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).
"""
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)
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*.
"""
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 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