From 812c546d9520a83ccd8b71e68f0bbb53a701f951 Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Wed, 24 Jan 2018 14:13:29 -0700 Subject: [PATCH] Changed vertcross and interpline to use a CoordPair for the lower left corner point. Also added some helper functions to extract the lower left points from 2D lat/lon arrays and extract lat/lons from sequences of CoordPair objects. Updated documentation and unit tests. --- doc/source/user_api/index.rst | 2 ++ src/wrf/api.py | 4 +-- src/wrf/coordpair.py | 4 +-- src/wrf/interp.py | 44 ++++++++++--------------- src/wrf/interputils.py | 35 ++++++++------------ src/wrf/latlonutils.py | 6 ++-- src/wrf/metadecorators.py | 24 ++++++-------- src/wrf/util.py | 62 +++++++++++++++++++++++++++++++++++ test/utests.py | 12 +++---- 9 files changed, 117 insertions(+), 76 deletions(-) diff --git a/doc/source/user_api/index.rst b/doc/source/user_api/index.rst index 0832ca5..5d5d77d 100644 --- a/doc/source/user_api/index.rst +++ b/doc/source/user_api/index.rst @@ -292,6 +292,8 @@ them helpful for other purposes. wrf.getproj wrf.cache_item wrf.get_cached_item + wrf.ll_points + wrf.pairs_to_latlon ------------------------ diff --git a/src/wrf/api.py b/src/wrf/api.py index 660c9d6..46f7d46 100644 --- a/src/wrf/api.py +++ b/src/wrf/api.py @@ -45,7 +45,7 @@ from .util import (to_np, extract_global_attrs, is_standard_wrf_var, has_time_coord, is_multi_file, is_multi_time_req, get_coord_pairs, is_time_coord_var, geo_bounds, get_cartopy, get_basemap, get_pyngl, cartopy_xlim, - cartopy_ylim, latlon_coords) + cartopy_ylim, latlon_coords, ll_points, pairs_to_latlon) from .geobnds import GeoBounds, NullGeoBounds from .projection import (WrfProj, NullProjection, LambertConformal, Mercator, PolarStereographic, LatLon, RotatedLatLon, @@ -102,7 +102,7 @@ __all__ += ["to_np", "extract_global_attrs", "is_standard_wrf_var", "has_time_coord", "is_multi_file", "is_multi_time_req", "get_coord_pairs", "is_time_coord_var", "geo_bounds", "get_cartopy", "get_basemap", "get_pyngl", "cartopy_xlim", - "cartopy_ylim", "latlon_coords"] + "cartopy_ylim", "latlon_coords", "ll_points", "pairs_to_latlon"] __all__ += ["GeoBounds", "NullGeoBounds"] __all__ += ["WrfProj", "NullProjection", "LambertConformal", "Mercator", "PolarStereographic", "LatLon", "RotatedLatLon", "getproj"] diff --git a/src/wrf/coordpair.py b/src/wrf/coordpair.py index ad42332..8b3121a 100644 --- a/src/wrf/coordpair.py +++ b/src/wrf/coordpair.py @@ -3,7 +3,7 @@ from __future__ import (absolute_import, division, print_function, from .py3compat import py2round - + def _binary_operator(operator): """Function wrapper for binary operators. @@ -262,5 +262,5 @@ for operator in ("__neg__", "__pos__", "__abs__", "__invert__"): for operator in ("__lt__", "__le__", "__eq__", "__ne__", "__gt__", "__ge__"): setattr(CoordPair, operator, _cmp_operator(operator)) - + \ No newline at end of file diff --git a/src/wrf/interp.py b/src/wrf/interp.py index 9e5ead6..1892947 100755 --- a/src/wrf/interp.py +++ b/src/wrf/interp.py @@ -92,7 +92,7 @@ def interplevel(field3d, vert, desiredlev, missing=default_fill(np.float64), @set_interp_metadata("cross") def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64), wrfin=None, timeidx=0, stagger=None, projection=None, - ll_lat=None, ll_lon=None, + ll_point=None, pivot_point=None, angle=None, start_point=None, end_point=None, latlon=False, cache=None, meta=True): @@ -166,15 +166,11 @@ def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64), coordinates, and must be specified if *wrfin* is None. Default is None. - ll_lat (:obj:`float`, sequence of :obj:`float`, optional): The lower - left latitude(s) for your domain, and must be specified if - *wrfin* is None. If the domain is a moving nest, this should be a - sequence of lower left latitudes. Default is None. - - ll_lon (:obj:`float`, sequence of :obj:`float`, optional): The lower - left longitude(s) for your domain, and must be specified if - *wrfin* is None. If the domain is a moving nest, this should be - a sequence of lower left longitudes. Default is None. + ll_point (:class:`wrf.CoordPair`, sequence of :class:`wrf.CoordPair`, \ + optional): The lower left latitude, longitude point for your domain, + and must be specified + if *wrfin* is None. If the domain is a moving nest, this should be + a sequence of :class:`wrf.CoordPair`. Default is None. pivot_point (:class:`wrf.CoordPair`, optional): A coordinate pair for the pivot point, which indicates the location through which @@ -262,7 +258,7 @@ def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64), 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) + stagger, projection, ll_point) pivot_point_xy = (xy_coords.x, xy_coords.y) else: pivot_point_xy = (pivot_point.x, pivot_point.y) @@ -270,14 +266,14 @@ def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64), 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) + stagger, projection, ll_point) 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) + stagger, projection, ll_point) end_point_xy = (xy_coords.x, xy_coords.y) else: end_point_xy = (end_point.x, end_point.y) @@ -302,7 +298,7 @@ def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64), @set_interp_metadata("line") def interpline(field2d, pivot_point=None, wrfin=None, timeidx=0, stagger=None, projection=None, - ll_lat=None, ll_lon=None, + ll_point=None, angle=None, start_point=None, end_point=None, latlon=False, cache=None, meta=True): @@ -352,15 +348,11 @@ def interpline(field2d, pivot_point=None, not be used when working with x,y coordinates. Default is None. - ll_lat (:obj:`float`, sequence of :obj:`float`, optional): The lower - left latitude(s) for your domain, and must be specified if - *wrfin* is None. If the domain is a moving nest, this should be a - sequence of lower left latitudes. Default is None. - - ll_lon (:obj:`float`, sequence of :obj:`float`, optional): The lower - left longitude(s) for your domain, and must be specified if - *wrfin* is None. If the domain is a moving nest, this should be - a sequence of lower left longitudes. Default is None. + ll_point (:class:`wrf.CoordPair`, sequence of :class:`wrf.CoordPair`, \ + optional): The lower left latitude, longitude point for your domain, + and must be specified + if *wrfin* is None. If the domain is a moving nest, this should be + a sequence of :class:`wrf.CoordPair`. Default is None. pivot_point (:class:`wrf.CoordPair`, optional): A coordinate pair for the pivot point, which indicates the location through which @@ -442,7 +434,7 @@ def interpline(field2d, pivot_point=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, ll_lat, ll_lon) + stagger, projection, ll_point) pivot_point_xy = (xy_coords.x, xy_coords.y) else: pivot_point_xy = (pivot_point.x, pivot_point.y) @@ -450,14 +442,14 @@ def interpline(field2d, pivot_point=None, 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, ll_lat, ll_lon) + stagger, projection, ll_point) 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, ll_lat, ll_lon) + stagger, projection, ll_point) end_point_xy = (xy_coords.x, xy_coords.y) else: end_point_xy = (end_point.x, end_point.y) diff --git a/src/wrf/interputils.py b/src/wrf/interputils.py index 5289f08..cd4ddc1 100644 --- a/src/wrf/interputils.py +++ b/src/wrf/interputils.py @@ -2,7 +2,6 @@ from __future__ import (absolute_import, division, print_function, unicode_literals) from math import floor, ceil -from collections import Iterable import numpy as np @@ -11,6 +10,7 @@ from .py3compat import py3range from .coordpair import CoordPair from .constants import Constants, ProjectionTypes from .latlonutils import _ll_to_xy +from .util import pairs_to_latlon def to_positive_idxs(shape, coord): @@ -330,7 +330,7 @@ def get_xy(var, pivot_point=None, angle=None, return xy def to_xy_coords(pairs, wrfin=None, timeidx=0, stagger=None, projection=None, - ll_lat=None, ll_lon=None): + ll_point=None): """Return the coordinate pairs in grid space. This function converts latitude,longitude coordinate pairs to @@ -373,15 +373,11 @@ def to_xy_coords(pairs, wrfin=None, timeidx=0, stagger=None, projection=None, coordinates, and must be specified if *wrfin* is None. Default is None. - ll_lat (:obj:`float`, sequence of :obj:`float`, optional): The lower - left latitude(s) for your domain, and must be specified if - *wrfin* is None. If the domain is a moving nest, this should be a - sequence of lower left latitudes. Default is None. - - ll_lon (:obj:`float`, sequence of :obj:`float`, optional): The lower - left longitude(s) for your domain, and must be specified if - *wrfin* is None. If the domain is a moving nest, this should be - a sequence of lower left longitudes. Default is None. + ll_point (:class:`wrf.CoordPair`, sequence of :class:`wrf.CoordPair`, \ + optional): The lower left latitude, longitude point for your domain, + and must be specified + if *wrfin* is None. If the domain is a moving nest, this should be + a sequence of :class:`wrf.CoordPair`. Default is None. Returns: @@ -390,19 +386,13 @@ def to_xy_coords(pairs, wrfin=None, timeidx=0, stagger=None, projection=None, """ - if (wrfin is None and - (projection is None or ll_lat is None or ll_lon is None)): + if (wrfin is None and (projection is None or ll_point is None)): raise ValueError ("'wrfin' parameter or " - "'projection', 'll_lat', and 'll_lon' parameters " + "'projection' and 'll_point' parameters " "are 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 - + lat, lon = pairs_to_latlon(pairs) + 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) @@ -422,7 +412,8 @@ def to_xy_coords(pairs, wrfin=None, timeidx=0, stagger=None, projection=None, pole_lon = 0.0 latinc = 0.0 loninc = 0.0 - + + ll_lat, ll_lon = pairs_to_latlon(ll_point) xy_vals = _ll_to_xy(lat, lon, meta=False, squeeze=True, as_int=True, map_proj=projection.map_proj, diff --git a/src/wrf/latlonutils.py b/src/wrf/latlonutils.py index b050ae1..ccd5221 100644 --- a/src/wrf/latlonutils.py +++ b/src/wrf/latlonutils.py @@ -242,9 +242,9 @@ def _kwarg_proj_params(**projparams): if var is None: raise ValueError("'{}' argument required".format(name)) - # ref_lat and ref_lon are expectd to be lists - ref_lat = np.asarray([ref_lat]) - ref_lon = np.asarray([ref_lon]) + # ref_lat and ref_lon are expected to be lists + ref_lat = np.ravel(np.asarray([ref_lat])) + ref_lon = np.ravel(np.asarray([ref_lon])) # Fortran wants 1-based indexing known_x = known_x + 1 diff --git a/src/wrf/metadecorators.py b/src/wrf/metadecorators.py index 2abcb63..27425c1 100644 --- a/src/wrf/metadecorators.py +++ b/src/wrf/metadecorators.py @@ -887,8 +887,7 @@ def _set_cross_meta(wrapped, instance, args, kwargs): argvars = from_args(wrapped, ("field3d", "vert", "levels", "latlon", "missing", "wrfin", "timeidx", "stagger", "projection", - "ll_lat", "ll_lon", - "pivot_point", "angle", + "ll_point", "pivot_point", "angle", "start_point", "end_point", "cache"), *args, **kwargs) @@ -902,8 +901,7 @@ def _set_cross_meta(wrapped, instance, args, kwargs): timeidx = argvars["timeidx"] stagger = argvars["stagger"] projection = argvars["projection"] - ll_lat = argvars["ll_lat"] - ll_lon = argvars["ll_lon"] + ll_point = argvars["ll_point"] pivot_point = argvars["pivot_point"] angle = argvars["angle"] start_point = argvars["start_point"] @@ -917,7 +915,7 @@ def _set_cross_meta(wrapped, instance, args, kwargs): 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, ll_lat, ll_lon) + stagger, projection, ll_point) pivot_point_xy = (xy_coords.x, xy_coords.y) else: pivot_point_xy = (pivot_point.x, pivot_point.y) @@ -925,14 +923,14 @@ def _set_cross_meta(wrapped, instance, args, kwargs): 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, ll_lat, ll_lon) + stagger, projection, ll_point) 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, ll_lat, ll_lon) + stagger, projection, ll_point) end_point_xy = (xy_coords.x, xy_coords.y) else: end_point_xy = (end_point.x, end_point.y) @@ -1107,8 +1105,7 @@ def _set_line_meta(wrapped, instance, args, kwargs): """ argvars = from_args(wrapped, ("field2d", "wrfin", "timeidx", "stagger", "projection", - "ll_lat", "ll_lon", - "pivot_point", "angle", + "ll_point", "pivot_point", "angle", "start_point", "end_point", "latlon", "cache"), *args, **kwargs) @@ -1118,8 +1115,7 @@ def _set_line_meta(wrapped, instance, args, kwargs): timeidx = argvars["timeidx"] stagger = argvars["stagger"] projection = argvars["projection"] - ll_lat = argvars["ll_lat"] - ll_lon = argvars["ll_lon"] + ll_point = argvars["ll_point"] pivot_point = argvars["pivot_point"] angle = argvars["angle"] start_point = argvars["start_point"] @@ -1137,7 +1133,7 @@ def _set_line_meta(wrapped, instance, args, kwargs): 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, ll_lat, ll_lon) + stagger, projection, ll_point) pivot_point_xy = (xy_coords.x, xy_coords.y) else: pivot_point_xy = (pivot_point.x, pivot_point.y) @@ -1146,14 +1142,14 @@ def _set_line_meta(wrapped, instance, args, kwargs): 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, ll_lat, ll_lon) + stagger, projection, ll_point) 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, ll_lat, ll_lon) + stagger, projection, ll_point) end_point_xy = (xy_coords.x, xy_coords.y) else: end_point_xy = (end_point.x, end_point.y) diff --git a/src/wrf/util.py b/src/wrf/util.py index 22cc2d6..3aa9583 100644 --- a/src/wrf/util.py +++ b/src/wrf/util.py @@ -34,6 +34,7 @@ from .constants import default_fill, ALL_TIMES from .py3compat import (viewitems, viewkeys, isstr, py3range, ucode) from .cache import cache_item, get_cached_item from .geobnds import GeoBounds, NullGeoBounds +from .coordpair import CoordPair from .projection import getproj @@ -3813,7 +3814,68 @@ def cartopy_ylim(var=None, geobounds=None, wrfin=None, varname=None, timeidx=0, return wrf_proj.cartopy_ylim(native_geobnds) +def ll_points(lat, lon): + """Return the lower left latitude and longitude point(s). + + This functions extracts the lower left corner points and returns the result + as either a single :class:`CoordPair` object, or a list of + :class:`CoordPair` objects. + + This is primarily used for testing or constructing the corner point objects + from the XLAT and XLONG variables. + + Args: + + lat (:class:`xarray.DataArray` or :class:`numpy.ndarray`): The latitude + array. Must be at least two dimensions. + + lon (:class:`xarray.DataArray` or :class:`numpy.ndarray`): The + longitude array. Must be at least two dimensions. + + Returns: + + :class:`wrf.CoordPair` or :obj:`list`: A single :class:`wrf.CoordPair` + object or a list of :class:`wrf.CoordPair` objects. + + """ + latvals = np.ravel(to_np(lat)[...,0,0]) + lonvals = np.ravel(to_np(lon)[...,0,0]) + + if latvals.shape[0] == 1: + return CoordPair(lat=float(latvals), lon=float(lonvals)) + else: + return [CoordPair(lat=latvals[i], lon=lonvals[i]) + for i in py3range(latvals.shape[0])] + +def pairs_to_latlon(pairs): + """Return latitude and longitude arrays from a sequence of \ + :class:`wrf.CoordPair` objects. + + This function converts a sequence of :class:`wrf.CoordPair` objects into + lists of latitude and longitude points. If the *pairs* argument is a + single :class:`wrf.CoordPair` object, then a single latitude and + longitude value is returned. + + Args: + + pairs (:class:`wrf.CoordPair` or sequence): A single + :class:`wrf.CoordPair` or sequence of :class:`wrf.CoordPair`. + + Returns: + + :obj:`tuple`: A tuple of (lat, lon), where lat and lon are single + values or lists of values. + + """ + + if isinstance(pairs, CoordPair): + return (pairs.lat, pairs.lon) + else: + lats = [pair.lat for pair in pairs] + lons = [pair.lon for pair in pairs] + + return lats, lons diff --git a/test/utests.py b/test/utests.py index d00fef8..a4e3013 100644 --- a/test/utests.py +++ b/test/utests.py @@ -8,7 +8,7 @@ import subprocess from wrf import (getvar, interplevel, interpline, vertcross, vinterp, disable_xarray, xarray_enabled, to_np, xy_to_ll, ll_to_xy, xy_to_ll_proj, ll_to_xy_proj, - extract_global_attrs, viewitems, CoordPair) + extract_global_attrs, viewitems, CoordPair, ll_points) from wrf.util import is_multi_file NCL_EXE = "/Users/ladwig/nclbuild/6.3.0/bin/ncl" @@ -299,14 +299,13 @@ def make_interp_test(varname, wrf_in, referent, multi=False, # Test the manual projection method with lat/lon lats = hts.coords["XLAT"] lons = hts.coords["XLONG"] - ll_lat = lats[0,0] - ll_lon = lons[0,0] + ll_point = ll_points(lats, lons) pivot = CoordPair(lat=lats[lats.shape[-2]/2, lats.shape[-1]/2], lon=lons[lons.shape[-2]/2, lons.shape[-1]/2]) v1 = vertcross(hts,p,wrfin=in_wrfnc,pivot_point=pivot_point, angle=90.0) v2 = vertcross(hts,p,projection=hts.attrs["projection"], - ll_lat=ll_lat, ll_lon=ll_lon, + ll_point=ll_point, pivot_point=pivot_point, angle=90.) nt.assert_allclose(to_np(v1), to_np(v2), rtol=.01) @@ -348,14 +347,13 @@ def make_interp_test(varname, wrf_in, referent, multi=False, # Test the manual projection method with lat/lon lats = t2.coords["XLAT"] lons = t2.coords["XLONG"] - ll_lat = lats[0,0] - ll_lon = lons[0,0] + ll_point = ll_points(lats, lons) pivot = CoordPair(lat=lats[lats.shape[-2]/2, lats.shape[-1]/2], lon=lons[lons.shape[-2]/2, lons.shape[-1]/2]) l1 = interpline(t2,wrfin=in_wrfnc,pivot_point=pivot_point, angle=90.0) l2 = interpline(t2,projection=t2.attrs["projection"], - ll_lat=ll_lat, ll_lon=ll_lon, + ll_point=ll_point, pivot_point=pivot_point, angle=90.) nt.assert_allclose(to_np(l1), to_np(l2), rtol=.01)