Fix issues with moving nests and line interpolation.

Updated unit tests.

src/wrf/interp.py

@@ -106,52 +106,6 @@ def interplevel(field3d, vert, desiredlev, missing=default_fill(np.float64),
     return masked
 
 
-def _vertcross_nest_alltimes(field3d, vert, levels=None, 
-              missing=default_fill(np.float64),
-              wrfin=None, timeidx=0, stagger=None, projection=None,
-              ll_point=None,
-              pivot_point=None, angle=None,
-              start_point=None, end_point=None,
-              latlon=False, autolevels=100, cache=None, meta=True):
-    
-    # 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
-    
-    # Check if we have a wrfin file, or this is a no go.
-    if wrfin is None:
-        raise ValueError("'wrfin' is required when using all times "
-                         "from a moving nest with lat/lon coords")
-        
-    if multi:
-        if field3d.ndim == 4:
-            raise ValueError("all times requested for a moving nest, "
-                             "but no time dimension found for "
-                             'field3d')
-    else:
-        if field3d.ndim < 4:
-            raise ValueError("all times requested for a moving nest, "
-                             "but no time dimension found for "
-                             'field3d')
-            
-    numtimes = field3d.shape[-4]
-    
-    for t in py3range(numtimes):
-        #_meta = True if t == 0 else False
-        _meta = True
-        
-        v = vertcross(field3d, vert, levels, missing, wrfin, t, stagger, 
-                      projection, ll_point, pivot_point, angle, start_point, 
-                      end_point, latlon, autolevels, cache, _meta)
-        
-        print(v.attrs)
-        
-        
-        
-        
-    
-
-
 @set_interp_metadata("cross")
 def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64),
               wrfin=None, timeidx=0, stagger=None, projection=None,
@@ -308,6 +262,15 @@ def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64),
     # type, we'll handle that iteration here.
     multi = True if field3d.ndim - vert.ndim == 1 else False
     
+    try:
+        xy = cache["xy"]
+        var2dz = cache["var2dz"]
+        z_var2d = cache["z_var2d"]
+    except (KeyError, TypeError):
+        # Convert Lat/Lon points to grid points
+        start_point_xy = None
+        end_point_xy = None
+        pivot_point_xy = None
         
         if timeidx is None:
             if (latlon is True or is_latlon_pair(start_point) or 
@@ -334,17 +297,6 @@ def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64),
         else:
             _timeidx = timeidx
         
-    
-    try:
-        xy = cache["xy"]
-        var2dz = cache["var2dz"]
-        z_var2d = cache["z_var2d"]
-    except (KeyError, TypeError):
-        # 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, 
@@ -512,6 +464,14 @@ def interpline(field2d, pivot_point=None,
         
         
     """
+    
+    try:
+        xy = cache["xy"]
+    except (KeyError, TypeError):
+        start_point_xy = None
+        end_point_xy = None
+        pivot_point_xy = None
+        
         if timeidx is None:
             if (latlon is True or is_latlon_pair(start_point) or 
                 is_latlon_pair(pivot_point)):
@@ -528,6 +488,7 @@ def interpline(field2d, pivot_point=None,
                                          "You must request each time "
                                          "individually.")
                     else:
+                        # Domain not moving, just use 0
                         _timeidx = 0
                 
             # If using grid coordinates, then don't care about lat/lon
@@ -537,13 +498,6 @@ def interpline(field2d, pivot_point=None,
         else:
             _timeidx = timeidx
         
-    try:
-        xy = cache["xy"]
-    except (KeyError, TypeError):
-        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, 

src/wrf/metadecorators.py

@@ -9,7 +9,8 @@ import numpy.ma as ma
 from .extension import _interpline
 from .util import (extract_vars, either, from_args, arg_location,
                    is_coordvar, latlon_coordvars, to_np, 
-                   from_var, iter_left_indexes, is_mapping)
+                   from_var, iter_left_indexes, is_mapping,
+                   is_moving_domain, is_latlon_pair)
 from .coordpair import CoordPair
 from .py3compat import viewkeys, viewitems, py3range
 from .interputils import get_xy_z_params, get_xy, to_xy_coords
@@ -932,9 +933,34 @@ def _set_cross_meta(wrapped, instance, args, kwargs):
     end_point_xy = None
     pivot_point_xy = None
     
+    if timeidx is None:
+        if (inc_latlon is True or is_latlon_pair(start_point) or 
+            is_latlon_pair(pivot_point)):
+            
+            if wrfin is not None:
+                # Moving nests aren't supported with ALL_TIMES because the 
+                # domain could move outside of the cross section, which causes 
+                # crashes or different line lengths.
+                if is_moving_domain(wrfin):
+                    raise ValueError("Requesting all times with a moving nest "
+                                     "is not supported when using lat/lon "
+                                     "cross sections because the domain could "
+                                     "move outside of the cross section. "
+                                     "You must request each time "
+                                     "individually.")
+                else:
+                    _timeidx = 0
+            
+        # If using grid coordinates, then don't care about lat/lon
+        # coordinates. Just use 0.
+        else:
+            _timeidx = 0
+    else:
+        _timeidx = timeidx
+    
     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, 
+            xy_coords = to_xy_coords(pivot_point, wrfin, _timeidx, 
                                      stagger, projection, ll_point)
             pivot_point_xy = (xy_coords.x, xy_coords.y)
         else:
@@ -942,14 +968,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, 
+            xy_coords = to_xy_coords(start_point, wrfin, _timeidx, 
                                      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, 
+            xy_coords = to_xy_coords(end_point, wrfin, _timeidx, 
                                      stagger, projection, ll_point)
             end_point_xy = (xy_coords.x, xy_coords.y)
         else:
@@ -1150,9 +1176,35 @@ def _set_line_meta(wrapped, instance, args, kwargs):
     end_point_xy = None
     pivot_point_xy = None
     
+    if timeidx is None:
+        if (inc_latlon is True or is_latlon_pair(start_point) or 
+            is_latlon_pair(pivot_point)):
+            
+            if wrfin is not None:
+                # Moving nests aren't supported with ALL_TIMES because the 
+                # domain could move outside of the line, which causes 
+                # crashes or different line lengths.
+                if is_moving_domain(wrfin):
+                    raise ValueError("Requesting all times with a moving nest "
+                                     "is not supported when using a lat/lon "
+                                     "line because the domain could "
+                                     "move outside of line. "
+                                     "You must request each time "
+                                     "individually.")
+                else:
+                    # Domain not moving, just use 0
+                    _timeidx = 0
+            
+        # If using grid coordinates, then don't care about lat/lon
+        # coordinates. Just use 0.
+        else:
+            _timeidx = 0
+    else:
+        _timeidx = timeidx
+        
     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, 
+            xy_coords = to_xy_coords(pivot_point, wrfin, _timeidx, 
                                      stagger, projection, ll_point)
             pivot_point_xy = (xy_coords.x, xy_coords.y)
         else:
@@ -1161,19 +1213,20 @@ 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, 
+            xy_coords = to_xy_coords(start_point, wrfin, _timeidx, 
                                      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, 
+            xy_coords = to_xy_coords(end_point, wrfin, _timeidx, 
                                      stagger, projection, ll_point)
             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

test/ncl_get_var.ncl

@@ -236,18 +236,20 @@
     
     t2_line2 = wrf_user_interpline(t2, pivot, xopt)
     
+    fout->t2_line2 = t2_line2
+    
     lats = wrf_user_getvar(input_file, "lat", 0)
     lons = wrf_user_getvar(input_file, "lon", 0) 
         
     start_lat = min(lats) + .25d*(max(lats) - min(lats))
-    end_lat = min(lats) + .75d*(max(lats) - min(lats))
+    end_lat = min(lats) + .65d*(max(lats) - min(lats))
                 
     start_lon = min(lons) + .25d*(max(lons) - min(lons))
-    end_lon = min(lons) + .75d*(max(lons) - min(lons))
+    end_lon = min(lons) + .65d*(max(lons) - min(lons))
         
     start_end = (/ start_lon, start_lat, end_lon, end_lat /)
     
-    ; For the new cross section routine
+    ; For the new line routine
     xopt := True
     xopt@use_pivot = False
     xopt@latlon = True
@@ -256,6 +258,9 @@
     xopt@linecoords = True
     
     t2_line3 = wrf_user_interpline(t2, start_end, xopt)
+    t2_line3!1 = "line_idx_t2_line3"
+    
+    fout->t2_line3 = t2_line3
     
     times = wrf_user_getvar(input_file, "times", -1)
     ntimes = dimsizes(times)
@@ -270,6 +275,22 @@
       fout->$name$ = t2_line
     end do
     
+    ; Make sure the 1 time case still works
+    t2 := wrf_user_getvar(input_file, "T2", 0)
+    
+    ; For the new line routine
+    xopt := True
+    xopt@use_pivot = False
+    xopt@latlon = True
+    xopt@file_handle = input_file
+    xopt@timeidx = 0
+    xopt@linecoords = True
+    
+    t2_line4 = wrf_user_interpline(t2, start_end, xopt)
+    t2_line4!0 = "t2_line4_idx"
+    
+    fout->t2_line4 = t2_line4
+    
     
     ;;;;;;;;;;;;;;;;;;;;;;; 3D interpolation to new vertical coordinates
     time = -1

test/utests.py

@@ -197,6 +197,7 @@ def _get_refvals(referent, varname, multi):
     
     multi2d = ("uvmet10", "cape_2d", "cfrac")
     multi3d = ("uvmet", "cape_3d")
+    interpline = ("t2_line","t2_line2", "t2_line3")
     
     if not multi:
         if varname in multi2d:
@@ -206,6 +207,9 @@ def _get_refvals(referent, varname, multi):
         else:
             v = refnc.variables[varname][:]
             if v.ndim == 2:
+                if varname in interpline:
+                    ref_vals = v[0,:]
+                else:
                     ref_vals = v
             else:
                 ref_vals = v[0,:]
@@ -477,6 +481,8 @@ def make_interp_test(varname, dir, pattern, referent, multi=False,
         elif (varname == "interpline"):
             
             ref_t2_line = _get_refvals(referent, "t2_line", multi)
+            ref_t2_line2 = _get_refvals(referent, "t2_line2", multi)
+            ref_t2_line3 = _get_refvals(referent, "t2_line3", multi)
             
             t2 = getvar(wrfin, "T2", timeidx=timeidx)
             pivot_point = CoordPair(t2.shape[-1] / 2, t2.shape[-2] / 2)
@@ -488,6 +494,9 @@ def make_interp_test(varname, dir, pattern, referent, multi=False,
             
             nt.assert_allclose(to_np(t2_line1), ref_t2_line)
             
+            # Test the new NCL wrf_user_interplevel result
+            nt.assert_allclose(to_np(t2_line1), ref_t2_line2)
+            
             # Test the manual projection method with lat/lon
             lats = t2.coords["XLAT"]
             lons = t2.coords["XLONG"]
@@ -516,6 +525,63 @@ def make_interp_test(varname, dir, pattern, referent, multi=False,
                                   end_point=end_point)
             
             nt.assert_allclose(to_np(t2_line1), to_np(t2_line2))
+            
+            # Now test the start/end with lat/lon points
+            
+            start_lat = float(np.amin(lats) + .25*(np.amax(lats) 
+                                                   - np.amin(lats)))
+            end_lat = float(np.amin(lats) + .65*(np.amax(lats) 
+                                                 - np.amin(lats)))
+            
+            start_lon = float(np.amin(lons) + .25*(np.amax(lons) 
+                                                   - np.amin(lons)))
+            end_lon = float(np.amin(lons) + .65*(np.amax(lons) 
+                                                 - np.amin(lons)))
+            
+            start_point = CoordPair(lat=start_lat, lon=start_lon)
+            end_point = CoordPair(lat=end_lat, lon=end_lon)
+            
+            t2_line3 = interpline(t2,wrfin=wrfin,timeidx=0,
+                                  start_point=start_point,
+                                  end_point=end_point,latlon=True)
+            
+            
+            nt.assert_allclose(to_np(t2_line3), ref_t2_line3, rtol=.01)
+            
+            # Test all time steps
+            if multi:
+                refnc = NetCDF(referent)
+                ntimes = t2.shape[0]
+                
+                for t in range(ntimes):
+                    t2 = getvar(wrfin, "T2", timeidx=t)
+                    
+                    line = interpline(t2,wrfin=wrfin,timeidx=t,
+                                  start_point=start_point,
+                                  end_point=end_point,latlon=True)
+                    
+                    refname = "t2_line_t{}".format(t)
+                    refline = refnc.variables[refname][:]
+                    
+                    nt.assert_allclose(to_np(line), 
+                               to_np(refline),rtol=.005)
+            
+                refnc.close()
+            
+            # Test NCLs single time case
+            if not multi:
+                refnc = NetCDF(referent)
+                ref_t2_line4 = refnc.variables["t2_line4"][:]
+                
+                t2 = getvar(wrfin, "T2", timeidx=0)
+                line = interpline(t2,wrfin=wrfin,timeidx=0,
+                                  start_point=start_point,
+                                  end_point=end_point,latlon=True)
+                
+                nt.assert_allclose(to_np(line), 
+                               to_np(ref_t2_line4),rtol=.005)
+                refnc.close()
+
         elif (varname == "vinterp"):
             # Tk to theta
             fld_tk_theta = _get_refvals(referent, "fld_tk_theta", multi)
@@ -844,9 +910,9 @@ class WRFLatLonTest(ut.TestCase):
 if __name__ == "__main__":
     from wrf import (omp_set_num_threads, omp_set_schedule, omp_get_schedule, 
                      omp_set_dynamic, omp_get_num_procs, OMP_SCHED_STATIC)
-    #omp_set_num_threads(omp_get_num_procs()//2)
-    #omp_set_schedule(OMP_SCHED_STATIC, 0)
-    #omp_set_dynamic(False)
+    omp_set_num_threads(omp_get_num_procs()//2)
+    omp_set_schedule(OMP_SCHED_STATIC, 0)
+    omp_set_dynamic(False)
     
     ignore_vars = []  # Not testable yet
     wrf_vars = ["avo", "eth", "cape_2d", "cape_3d", "ctt", "dbz", "mdbz",