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

@@ -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
  
  
 ------------------------

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"]

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))
     
-
+    
     

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)

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, 

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

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)

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
         
             
 

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)