Browse Source

Now uses CoordPair for vertical cross sections so that lat/lon points can be used directly.

main
Bill Ladwig 9 years ago
parent
commit
9be5a86a1b
  1. 225
      src/wrf/interp.py
  2. 116
      src/wrf/interputils.py
  3. 33
      src/wrf/latlon.py
  4. 23
      src/wrf/latlonutils.py
  5. 86
      src/wrf/metadecorators.py
  6. 119
      src/wrf/projection.py
  7. 18
      src/wrf/projutils.py
  8. 2
      src/wrf/py3compat.py
  9. 10
      src/wrf/util.py

225
src/wrf/interp.py

@ -10,7 +10,7 @@ from .extension import (_interpz3d, _vertcross, _interpline, _smooth2d,
from .metadecorators import set_interp_metadata from .metadecorators import set_interp_metadata
from .util import extract_vars, is_staggered, get_id, npvalues from .util import extract_vars, is_staggered, get_id, npvalues
from .py3compat import py3range from .py3compat import py3range
from .interputils import get_xy, get_xy_z_params from .interputils import get_xy, get_xy_z_params, to_xy_coords
from .constants import Constants, ConversionFactors from .constants import Constants, ConversionFactors
from .terrain import get_terrain from .terrain import get_terrain
from .geoht import get_height from .geoht import get_height
@ -91,6 +91,7 @@ def interplevel(field3d, vert, desiredlev, missing=Constants.DEFAULT_FILL,
@set_interp_metadata("cross") @set_interp_metadata("cross")
def vertcross(field3d, vert, levels=None, missing=Constants.DEFAULT_FILL, def vertcross(field3d, vert, levels=None, missing=Constants.DEFAULT_FILL,
wrfin=None, timeidx=0, stagger=None, projection=None,
pivot_point=None, angle=None, pivot_point=None, angle=None,
start_point=None, end_point=None, start_point=None, end_point=None,
latlon=False, cache=None, meta=True): latlon=False, cache=None, meta=True):
@ -99,7 +100,10 @@ def vertcross(field3d, vert, levels=None, missing=Constants.DEFAULT_FILL,
The cross section is defined by a horizontal line through the domain. The cross section is defined by a horizontal line through the domain.
This horizontal line is defined by either including the This horizontal line is defined by either including the
*pivot_point* and *angle* parameters, or the *start_point* and *pivot_point* and *angle* parameters, or the *start_point* and
*end_point* parameters. *end_point* parameters. The *pivot_point*, *start_point*, and *end_point*
coordinates can be defined in either x,y or latitude,longitude space.
If latitude,longitude coordinates are used, then a WRF input file or
map projection must also be specified.
The vertical levels for the cross section are fixed if *levels* is not The vertical levels for the cross section are fixed if *levels* is not
specified, and are determined by dividing the vertical coordinate in to specified, and are determined by dividing the vertical coordinate in to
@ -124,32 +128,71 @@ def vertcross(field3d, vert, levels=None, missing=Constants.DEFAULT_FILL,
as *field3d* as *field3d*
levels (sequence, optional): A sequence of :obj:`float` for the desired levels (sequence, optional): A sequence of :obj:`float` for the desired
vertical levels in the output array. If None, a fixed set of vertical levels in the output array. Must be in the same units
vertical levels is provided. Default is None. as *vert*. If None, a fixed set of vertical levels is provided.
Default is None.
missing (:obj:`float`): The fill value to use for the output. missing (:obj:`float`): The fill value to use for the output.
Default is :data:`wrf.Constants.DEFAULT_FILL`. Default is :data:`wrf.Constants.DEFAULT_FILL`.
pivot_point (:obj:`tuple` or :obj:`list`, optional): A wrfin (:class:`netCDF4.Dataset`, :class:`Nio.NioFile`, or an \
:obj:`tuple` or :obj:`list` with two entries, iterable, optional): Input WRF ARW NetCDF
in the form of [x, y] (or [west_east, south_north]), which data as a :class:`netCDF4.Dataset`, :class:`Nio.NioFile`
indicates the x,y location through which the plane will pass. or an iterable sequence of the aforementioned types. This is used
Must also specify `angle`. to obtain the map projection when using latitude,longitude
coordinates. Default is None.
timeidx (:obj:`int`, optional): The
desired time index when obtaining map boundary information
from moving nests. This value can be a positive or negative integer.
Only required when *wrfin* is specified and the nest is moving.
Currently, :data:`wrf.ALL_TIMES` is not supported.
Default is 0.
stagger (:obj:`str`): If using latitude, longitude coordinate pairs
for *start_point*, *end_point*, or *pivot_point*,
set the appropriate grid staggering type for *field2d*. By default,
the mass grid is used. The options are:
- 'm': Use the mass grid (default).
- 'u': Use the same staggered grid as the u wind component,
which has a staggered west_east (x) dimension.
- 'v': Use the same staggered grid as the v wind component,
which has a staggered south_north (y) dimension.
projection (:class:`wrf.WrfProj` subclass, optional): The map
projection object to use when working with latitude, longitude
coordinates, and must be specified if *wrfin* is None. Default
is None.
pivot_point (:class`wrf.CoordPair`, optional): A coordinate pair for
the pivot point, which indicates the location through which
the plane will pass. Must also specify *angle*. The coordinate
pair can be in x,y grid coordinates or latitude, longitude
coordinates. If using latitude, longitude coordinates, then
either *wrfin* or *projection* must be specified to obtain the
map projection. Default is None.
angle (:obj:`float`, optional): Only valid for cross sections where angle (:obj:`float`, optional): Only valid for cross sections where
a plane will be plotted through a plane will be plotted through
a given point on the model domain. 0.0 represents a S-N cross a given point on the model domain. 0.0 represents a S-N cross
section. 90.0 is a W-E cross section. section. 90.0 is a W-E cross section.
start_point (:obj:`tuple` or :obj:`list`, optional): A start_point (:class:`wrf.CoordPair`, optional): A coordinate pair
:obj:`tuple` or :obj:`list` with two entries, in the form of which indicates the start location through which the plane will
[x, y] (or [west_east, south_north]), which indicates the start pass. Must also specify *end_point*. The coordinate
x,y location through which the plane will pass. pair can be in x,y grid coordinates or latitude, longitude
coordinates. If using latitude, longitude coordinates, then
end_point (:obj:`tuple` or :obj:`list`, optional): A either *wrfin* or *projection* must be specified to obtain the
:obj:`tuple` or :obj:`list` with two entries, in the form of map projection. Default is None.
[x, y] (or [west_east, south_north]), which indicates the end x,y
location through which the plane will pass. end_point (:class:`wrf.CoordPair`, optional): A coordinate pair
which indicates the end location through which the plane will
pass. Must also specify *end_point*. The coordinate
pair can be in x,y grid coordinates or latitude, longitude
coordinates. If using latitude, longitude coordinates, then
either *wrfin* or *projection* must be specified to obtain the
map projection. Default is None.
latlon (:obj:`bool`, optional): Set to True to also interpolate the latlon (:obj:`bool`, optional): Set to True to also interpolate the
two-dimensional latitude and longitude coordinates along the same two-dimensional latitude and longitude coordinates along the same
@ -178,6 +221,10 @@ def vertcross(field3d, vert, levels=None, missing=Constants.DEFAULT_FILL,
:class:`numpy.ndarray` object with no metadata. :class:`numpy.ndarray` object with no metadata.
""" """
if timeidx is None:
raise ValueError("'timeidx' must be a positive or negative integer")
# Some fields like uvmet have an extra left dimension for the product # Some fields like uvmet have an extra left dimension for the product
# type, we'll handle that iteration here. # type, we'll handle that iteration here.
multi = True if field3d.ndim - vert.ndim == 1 else False multi = True if field3d.ndim - vert.ndim == 1 else False
@ -187,9 +234,37 @@ def vertcross(field3d, vert, levels=None, missing=Constants.DEFAULT_FILL,
var2dz = cache["var2dz"] var2dz = cache["var2dz"]
z_var2d = cache["z_var2d"] z_var2d = cache["z_var2d"]
except (KeyError, TypeError): except (KeyError, TypeError):
xy, var2dz, z_var2d = get_xy_z_params(npvalues(vert), pivot_point, # Convert Lat/Lon points to grid points
angle, start_point, end_point, start_point_xy = None
levels) end_point_xy = None
pivot_point_xy = None
if pivot_point is not None:
if pivot_point.lat is not None and pivot_point.lon is not None:
xy_coords = to_xy_coords(pivot_point, wrfin, timeidx,
stagger, projection)
pivot_point_xy = (xy_coords.x, xy_coords.y)
else:
pivot_point_xy = (pivot_point.x, pivot_point.y)
if start_point is not None and end_point is not None:
if start_point.lat is not None and start_point.lon is not None:
xy_coords = to_xy_coords(start_point, wrfin, timeidx,
stagger, projection)
start_point_xy = (xy_coords.x, xy_coords.y)
else:
start_point_xy = (start_point.x, start_point.y)
if end_point.lat is not None and end_point.lon is not None:
xy_coords = to_xy_coords(end_point, wrfin, timeidx,
stagger, projection)
end_point_xy = (xy_coords.x, xy_coords.y)
else:
end_point_xy = (end_point.x, end_point.y)
xy, var2dz, z_var2d = get_xy_z_params(npvalues(vert), pivot_point_xy,
angle, start_point_xy,
end_point_xy, levels)
if not multi: if not multi:
result = _vertcross(field3d, xy, var2dz, z_var2d, missing) result = _vertcross(field3d, xy, var2dz, z_var2d, missing)
@ -206,36 +281,84 @@ def vertcross(field3d, vert, levels=None, missing=Constants.DEFAULT_FILL,
@set_interp_metadata("line") @set_interp_metadata("line")
def interpline(field2d, pivot_point=None, def interpline(field2d, pivot_point=None,
wrfin=None, timeidx=0, stagger=None, projection=None,
angle=None, start_point=None, angle=None, start_point=None,
end_point=None, latlon=False, end_point=None, latlon=False,
cache=None, meta=True): cache=None, meta=True):
"""Return the two-dimensional field interpolated along a line. """Return the two-dimensional field interpolated along a line.
This line is defined by either including the
*pivot_point* and *angle* parameters, or the *start_point* and
*end_point* parameters. The *pivot_point*, *start_point*, and *end_point*
coordinates can be defined in either x,y or latitude,longitude space.
If latitude,longitude coordinates are used, then a WRF input file or
map projection must also be specified.
Args: Args:
field2d (:class:`xarray.DataArray` or :class:`numpy.ndarray`): field2d (:class:`xarray.DataArray` or :class:`numpy.ndarray`):
A two-dimensional field. A two-dimensional field.
pivot_point (:obj:`tuple` or :obj:`list`, optional): A wrfin (:class:`netCDF4.Dataset`, :class:`Nio.NioFile`, or an \
:obj:`tuple` or :obj:`list` with two entries, iterable, optional): Input WRF ARW NetCDF
in the form of [x, y] (or [west_east, south_north]), which data as a :class:`netCDF4.Dataset`, :class:`Nio.NioFile`
indicates the x,y location through which the plane will pass. or an iterable sequence of the aforementioned types. This is used
Must also specify `angle`. to obtain the map projection when using latitude,longitude
coordinates. Should not be used when working with x,y
coordinates. Default is None.
timeidx (:obj:`int`, optional): The
desired time index when obtaining map boundary information
from moving nests. This value can be a positive or negative integer.
Only required when *wrfin* is specified and the nest is moving.
Currently, :data:`wrf.ALL_TIMES` is not supported.
Default is 0.
stagger (:obj:`str`): If using latitude, longitude coordinate pairs
for *start_point*, *end_point*, or *pivot_point*,
set the appropriate grid staggering type for *field2d*. By default,
the mass grid is used. The options are:
- 'm': Use the mass grid (default).
- 'u': Use the same staggered grid as the u wind component,
which has a staggered west_east (x) dimension.
- 'v': Use the same staggered grid as the v wind component,
which has a staggered south_north (y) dimension.
projection (:class:`wrf.WrfProj`, optional): The map
projection object to use when working with latitude, longitude
coordinates, and must be specified if *wrfin* is None. Should
not be used when working with x,y coordinates. Default
is None.
pivot_point (:class`wrf.CoordPair`, optional): A coordinate pair for
the pivot point, which indicates the location through which
the plane will pass. Must also specify *angle*. The coordinate
pair can be in x,y grid coordinates or latitude, longitude
coordinates. If using latitude, longitude coordinates, then
either *wrfin* or *projection* must be specified to obtain the
map projection. Default is None.
angle (:obj:`float`, optional): Only valid for cross sections where angle (:obj:`float`, optional): Only valid for cross sections where
a plane will be plotted through a plane will be plotted through
a given point on the model domain. 0.0 represents a S-N cross a given point on the model domain. 0.0 represents a S-N cross
section. 90.0 is a W-E cross section. section. 90.0 is a W-E cross section.
start_point (:obj:`tuple` or :obj:`list`, optional): A start_point (:class:`wrf.CoordPair`, optional): A coordinate pair
:obj:`tuple` or :obj:`list` with two entries, in the form of which indicates the start location through which the plane will
[x, y] (or [west_east, south_north]), which indicates the start pass. Must also specify *end_point*. The coordinate
x,y location through which the plane will pass. pair can be in x,y grid coordinates or latitude, longitude
coordinates. If using latitude, longitude coordinates, then
end_point (:obj:`tuple` or :obj:`list`, optional): A either *wrfin* or *projection* must be specified to obtain the
:obj:`tuple` or :obj:`list` with two entries, in the form of map projection. Default is None.
[x, y] (or [west_east, south_north]), which indicates the end x,y
location through which the plane will pass. end_point (:class:`wrf.CoordPair`, optional): A coordinate pair
which indicates the end location through which the plane will
pass. Must also specify *end_point*. The coordinate
pair can be in x,y grid coordinates or latitude, longitude
coordinates. If using latitude, longitude coordinates, then
either *wrfin* or *projection* must be specified to obtain the
map projection. Default is None.
latlon (:obj:`bool`, optional): Set to True to also interpolate the latlon (:obj:`bool`, optional): Set to True to also interpolate the
two-dimensional latitude and longitude coordinates along the same two-dimensional latitude and longitude coordinates along the same
@ -266,11 +389,41 @@ def interpline(field2d, pivot_point=None,
""" """
if timeidx is None:
raise ValueError("'timeidx' must be a positive or negative integer")
try: try:
xy = cache["xy"] xy = cache["xy"]
except (KeyError, TypeError): except (KeyError, TypeError):
xy = get_xy(field2d, pivot_point, angle, start_point, end_point) start_point_xy = None
end_point_xy = None
pivot_point_xy = None
if pivot_point is not None:
if pivot_point.lat is not None and pivot_point.lon is not None:
xy_coords = to_xy_coords(pivot_point, wrfin, timeidx,
stagger, projection)
pivot_point_xy = (xy_coords.x, xy_coords.y)
else:
pivot_point_xy = (pivot_point.x, pivot_point.y)
if start_point is not None and end_point is not None:
if start_point.lat is not None and start_point.lon is not None:
xy_coords = to_xy_coords(start_point, wrfin, timeidx,
stagger, projection)
start_point_xy = (xy_coords.x, xy_coords.y)
else:
start_point_xy = (start_point.x, start_point.y)
if end_point.lat is not None and end_point.lon is not None:
xy_coords = to_xy_coords(end_point, wrfin, timeidx,
stagger, projection)
end_point_xy = (xy_coords.x, xy_coords.y)
else:
end_point_xy = (end_point.x, end_point.y)
xy = get_xy(field2d, pivot_point_xy, angle, start_point_xy,
end_point_xy)
return _interpline(field2d, xy) return _interpline(field2d, xy)

116
src/wrf/interputils.py

@ -2,11 +2,15 @@ from __future__ import (absolute_import, division, print_function,
unicode_literals) unicode_literals)
from math import floor, ceil from math import floor, ceil
from collections import Iterable
import numpy as np import numpy as np
from .extension import _interp2dxy from .extension import _interp2dxy
from .py3compat import py3range from .py3compat import py3range
from .coordpair import CoordPair
from .constants import Constants, ProjectionTypes
from .latlonutils import _ll_to_xy
def to_positive_idxs(shape, coord): def to_positive_idxs(shape, coord):
@ -256,11 +260,11 @@ def get_xy_z_params(z, pivot_point=None, angle=None,
nlevels = int(z_max/dz) nlevels = int(z_max/dz)
z_var2d = np.zeros((nlevels), dtype=z.dtype) z_var2d = np.zeros((nlevels), dtype=z.dtype)
z_var2d[0] = z_min z_var2d[0] = z_min
else:
z_var2d = np.asarray(levels)
for i in py3range(1,nlevels): for i in py3range(1,nlevels):
z_var2d[i] = z_var2d[0] + i*dz z_var2d[i] = z_var2d[0] + i*dz
else:
z_var2d = np.asarray(levels, z.dtype)
return xy, var2dz, z_var2d return xy, var2dz, z_var2d
@ -323,3 +327,107 @@ def get_xy(var, pivot_point=None, angle=None,
xy = _calc_xy(xdim, ydim, pos_pivot, angle, pos_start, pos_end) xy = _calc_xy(xdim, ydim, pos_pivot, angle, pos_start, pos_end)
return xy return xy
def to_xy_coords(pairs, wrfin=None, timeidx=0, stagger=None, projection=None):
"""Return the coordinate pairs in grid space.
This function converts latitude,longitude coordinate pairs to
x,y coordinate pairs.
Args:
pairs (:class:`CoordPair` or sequence): A single coordinate pair or
a sequence of coordinate pairs to be converted.
wrfin (:class:`netCDF4.Dataset`, :class:`Nio.NioFile`, or an \
iterable, optional): Input WRF ARW NetCDF
data as a :class:`netCDF4.Dataset`, :class:`Nio.NioFile`
or an iterable sequence of the aforementioned types. This is used
to obtain the map projection when using latitude,longitude
coordinates. Should not be used when working with x,y
coordinates. Default is None.
timeidx (:obj:`int` or :data:`wrf.ALL_TIMES`, optional): The
desired time index when obtaining map boundary information
from moving nests. This value can be a positive integer,
negative integer, or
:data:`wrf.ALL_TIMES` (an alias for None) to return
all times in the file or sequence. Only required when
*wrfin* is specified and the nest is moving. Default is 0.
stagger (:obj:`str`): If using latitude, longitude coordinate pairs
for *start_point*, *end_point*, or *pivot_point*,
set the appropriate grid staggering type for *field2d*. By default,
the mass grid is used. The options are:
- 'm': Use the mass grid (default).
- 'u': Use the same staggered grid as the u wind component,
which has a staggered west_east (x) dimension.
- 'v': Use the same staggered grid as the v wind component,
which has a staggered south_north (y) dimension.
projection (:class:`wrf.WrfProj`, optional): The map
projection object to use when working with latitude, longitude
coordinates, and must be specified if *wrfin* is None. Default
is None.
Returns:
:class:`wrf.CoordPair` or sequence: The coordinate pair(s) in
x,y grid coordinates.
"""
if wrfin is None and projection is None:
raise ValueError ("'wrfin' or 'projection' parameter is required")
if isinstance(pairs, Iterable):
lat = [pair.lat for pair in pairs]
lon = [pair.lon for pair in pairs]
else:
lat = pairs.lat
lon = pairs.lon
if wrfin is not None:
xy_vals = _ll_to_xy(lat, lon, wrfin=wrfin, timeidx=timeidx,
squeeze=True, meta=False, stagger=stagger, as_int=True)
else:
map_proj = projection.map_proj
if map_proj == ProjectionTypes.LAT_LON:
pole_lat = projection.pole_lat
pole_lon = projection.pole_lon
latinc = ((projection.dy*360.0)/2.0 /
Constants.PI/Constants.WRF_EARTH_RADIUS)
loninc = ((projection.dx*360.0)/2.0 /
Constants.PI/Constants.WRF_EARTH_RADIUS)
else:
pole_lat = 90.0
pole_lon = 0.0
latinc = 0.0
loninc = 0.0
xy_vals = _ll_to_xy(lat, lon, meta=False, squeeze=True,
as_int=True,
map_proj=projection.map_proj,
truelat1=projection.truelat1,
truelat2=projection.truelat2,
stand_lon=projection.stand_lon,
ref_lat=projection.ll_lat,
ref_lon=projection.ll_lon,
pole_lat=pole_lat,
pole_lon=pole_lon,
known_x=0,
known_y=0,
dx=projection.dx,
dy=projection.dy,
latinc=latinc,
loninc=loninc)
if xy_vals.ndim == 1:
return CoordPair(x=xy_vals[0], y=xy_vals[1])
else:
return [CoordPair(x=xy_vals[0,i], y=xy_vals[1,i])
for i in py3range(xy_vals.shape[1])]

33
src/wrf/latlon.py

@ -290,11 +290,17 @@ def ll_to_xy_proj(latitude, longitude, meta=True, squeeze=True, as_int=True,
dy (:obj:`float`) - The y spacing in meters at the true latitude. dy (:obj:`float`) - The y spacing in meters at the true latitude.
Required for *map_proj* = 1, 2, 3 (defaults to 0 otherwise). Required for *map_proj* = 1, 2, 3 (defaults to 0 otherwise).
latinc (:obj:`float`): The lower left corner latitude at (0,0). latinc (:obj:`float`): Required for *map_proj* = 6. Defined as:
Required for *map_proj* = 6.
loninc (:obj:`float`) - The lower left corner longitude at (0,0). .. code-block:: python
Required for *map_proj* = 6.
latinc = (dy*360.0)/2.0/Constants.PI/Constants.WRF_EARTH_RADIUS
loninc (:obj:`float`): Required for *map_proj* = 6. Defined as:
.. code-block:: python
loninc = (dx*360.0)/2.0/Constants.PI/Constants.WRF_EARTH_RADIUS
Returns: Returns:
:class:`xarray.DataArray` or :class:`numpy.ndarray`: The :class:`xarray.DataArray` or :class:`numpy.ndarray`: The
@ -450,11 +456,17 @@ def xy_to_ll_proj(x, y, meta=True, squeeze=True, map_proj=None, truelat1=None,
dy (:obj:`float`) - The y spacing in meters at the true latitude. dy (:obj:`float`) - The y spacing in meters at the true latitude.
Required for *map_proj* = 1, 2, 3 (defaults to 0 otherwise). Required for *map_proj* = 1, 2, 3 (defaults to 0 otherwise).
latinc (:obj:`float`): The lower left corner latitude at (0,0). latinc (:obj:`float`): Required for *map_proj* = 6. Defined as:
Required for *map_proj* = 6.
.. code-block:: python
latinc = (dy*360.0)/2.0/Constants.PI/Constants.WRF_EARTH_RADIUS
loninc (:obj:`float`): Required for *map_proj* = 6. Defined as:
loninc (:obj:`float`) - The lower left corner longitude at (0,0). .. code-block:: python
Required for *map_proj* = 6.
loninc = (dx*360.0)/2.0/Constants.PI/Constants.WRF_EARTH_RADIUS
Returns: Returns:
:class:`xarray.DataArray` or :class:`numpy.ndarray`: The :class:`xarray.DataArray` or :class:`numpy.ndarray`: The
@ -475,3 +487,8 @@ def xy_to_ll_proj(x, y, meta=True, squeeze=True, map_proj=None, truelat1=None,

23
src/wrf/latlonutils.py

@ -11,6 +11,8 @@ from .util import (extract_vars, extract_global_attrs,
either, is_moving_domain, is_multi_time_req, either, is_moving_domain, is_multi_time_req,
iter_left_indexes, is_mapping, is_multi_file) iter_left_indexes, is_mapping, is_multi_file)
from .py3compat import viewkeys, viewitems from .py3compat import viewkeys, viewitems
from .projutils import dict_keys_to_upper
def _lat_varname(wrfin, stagger): def _lat_varname(wrfin, stagger):
"""Return the latitude variable name for the specified stagger type. """Return the latitude variable name for the specified stagger type.
@ -194,20 +196,6 @@ def _get_proj_params(wrfin, timeidx, stagger, method, squeeze, cache, _key):
pole_lat, pole_lon, known_x, known_y, dx, dy, latinc, loninc) pole_lat, pole_lon, known_x, known_y, dx, dy, latinc, loninc)
def _dict_keys_to_upper(d):
"""Return a dictionary with the keys changed to uppercase.
Args:
d (:obj:`dict`): A dictionary.
Returns:
:obj:`dict`: A dictionary with uppercase keys.
"""
return {key.upper(): val for key, val in viewitems(d)}
# known_x and known_y are 0-based # known_x and known_y are 0-based
def _kwarg_proj_params(**projparams): def _kwarg_proj_params(**projparams):
"""Return the map projection parameters. """Return the map projection parameters.
@ -224,7 +212,7 @@ def _kwarg_proj_params(**projparams):
:obj:`tuple`: The map projection parameters. :obj:`tuple`: The map projection parameters.
""" """
projparams = _dict_keys_to_upper(projparams) projparams = dict_keys_to_upper(projparams)
map_proj = projparams.get("MAP_PROJ") map_proj = projparams.get("MAP_PROJ")
truelat1 = projparams.get("TRUELAT1") truelat1 = projparams.get("TRUELAT1")
@ -609,3 +597,8 @@ def _xy_to_ll(x, y, wrfin=None, timeidx=0, stagger=None,
loninc, x_val, y_val) loninc, x_val, y_val)
return result return result

86
src/wrf/metadecorators.py

@ -12,7 +12,7 @@ from .util import (extract_vars, either, from_args, arg_location,
from_var, iter_left_indexes) from_var, iter_left_indexes)
from .coordpair import CoordPair from .coordpair import CoordPair
from .py3compat import viewkeys, viewitems, py3range, ucode from .py3compat import viewkeys, viewitems, py3range, ucode
from .interputils import get_xy_z_params, get_xy from .interputils import get_xy_z_params, get_xy, to_xy_coords
from .config import xarray_enabled from .config import xarray_enabled
if xarray_enabled(): if xarray_enabled():
@ -830,7 +830,9 @@ def _set_cross_meta(wrapped, instance, args, kwargs):
:mod:`wrapt` :mod:`wrapt`
""" """
argvars = from_args(wrapped, ("field3d", "vert", "latlon", "missing", argvars = from_args(wrapped, ("field3d", "vert", "levels",
"latlon", "missing",
"wrfin", "timeidx", "stagger", "projection",
"pivot_point", "angle", "pivot_point", "angle",
"start_point", "end_point", "start_point", "end_point",
"cache"), "cache"),
@ -838,16 +840,48 @@ def _set_cross_meta(wrapped, instance, args, kwargs):
field3d = argvars["field3d"] field3d = argvars["field3d"]
z = argvars["vert"] z = argvars["vert"]
levels = argvars["levels"]
inc_latlon = argvars["latlon"] inc_latlon = argvars["latlon"]
missingval = argvars["missing"] missingval = argvars["missing"]
wrfin = argvars["wrfin"]
timeidx = argvars["timeidx"]
stagger = argvars["stagger"]
projection = argvars["projection"]
pivot_point = argvars["pivot_point"] pivot_point = argvars["pivot_point"]
angle = argvars["angle"] angle = argvars["angle"]
start_point = argvars["start_point"] start_point = argvars["start_point"]
end_point = argvars["end_point"] end_point = argvars["end_point"]
cache = argvars["cache"] cache = argvars["cache"]
xy, var2dz, z_var2d = get_xy_z_params(npvalues(z), pivot_point, angle, start_point_xy = None
start_point, end_point) end_point_xy = None
pivot_point_xy = None
if pivot_point is not None:
if pivot_point.lat is not None and pivot_point.lon is not None:
xy_coords = to_xy_coords(pivot_point, wrfin, timeidx,
stagger, projection)
pivot_point_xy = (xy_coords.x, xy_coords.y)
else:
pivot_point_xy = (pivot_point.x, pivot_point.y)
if start_point is not None and end_point is not None:
if start_point.lat is not None and start_point.lon is not None:
xy_coords = to_xy_coords(start_point, wrfin, timeidx,
stagger, projection)
start_point_xy = (xy_coords.x, xy_coords.y)
else:
start_point_xy = (start_point.x, start_point.y)
if end_point.lat is not None and end_point.lon is not None:
xy_coords = to_xy_coords(end_point, wrfin, timeidx,
stagger, projection)
end_point_xy = (xy_coords.x, xy_coords.y)
else:
end_point_xy = (end_point.x, end_point.y)
xy, var2dz, z_var2d = get_xy_z_params(npvalues(z), pivot_point_xy, angle,
start_point_xy, end_point_xy, levels)
# Make a copy so we don't modify a user supplied cache # Make a copy so we don't modify a user supplied cache
if cache is not None: if cache is not None:
@ -877,8 +911,7 @@ def _set_cross_meta(wrapped, instance, args, kwargs):
ed_x = xy[-1,0] ed_x = xy[-1,0]
ed_y = xy[-1,1] ed_y = xy[-1,1]
cross_str = "({0}, {1}) to ({2}, {3})".format(st_x, st_y, cross_str = "({0}, {1}) to ({2}, {3})".format(st_x, st_y, ed_x, ed_y)
ed_x, ed_y)
if angle is not None: if angle is not None:
cross_str += " ; center={0} ; angle={1}".format(pivot_point, cross_str += " ; center={0} ; angle={1}".format(pivot_point,
angle) angle)
@ -1009,12 +1042,18 @@ def _set_line_meta(wrapped, instance, args, kwargs):
:mod:`wrapt` :mod:`wrapt`
""" """
argvars = from_args(wrapped, ("field2d", "pivot_point", "angle", argvars = from_args(wrapped, ("field2d",
"wrfin", "timeidx", "stagger", "projection",
"pivot_point", "angle",
"start_point", "end_point", "latlon", "start_point", "end_point", "latlon",
"cache"), "cache"),
*args, **kwargs) *args, **kwargs)
field2d = argvars["field2d"] field2d = argvars["field2d"]
wrfin = argvars["wrfin"]
timeidx = argvars["timeidx"]
stagger = argvars["stagger"]
projection = argvars["projection"]
pivot_point = argvars["pivot_point"] pivot_point = argvars["pivot_point"]
angle = argvars["angle"] angle = argvars["angle"]
start_point = argvars["start_point"] start_point = argvars["start_point"]
@ -1025,7 +1064,35 @@ def _set_line_meta(wrapped, instance, args, kwargs):
if cache is None: if cache is None:
cache = {} cache = {}
xy = get_xy(field2d, pivot_point, angle, start_point, end_point) start_point_xy = None
end_point_xy = None
pivot_point_xy = None
if pivot_point is not None:
if pivot_point.lat is not None and pivot_point.lon is not None:
xy_coords = to_xy_coords(pivot_point, wrfin, timeidx,
stagger, projection)
pivot_point_xy = (xy_coords.x, xy_coords.y)
else:
pivot_point_xy = (pivot_point.x, pivot_point.y)
if start_point is not None and end_point is not None:
if start_point.lat is not None and start_point.lon is not None:
xy_coords = to_xy_coords(start_point, wrfin, timeidx,
stagger, projection)
start_point_xy = (xy_coords.x, xy_coords.y)
else:
start_point_xy = (start_point.x, start_point.y)
if end_point.lat is not None and end_point.lon is not None:
xy_coords = to_xy_coords(end_point, wrfin, timeidx,
stagger, projection)
end_point_xy = (xy_coords.x, xy_coords.y)
else:
end_point_xy = (end_point.x, end_point.y)
xy = get_xy(field2d, pivot_point_xy, angle, start_point_xy, end_point_xy)
# Make a copy so we don't modify a user supplied cache # Make a copy so we don't modify a user supplied cache
new_cache = dict(cache) new_cache = dict(cache)
@ -1050,8 +1117,7 @@ def _set_line_meta(wrapped, instance, args, kwargs):
ed_x = xy[-1,0] ed_x = xy[-1,0]
ed_y = xy[-1,1] ed_y = xy[-1,1]
cross_str = "({0}, {1}) to ({2}, {3})".format(st_x, st_y, cross_str = "({0}, {1}) to ({2}, {3})".format(st_x, st_y, ed_x, ed_y)
ed_x, ed_y)
if angle is not None: if angle is not None:
cross_str += " ; center={0} ; angle={1}".format(pivot_point, cross_str += " ; center={0} ; angle={1}".format(pivot_point,
angle) angle)

119
src/wrf/projection.py

@ -3,8 +3,10 @@ from __future__ import (absolute_import, division, print_function,
import numpy as np import numpy as np
import math import math
from .coordpair import CoordPair
from .config import basemap_enabled, cartopy_enabled, pyngl_enabled from .config import basemap_enabled, cartopy_enabled, pyngl_enabled
from .constants import Constants, ProjectionTypes from .constants import Constants, ProjectionTypes
from .projutils import dict_keys_to_upper
if cartopy_enabled(): if cartopy_enabled():
from cartopy import crs from cartopy import crs
@ -78,6 +80,7 @@ if cartopy_enabled():
self._threshold = np.diff(self.x_limits)[0] / 720 self._threshold = np.diff(self.x_limits)[0] / 720
def _ismissing(val): def _ismissing(val):
"""Return True if a value is None, greater than 90.0, or less than -90. """Return True if a value is None, greater than 90.0, or less than -90.
@ -95,6 +98,7 @@ def _ismissing(val):
""" """
return val is None or val > 90. or val < -90. return val is None or val > 90. or val < -90.
class WrfProj(object): class WrfProj(object):
"""A base class for storing map projection information from WRF data. """A base class for storing map projection information from WRF data.
@ -117,6 +121,25 @@ class WrfProj(object):
top_right (indexable sequence): A pair of (ur_lat, ur_lon). top_right (indexable sequence): A pair of (ur_lat, ur_lon).
map_proj (:obj:`int`): Model projection integer id.
truelat1 (:obj:`float`): True latitude 1.
truelat2 (:obj:`float`): True latitude 2.
moad_cen_lat (:obj:`float`): Mother of all domains center latitude.
stand_lon (:obj:`float`): Standard longitude.
pole_lat (:obj:`float`): The pole latitude.
pole_lon (:obj:`float`): The pole longitude.
dx (:obj:`float`): The x grid spacing.
dy (:obj:`float`): The y grid spacing.
""" """
def __init__(self, bottom_left=None, top_right=None, def __init__(self, bottom_left=None, top_right=None,
lats=None, lons=None, **proj_params): lats=None, lons=None, **proj_params):
@ -124,13 +147,13 @@ class WrfProj(object):
Args: Args:
bottom_left (indexable sequence, optional): The lower left corner bottom_left (:class:`wrf.CoordPair`, optional): The lower left
as a (latitude, longitude) pair. Must also specify *top_right* corner. Must also specify *top_right* if used.
if used. Default is None. Default is None.
top_right (indexable sequence): The upper right corner as a top_right (:class:`wrf.CoordPair`, optional): The upper right
(latitude, longitude) pair. Must also specify *bottom_left* corner. Must also specify *bottom_left* if used.
if used. Default is None. Default is None.
lats (:class:`numpy.ndarray`, optional): An array of at least lats (:class:`numpy.ndarray`, optional): An array of at least
two dimensions containing all of the latitude values. Must two dimensions containing all of the latitude values. Must
@ -142,11 +165,9 @@ class WrfProj(object):
**proj_params: Map projection optional keyword arguments, that **proj_params: Map projection optional keyword arguments, that
have the same names as found in WRF output NetCDF global have the same names as found in WRF output NetCDF global
attributes: attributes (case insensitive):
- 'MAP_PROJ': The map projection type as an integer. - 'MAP_PROJ': The map projection type as an integer.
- 'CEN_LAT': Center latitude.
- 'CEN_LON': Center longitude.
- 'TRUELAT1': True latitude 1. - 'TRUELAT1': True latitude 1.
- 'TRUELAT2': True latitude 2. - 'TRUELAT2': True latitude 2.
- 'MOAD_CEN_LAT': Mother of all domains center latitude. - 'MOAD_CEN_LAT': Mother of all domains center latitude.
@ -157,10 +178,10 @@ class WrfProj(object):
""" """
if bottom_left is not None and top_right is not None: if bottom_left is not None and top_right is not None:
self.ll_lat = bottom_left[0] self.ll_lat = bottom_left.lat
self.ll_lon = bottom_left[1] self.ll_lon = bottom_left.lon
self.ur_lat = top_right[0] self.ur_lat = top_right.lat
self.ur_lon = top_right[1] self.ur_lon = top_right.lon
self.bottom_left = bottom_left self.bottom_left = bottom_left
self.top_right = top_right self.top_right = top_right
elif lats is not None and lons is not None: elif lats is not None and lons is not None:
@ -168,24 +189,33 @@ class WrfProj(object):
self.ur_lat = lats[-1,-1] self.ur_lat = lats[-1,-1]
self.ll_lon = lons[0,0] self.ll_lon = lons[0,0]
self.ur_lon = lons[-1,-1] self.ur_lon = lons[-1,-1]
self.bottom_left = [self.ll_lat, self.ll_lon] self.bottom_left = CoordPair(lat=self.ll_lat, lon=self.ll_lon)
self.top_right = [self.ur_lat, self.ur_lon] self.top_right = CoordPair(self.ur_lat, self.ur_lon)
else: else:
raise ValueError("invalid corner point arguments") raise ValueError("invalid corner point arguments")
up_proj_params = dict_keys_to_upper(proj_params)
self.map_proj = up_proj_params.get("MAP_PROJ", None)
# These indicate the center of the nest/domain, not necessarily the # These indicate the center of the nest/domain, not necessarily the
# center of the projection # center of the projection
self._cen_lat = proj_params.get("CEN_LAT", None) self._cen_lat = up_proj_params.get("CEN_LAT", None)
self._cen_lon = proj_params.get("CEN_LON", None) self._cen_lon = up_proj_params.get("CEN_LON", None)
self.truelat1 = proj_params.get("TRUELAT1", None) self.truelat1 = up_proj_params.get("TRUELAT1", None)
self.truelat2 = (proj_params.get("TRUELAT2", None) self.truelat2 = (up_proj_params.get("TRUELAT2", None)
if not _ismissing(proj_params.get("TRUELAT2", None)) if not _ismissing(up_proj_params.get("TRUELAT2",
None))
else None) else None)
self.moad_cen_lat = proj_params.get("MOAD_CEN_LAT", None) self.moad_cen_lat = up_proj_params.get("MOAD_CEN_LAT", None)
self.stand_lon = proj_params.get("STAND_LON", None) self.stand_lon = up_proj_params.get("STAND_LON", None)
self.pole_lat = proj_params.get("POLE_LAT", None) self.pole_lat = up_proj_params.get("POLE_LAT", None)
self.pole_lon = proj_params.get("POLE_LON", None) self.pole_lon = up_proj_params.get("POLE_LON", None)
self.dx = up_proj_params.get("DX", None)
self.dy = up_proj_params.get("DY", None)
# Just in case... # Just in case...
if self.moad_cen_lat is None: if self.moad_cen_lat is None:
@ -194,6 +224,7 @@ class WrfProj(object):
if self.stand_lon is None: if self.stand_lon is None:
self.stand_lon = self._cen_lon self.stand_lon = self._cen_lon
def _basemap(self, resolution='l'): def _basemap(self, resolution='l'):
return None return None
@ -396,8 +427,6 @@ class LambertConformal(WrfProj):
have the same names as found in WRF output NetCDF global have the same names as found in WRF output NetCDF global
attributes: attributes:
- 'CEN_LAT': Center latitude.
- 'CEN_LON': Center longitude.
- 'TRUELAT1': True latitude 1. - 'TRUELAT1': True latitude 1.
- 'TRUELAT2': True latitude 2. - 'TRUELAT2': True latitude 2.
- 'MOAD_CEN_LAT': Mother of all domains center latitude. - 'MOAD_CEN_LAT': Mother of all domains center latitude.
@ -542,8 +571,6 @@ class Mercator(WrfProj):
have the same names as found in WRF output NetCDF global have the same names as found in WRF output NetCDF global
attributes: attributes:
- 'CEN_LAT': Center latitude.
- 'CEN_LON': Center longitude.
- 'TRUELAT1': True latitude 1. - 'TRUELAT1': True latitude 1.
- 'TRUELAT2': True latitude 2. - 'TRUELAT2': True latitude 2.
- 'MOAD_CEN_LAT': Mother of all domains center latitude. - 'MOAD_CEN_LAT': Mother of all domains center latitude.
@ -686,8 +713,6 @@ class PolarStereographic(WrfProj):
have the same names as found in WRF output NetCDF global have the same names as found in WRF output NetCDF global
attributes: attributes:
- 'CEN_LAT': Center latitude.
- 'CEN_LON': Center longitude.
- 'TRUELAT1': True latitude 1. - 'TRUELAT1': True latitude 1.
- 'TRUELAT2': True latitude 2. - 'TRUELAT2': True latitude 2.
- 'MOAD_CEN_LAT': Mother of all domains center latitude. - 'MOAD_CEN_LAT': Mother of all domains center latitude.
@ -828,8 +853,6 @@ class LatLon(WrfProj):
have the same names as found in WRF output NetCDF global have the same names as found in WRF output NetCDF global
attributes: attributes:
- 'CEN_LAT': Center latitude.
- 'CEN_LON': Center longitude.
- 'TRUELAT1': True latitude 1. - 'TRUELAT1': True latitude 1.
- 'TRUELAT2': True latitude 2. - 'TRUELAT2': True latitude 2.
- 'MOAD_CEN_LAT': Mother of all domains center latitude. - 'MOAD_CEN_LAT': Mother of all domains center latitude.
@ -971,8 +994,6 @@ class RotatedLatLon(WrfProj):
have the same names as found in WRF output NetCDF global have the same names as found in WRF output NetCDF global
attributes: attributes:
- 'CEN_LAT': Center latitude.
- 'CEN_LON': Center longitude.
- 'TRUELAT1': True latitude 1. - 'TRUELAT1': True latitude 1.
- 'TRUELAT2': True latitude 2. - 'TRUELAT2': True latitude 2.
- 'MOAD_CEN_LAT': Mother of all domains center latitude. - 'MOAD_CEN_LAT': Mother of all domains center latitude.
@ -1120,13 +1141,13 @@ def getproj(bottom_left=None, top_right=None,
Args: Args:
bottom_left (indexable sequence, optional): The lower left corner as bottom_left (:class:`wrf.CoordPair`, optional): The lower left
a (latitude, longitude) pair. Must also specify *top_right* if corner. Must also specify *top_right* if used.
used. Default is None. Default is None.
top_right (indexable sequence): The upper right corner as a top_right (:class:`wrf.CoordPair`, optional): The upper right
(latitude, longitude) pair. Must also specify *bottom_left* corner. Must also specify *bottom_left* if used.
if used. Default is None. Default is None.
lats (:class:`numpy.ndarray`, optional): An array of at least lats (:class:`numpy.ndarray`, optional): An array of at least
two dimensions containing all of the latitude values. Must two dimensions containing all of the latitude values. Must
@ -1141,8 +1162,6 @@ def getproj(bottom_left=None, top_right=None,
attributes: attributes:
- 'MAP_PROJ': The map projection type as an integer. - 'MAP_PROJ': The map projection type as an integer.
- 'CEN_LAT': Center latitude.
- 'CEN_LON': Center longitude.
- 'TRUELAT1': True latitude 1. - 'TRUELAT1': True latitude 1.
- 'TRUELAT2': True latitude 2. - 'TRUELAT2': True latitude 2.
- 'MOAD_CEN_LAT': Mother of all domains center latitude. - 'MOAD_CEN_LAT': Mother of all domains center latitude.
@ -1157,22 +1176,24 @@ def getproj(bottom_left=None, top_right=None,
""" """
proj_type = proj_params.get("MAP_PROJ", 0) up_proj_params = dict_keys_to_upper(proj_params)
proj_type = up_proj_params.get("MAP_PROJ", 0)
if proj_type == ProjectionTypes.LAMBERT_CONFORMAL: if proj_type == ProjectionTypes.LAMBERT_CONFORMAL:
return LambertConformal(bottom_left, top_right, return LambertConformal(bottom_left, top_right,
lats, lons, **proj_params) lats, lons, **proj_params)
elif proj_type == ProjectionTypes.POLAR_STEREOGRAPHIC: elif proj_type == ProjectionTypes.POLAR_STEREOGRAPHIC:
return PolarStereographic(bottom_left, top_right, return PolarStereographic(bottom_left, top_right,
lats, lons, **proj_params) lats, lons, **proj_params)
elif proj_type == ProjectionTypes.MERCATOR: elif proj_type == ProjectionTypes.MERCATOR:
return Mercator(bottom_left, top_right, return Mercator(bottom_left, top_right,
lats, lons, **proj_params) lats, lons, **proj_params)
elif (proj_type == ProjectionTypes.ZERO or elif (proj_type == ProjectionTypes.ZERO or
proj_type == ProjectionTypes.LAT_LON): proj_type == ProjectionTypes.LAT_LON):
if (proj_params.get("POLE_LAT", None) == 90. if (up_proj_params.get("POLE_LAT", None) == 90.
and proj_params.get("POLE_LON", None) == 0.): and up_proj_params.get("POLE_LON", None) == 0.):
return LatLon(bottom_left, top_right, return LatLon(bottom_left, top_right,
lats, lons, **proj_params) lats, lons, **proj_params)
else: else:
return RotatedLatLon(bottom_left, top_right, return RotatedLatLon(bottom_left, top_right,
lats, lons, **proj_params) lats, lons, **proj_params)

18
src/wrf/projutils.py

@ -0,0 +1,18 @@
from __future__ import (absolute_import, division, print_function,
unicode_literals)
from .py3compat import viewitems
def dict_keys_to_upper(d):
"""Return a dictionary with the keys changed to uppercase.
Args:
d (:obj:`dict`): A dictionary.
Returns:
:obj:`dict`: A dictionary with uppercase keys.
"""
return {key.upper() : val for key, val in viewitems(d)}

2
src/wrf/py3compat.py

@ -166,3 +166,5 @@ def ucode(*args, **kwargs):
return str(*args, **kwargs) return str(*args, **kwargs)
return unicode(*args, **kwargs) return unicode(*args, **kwargs)

10
src/wrf/util.py

@ -2596,7 +2596,8 @@ def get_proj_params(wrfin, timeidx=0, varname=None):
"CEN_LAT", "CEN_LON", "CEN_LAT", "CEN_LON",
"TRUELAT1", "TRUELAT2", "TRUELAT1", "TRUELAT2",
"MOAD_CEN_LAT", "STAND_LON", "MOAD_CEN_LAT", "STAND_LON",
"POLE_LAT", "POLE_LON")) "POLE_LAT", "POLE_LON",
"DX", "DY"))
multitime = is_multi_time_req(timeidx) multitime = is_multi_time_req(timeidx)
if not multitime: if not multitime:
time_idx_or_slice = timeidx time_idx_or_slice = timeidx
@ -2947,6 +2948,13 @@ def get_id(obj):

Loading…
Cancel
Save