From 9be5a86a1b2c0a92d569e1306b5d57e43535c0d7 Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Thu, 27 Oct 2016 17:03:13 -0600 Subject: [PATCH] Now uses CoordPair for vertical cross sections so that lat/lon points can be used directly. --- src/wrf/interp.py | 225 ++++++++++++++++++++++++++++++++------ src/wrf/interputils.py | 116 +++++++++++++++++++- src/wrf/latlon.py | 33 ++++-- src/wrf/latlonutils.py | 25 ++--- src/wrf/metadecorators.py | 86 +++++++++++++-- src/wrf/projection.py | 121 +++++++++++--------- src/wrf/projutils.py | 18 +++ src/wrf/py3compat.py | 2 + src/wrf/util.py | 10 +- 9 files changed, 511 insertions(+), 125 deletions(-) create mode 100644 src/wrf/projutils.py diff --git a/src/wrf/interp.py b/src/wrf/interp.py index 4a92b76..e69787b 100755 --- a/src/wrf/interp.py +++ b/src/wrf/interp.py @@ -10,7 +10,7 @@ from .extension import (_interpz3d, _vertcross, _interpline, _smooth2d, from .metadecorators import set_interp_metadata from .util import extract_vars, is_staggered, get_id, npvalues 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 .terrain import get_terrain from .geoht import get_height @@ -90,7 +90,8 @@ def interplevel(field3d, vert, desiredlev, missing=Constants.DEFAULT_FILL, @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, start_point=None, end_point=None, 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. This horizontal line is defined by either including the *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 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* 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 is provided. Default is None. + vertical levels in the output array. Must be in the same units + 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. Default is :data:`wrf.Constants.DEFAULT_FILL`. - - pivot_point (:obj:`tuple` or :obj:`list`, optional): A - :obj:`tuple` or :obj:`list` with two entries, - in the form of [x, y] (or [west_east, south_north]), which - indicates the x,y location through which the plane will pass. - Must also specify `angle`. + + 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. 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 a plane will be plotted through a given point on the model domain. 0.0 represents a S-N cross section. 90.0 is a W-E cross section. - start_point (:obj:`tuple` or :obj:`list`, optional): A - :obj:`tuple` or :obj:`list` with two entries, in the form of - [x, y] (or [west_east, south_north]), which indicates the start - x,y location through which the plane will pass. + start_point (:class:`wrf.CoordPair`, optional): A coordinate pair + which indicates the start 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. - end_point (:obj:`tuple` or :obj:`list`, optional): A - :obj:`tuple` or :obj:`list` with two entries, in the form of - [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 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. """ + 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 # type, we'll handle that iteration here. 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"] z_var2d = cache["z_var2d"] except (KeyError, TypeError): - xy, var2dz, z_var2d = get_xy_z_params(npvalues(vert), pivot_point, - angle, start_point, end_point, - levels) + # Convert Lat/Lon points to grid points + 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, var2dz, z_var2d = get_xy_z_params(npvalues(vert), pivot_point_xy, + angle, start_point_xy, + end_point_xy, levels) if not multi: 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") def interpline(field2d, pivot_point=None, + wrfin=None, timeidx=0, stagger=None, projection=None, angle=None, start_point=None, end_point=None, latlon=False, cache=None, meta=True): """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: field2d (:class:`xarray.DataArray` or :class:`numpy.ndarray`): A two-dimensional field. - pivot_point (:obj:`tuple` or :obj:`list`, optional): A - :obj:`tuple` or :obj:`list` with two entries, - in the form of [x, y] (or [west_east, south_north]), which - indicates the x,y location through which the plane will pass. - Must also specify `angle`. + 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`, 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 a plane will be plotted through a given point on the model domain. 0.0 represents a S-N cross section. 90.0 is a W-E cross section. - start_point (:obj:`tuple` or :obj:`list`, optional): A - :obj:`tuple` or :obj:`list` with two entries, in the form of - [x, y] (or [west_east, south_north]), which indicates the start - x,y location through which the plane will pass. + start_point (:class:`wrf.CoordPair`, optional): A coordinate pair + which indicates the start 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. - end_point (:obj:`tuple` or :obj:`list`, optional): A - :obj:`tuple` or :obj:`list` with two entries, in the form of - [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 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: xy = cache["xy"] 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) diff --git a/src/wrf/interputils.py b/src/wrf/interputils.py index 2f866bf..0faa989 100644 --- a/src/wrf/interputils.py +++ b/src/wrf/interputils.py @@ -2,11 +2,15 @@ from __future__ import (absolute_import, division, print_function, unicode_literals) from math import floor, ceil +from collections import Iterable import numpy as np from .extension import _interp2dxy 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): @@ -256,11 +260,11 @@ def get_xy_z_params(z, pivot_point=None, angle=None, nlevels = int(z_max/dz) z_var2d = np.zeros((nlevels), dtype=z.dtype) z_var2d[0] = z_min + + for i in py3range(1,nlevels): + z_var2d[i] = z_var2d[0] + i*dz else: - z_var2d = np.asarray(levels) - - for i in py3range(1,nlevels): - z_var2d[i] = z_var2d[0] + i*dz + z_var2d = np.asarray(levels, z.dtype) 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) 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])] diff --git a/src/wrf/latlon.py b/src/wrf/latlon.py index c333ac9..3cbb527 100755 --- a/src/wrf/latlon.py +++ b/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. Required for *map_proj* = 1, 2, 3 (defaults to 0 otherwise). - latinc (:obj:`float`): The lower left corner latitude at (0,0). - Required for *map_proj* = 6. + latinc (:obj:`float`): Required for *map_proj* = 6. Defined as: + + .. code-block:: python + + latinc = (dy*360.0)/2.0/Constants.PI/Constants.WRF_EARTH_RADIUS - loninc (:obj:`float`) - The lower left corner longitude at (0,0). - Required for *map_proj* = 6. + 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: :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. Required for *map_proj* = 1, 2, 3 (defaults to 0 otherwise). - latinc (:obj:`float`): The lower left corner latitude at (0,0). - Required for *map_proj* = 6. + latinc (:obj:`float`): Required for *map_proj* = 6. Defined as: + + .. code-block:: python + + latinc = (dy*360.0)/2.0/Constants.PI/Constants.WRF_EARTH_RADIUS - loninc (:obj:`float`) - The lower left corner longitude at (0,0). - Required for *map_proj* = 6. + 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: :class:`xarray.DataArray` or :class:`numpy.ndarray`: The @@ -472,6 +484,11 @@ def xy_to_ll_proj(x, y, meta=True, squeeze=True, map_proj=None, truelat1=None, "latinc", "loninc")} return _xy_to_ll(x, y, None, 0, None, "cat", squeeze, None, None, **projparams) + + + + + \ No newline at end of file diff --git a/src/wrf/latlonutils.py b/src/wrf/latlonutils.py index 527f445..00aee6b 100644 --- a/src/wrf/latlonutils.py +++ b/src/wrf/latlonutils.py @@ -11,6 +11,8 @@ from .util import (extract_vars, extract_global_attrs, either, is_moving_domain, is_multi_time_req, iter_left_indexes, is_mapping, is_multi_file) from .py3compat import viewkeys, viewitems +from .projutils import dict_keys_to_upper + def _lat_varname(wrfin, stagger): """Return the latitude variable name for the specified stagger type. @@ -192,21 +194,7 @@ def _get_proj_params(wrfin, timeidx, stagger, method, squeeze, cache, _key): return (map_proj, truelat1, truelat2, stdlon, ref_lat, ref_lon, 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 def _kwarg_proj_params(**projparams): @@ -224,7 +212,7 @@ def _kwarg_proj_params(**projparams): :obj:`tuple`: The map projection parameters. """ - projparams = _dict_keys_to_upper(projparams) + projparams = dict_keys_to_upper(projparams) map_proj = projparams.get("MAP_PROJ") 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) return result + + + + + diff --git a/src/wrf/metadecorators.py b/src/wrf/metadecorators.py index 8314349..0eb8e3a 100644 --- a/src/wrf/metadecorators.py +++ b/src/wrf/metadecorators.py @@ -12,7 +12,7 @@ from .util import (extract_vars, either, from_args, arg_location, from_var, iter_left_indexes) from .coordpair import CoordPair 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 if xarray_enabled(): @@ -830,7 +830,9 @@ def _set_cross_meta(wrapped, instance, args, kwargs): :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", "start_point", "end_point", "cache"), @@ -838,16 +840,48 @@ def _set_cross_meta(wrapped, instance, args, kwargs): field3d = argvars["field3d"] z = argvars["vert"] + levels = argvars["levels"] inc_latlon = argvars["latlon"] missingval = argvars["missing"] + wrfin = argvars["wrfin"] + timeidx = argvars["timeidx"] + stagger = argvars["stagger"] + projection = argvars["projection"] pivot_point = argvars["pivot_point"] angle = argvars["angle"] start_point = argvars["start_point"] end_point = argvars["end_point"] cache = argvars["cache"] - xy, var2dz, z_var2d = get_xy_z_params(npvalues(z), 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, 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 if cache is not None: @@ -877,8 +911,7 @@ def _set_cross_meta(wrapped, instance, args, kwargs): ed_x = xy[-1,0] ed_y = xy[-1,1] - cross_str = "({0}, {1}) to ({2}, {3})".format(st_x, st_y, - ed_x, ed_y) + cross_str = "({0}, {1}) to ({2}, {3})".format(st_x, st_y, ed_x, ed_y) if angle is not None: cross_str += " ; center={0} ; angle={1}".format(pivot_point, angle) @@ -1009,12 +1042,18 @@ def _set_line_meta(wrapped, instance, args, kwargs): :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", "cache"), *args, **kwargs) field2d = argvars["field2d"] + wrfin = argvars["wrfin"] + timeidx = argvars["timeidx"] + stagger = argvars["stagger"] + projection = argvars["projection"] pivot_point = argvars["pivot_point"] angle = argvars["angle"] start_point = argvars["start_point"] @@ -1024,8 +1063,36 @@ def _set_line_meta(wrapped, instance, args, kwargs): if cache is None: cache = {} + + 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) - xy = get_xy(field2d, pivot_point, angle, start_point, end_point) + + 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 new_cache = dict(cache) @@ -1050,8 +1117,7 @@ def _set_line_meta(wrapped, instance, args, kwargs): ed_x = xy[-1,0] ed_y = xy[-1,1] - cross_str = "({0}, {1}) to ({2}, {3})".format(st_x, st_y, - ed_x, ed_y) + cross_str = "({0}, {1}) to ({2}, {3})".format(st_x, st_y, ed_x, ed_y) if angle is not None: cross_str += " ; center={0} ; angle={1}".format(pivot_point, angle) diff --git a/src/wrf/projection.py b/src/wrf/projection.py index ab286a9..d706181 100644 --- a/src/wrf/projection.py +++ b/src/wrf/projection.py @@ -3,8 +3,10 @@ from __future__ import (absolute_import, division, print_function, import numpy as np import math +from .coordpair import CoordPair from .config import basemap_enabled, cartopy_enabled, pyngl_enabled from .constants import Constants, ProjectionTypes +from .projutils import dict_keys_to_upper if cartopy_enabled(): from cartopy import crs @@ -78,6 +80,7 @@ if cartopy_enabled(): self._threshold = np.diff(self.x_limits)[0] / 720 + def _ismissing(val): """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. + class WrfProj(object): """A base class for storing map projection information from WRF data. @@ -116,6 +120,25 @@ class WrfProj(object): bottom_left (indexable sequence): A pair of (ll_lat, ll_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, @@ -124,13 +147,13 @@ class WrfProj(object): Args: - bottom_left (indexable sequence, optional): The lower left corner - as a (latitude, longitude) pair. Must also specify *top_right* - if used. Default is None. + bottom_left (:class:`wrf.CoordPair`, optional): The lower left + corner. Must also specify *top_right* if used. + Default is None. - top_right (indexable sequence): The upper right corner as a - (latitude, longitude) pair. Must also specify *bottom_left* - if used. Default is None. + top_right (:class:`wrf.CoordPair`, optional): The upper right + corner. Must also specify *bottom_left* if used. + Default is None. lats (:class:`numpy.ndarray`, optional): An array of at least two dimensions containing all of the latitude values. Must @@ -142,11 +165,9 @@ class WrfProj(object): **proj_params: Map projection optional keyword arguments, that have the same names as found in WRF output NetCDF global - attributes: + attributes (case insensitive): - 'MAP_PROJ': The map projection type as an integer. - - 'CEN_LAT': Center latitude. - - 'CEN_LON': Center longitude. - 'TRUELAT1': True latitude 1. - 'TRUELAT2': True latitude 2. - '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: - self.ll_lat = bottom_left[0] - self.ll_lon = bottom_left[1] - self.ur_lat = top_right[0] - self.ur_lon = top_right[1] + self.ll_lat = bottom_left.lat + self.ll_lon = bottom_left.lon + self.ur_lat = top_right.lat + self.ur_lon = top_right.lon self.bottom_left = bottom_left self.top_right = top_right elif lats is not None and lons is not None: @@ -168,24 +189,33 @@ class WrfProj(object): self.ur_lat = lats[-1,-1] self.ll_lon = lons[0,0] self.ur_lon = lons[-1,-1] - self.bottom_left = [self.ll_lat, self.ll_lon] - self.top_right = [self.ur_lat, self.ur_lon] + self.bottom_left = CoordPair(lat=self.ll_lat, lon=self.ll_lon) + self.top_right = CoordPair(self.ur_lat, self.ur_lon) else: 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 # center of the projection - self._cen_lat = proj_params.get("CEN_LAT", None) - self._cen_lon = proj_params.get("CEN_LON", None) + self._cen_lat = up_proj_params.get("CEN_LAT", None) + self._cen_lon = up_proj_params.get("CEN_LON", None) - self.truelat1 = proj_params.get("TRUELAT1", None) - self.truelat2 = (proj_params.get("TRUELAT2", None) - if not _ismissing(proj_params.get("TRUELAT2", None)) + self.truelat1 = up_proj_params.get("TRUELAT1", None) + self.truelat2 = (up_proj_params.get("TRUELAT2", None) + if not _ismissing(up_proj_params.get("TRUELAT2", + None)) else None) - self.moad_cen_lat = proj_params.get("MOAD_CEN_LAT", None) - self.stand_lon = proj_params.get("STAND_LON", None) - self.pole_lat = proj_params.get("POLE_LAT", None) - self.pole_lon = proj_params.get("POLE_LON", None) + self.moad_cen_lat = up_proj_params.get("MOAD_CEN_LAT", None) + self.stand_lon = up_proj_params.get("STAND_LON", None) + self.pole_lat = up_proj_params.get("POLE_LAT", 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... if self.moad_cen_lat is None: @@ -194,6 +224,7 @@ class WrfProj(object): if self.stand_lon is None: self.stand_lon = self._cen_lon + def _basemap(self, resolution='l'): return None @@ -396,8 +427,6 @@ class LambertConformal(WrfProj): have the same names as found in WRF output NetCDF global attributes: - - 'CEN_LAT': Center latitude. - - 'CEN_LON': Center longitude. - 'TRUELAT1': True latitude 1. - 'TRUELAT2': True latitude 2. - '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 attributes: - - 'CEN_LAT': Center latitude. - - 'CEN_LON': Center longitude. - 'TRUELAT1': True latitude 1. - 'TRUELAT2': True latitude 2. - '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 attributes: - - 'CEN_LAT': Center latitude. - - 'CEN_LON': Center longitude. - 'TRUELAT1': True latitude 1. - 'TRUELAT2': True latitude 2. - '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 attributes: - - 'CEN_LAT': Center latitude. - - 'CEN_LON': Center longitude. - 'TRUELAT1': True latitude 1. - 'TRUELAT2': True latitude 2. - '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 attributes: - - 'CEN_LAT': Center latitude. - - 'CEN_LON': Center longitude. - 'TRUELAT1': True latitude 1. - 'TRUELAT2': True latitude 2. - 'MOAD_CEN_LAT': Mother of all domains center latitude. @@ -1120,13 +1141,13 @@ def getproj(bottom_left=None, top_right=None, Args: - bottom_left (indexable sequence, optional): The lower left corner as - a (latitude, longitude) pair. Must also specify *top_right* if - used. Default is None. - - top_right (indexable sequence): The upper right corner as a - (latitude, longitude) pair. Must also specify *bottom_left* - if used. Default is None. + bottom_left (:class:`wrf.CoordPair`, optional): The lower left + corner. Must also specify *top_right* if used. + Default is None. + + top_right (:class:`wrf.CoordPair`, optional): The upper right + corner. Must also specify *bottom_left* if used. + Default is None. lats (:class:`numpy.ndarray`, optional): An array of at least two dimensions containing all of the latitude values. Must @@ -1141,8 +1162,6 @@ def getproj(bottom_left=None, top_right=None, attributes: - 'MAP_PROJ': The map projection type as an integer. - - 'CEN_LAT': Center latitude. - - 'CEN_LON': Center longitude. - 'TRUELAT1': True latitude 1. - 'TRUELAT2': True latitude 2. - '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: return LambertConformal(bottom_left, top_right, - lats, lons, **proj_params) + lats, lons, **proj_params) elif proj_type == ProjectionTypes.POLAR_STEREOGRAPHIC: return PolarStereographic(bottom_left, top_right, - lats, lons, **proj_params) + lats, lons, **proj_params) elif proj_type == ProjectionTypes.MERCATOR: return Mercator(bottom_left, top_right, - lats, lons, **proj_params) + lats, lons, **proj_params) elif (proj_type == ProjectionTypes.ZERO or proj_type == ProjectionTypes.LAT_LON): - if (proj_params.get("POLE_LAT", None) == 90. - and proj_params.get("POLE_LON", None) == 0.): + if (up_proj_params.get("POLE_LAT", None) == 90. + and up_proj_params.get("POLE_LON", None) == 0.): return LatLon(bottom_left, top_right, - lats, lons, **proj_params) + lats, lons, **proj_params) else: return RotatedLatLon(bottom_left, top_right, lats, lons, **proj_params) diff --git a/src/wrf/projutils.py b/src/wrf/projutils.py new file mode 100644 index 0000000..c7ff210 --- /dev/null +++ b/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)} diff --git a/src/wrf/py3compat.py b/src/wrf/py3compat.py index 84b8ef8..c74ce79 100644 --- a/src/wrf/py3compat.py +++ b/src/wrf/py3compat.py @@ -166,3 +166,5 @@ def ucode(*args, **kwargs): return str(*args, **kwargs) return unicode(*args, **kwargs) + + diff --git a/src/wrf/util.py b/src/wrf/util.py index ce4fa5e..b4a1333 100644 --- a/src/wrf/util.py +++ b/src/wrf/util.py @@ -2596,7 +2596,8 @@ def get_proj_params(wrfin, timeidx=0, varname=None): "CEN_LAT", "CEN_LON", "TRUELAT1", "TRUELAT2", "MOAD_CEN_LAT", "STAND_LON", - "POLE_LAT", "POLE_LON")) + "POLE_LAT", "POLE_LON", + "DX", "DY")) multitime = is_multi_time_req(timeidx) if not multitime: time_idx_or_slice = timeidx @@ -2930,9 +2931,16 @@ def get_id(obj): # For each key in the mapping, recursively call get_id until # until a non-mapping is found return {key : get_id(val) for key,val in viewitems(obj)} + + + + + + +