Merge branch 'bugfix/new_ncl' into develop

src/wrf/interp.py

@@ -94,7 +94,7 @@ def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64),
               ll_point=None,
               pivot_point=None, angle=None,
               start_point=None, end_point=None,
-              latlon=False, cache=None, meta=True):
+              latlon=False, autolevels=100, cache=None, meta=True):
     """Return the vertical cross section for a three-dimensional field.
     
     The cross section is defined by a horizontal line through the domain.  
@@ -213,6 +213,10 @@ def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64),
                 *latlon* is set to True. Otherwise, a warning will be issued,
                 and the latitude and longitude information will not be 
                 present.
+                
+        autolevels(:obj:`int`, optional): The number of evenly spaced 
+            automatically chosen vertical levels to use when *levels* 
+            is None. Default is 100.
              
         cache (:obj:`dict`, optional): A dictionary of (varname, ndarray) 
             that can be used to supply pre-extracted NetCDF variables to the 
@@ -279,7 +283,8 @@ def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64),
                 
         xy, var2dz, z_var2d = get_xy_z_params(to_np(vert), pivot_point_xy, 
                                               angle, start_point_xy, 
-                                              end_point_xy, levels)
+                                              end_point_xy, levels,
+                                              autolevels)
     
     if not multi:
         result = _vertcross(field3d, xy, var2dz, z_var2d, missing)

src/wrf/interputils.py

@@ -185,7 +185,7 @@ def _calc_xy(xdim, ydim, pivot_point=None, angle=None,
 
 def get_xy_z_params(z, pivot_point=None, angle=None,
                     start_point=None, end_point=None,
-                    levels=None):
+                    levels=None, autolevels=100):
     """Return the cross section parameters.
     
     This function returns the xy horizontal cross section line coordinates, 
@@ -224,6 +224,10 @@ def get_xy_z_params(z, pivot_point=None, angle=None,
             vertical levels in the output array.  If None, a fixed set of 
             vertical levels is provided.  Default is None.
             
+        autolevels(:obj:`int`, optional): The number of evenly spaced 
+            automatically chosen vertical levels to use when *levels* 
+            is None. Default is 100.
+            
     Returns:
     
         :obj:`tuple`:  A tuple containing the xy horizontal cross section  
@@ -248,20 +252,18 @@ def get_xy_z_params(z, pivot_point=None, angle=None,
         if(var2dz[idx1] > var2dz[idx2]):  # monotonically decreasing coordinate
             z_max = floor(np.amax(z)/10) * 10     # bottom value
             z_min = ceil(np.amin(z)/10) * 10      # top value
-            dz = 10
-            nlevels = int((z_max - z_min)/dz)
-            z_var2d = np.zeros((nlevels), dtype=z.dtype)
+            dz = (1.0/autolevels) * (z_max - z_min)
+            z_var2d = np.zeros((autolevels), dtype=z.dtype)
             z_var2d[0] = z_max
             dz = -dz
         else:
             z_max = np.amax(z)
             z_min = 0.
-            dz = 0.01*z_max
-            nlevels = int(z_max/dz)
-            z_var2d = np.zeros((nlevels), dtype=z.dtype)
+            dz = (1.0/autolevels)*z_max
+            z_var2d = np.zeros((autolevels), dtype=z.dtype)
             z_var2d[0] = z_min
             
-        for i in py3range(1,nlevels):
+        for i in py3range(1,autolevels):
             z_var2d[i] = z_var2d[0] + i*dz
     else:
         z_var2d = np.asarray(levels, z.dtype)

src/wrf/metadecorators.py

@@ -889,7 +889,7 @@ def _set_cross_meta(wrapped, instance, args, kwargs):
                                   "latlon", "missing", 
                                   "wrfin", "timeidx", "stagger", "projection",
                                   "ll_point", "pivot_point", "angle",
-                                  "start_point", "end_point",
+                                  "start_point", "end_point", "autolevels",
                                   "cache"), 
                           *args, **kwargs)  
     
@@ -907,6 +907,7 @@ def _set_cross_meta(wrapped, instance, args, kwargs):
     angle = argvars["angle"]
     start_point = argvars["start_point"]
     end_point = argvars["end_point"]
+    autolevels = argvars["autolevels"]
     cache = argvars["cache"]
     
     start_point_xy = None
@@ -937,7 +938,7 @@ def _set_cross_meta(wrapped, instance, args, kwargs):
             end_point_xy = (end_point.x, end_point.y)
     
     xy, var2dz, z_var2d = get_xy_z_params(to_np(z), pivot_point_xy, angle,
-              start_point_xy, end_point_xy, levels)
+              start_point_xy, end_point_xy, levels, autolevels)
     
     # Make a copy so we don't modify a user supplied cache
     if cache is not None:

test/ncl_get_var.ncl

@@ -5,7 +5,7 @@
     ;system("printenv")
     
     if (.not. isvar("in_file")) then
-        in_file = "/Users/ladwig/Documents/wrf_files/wrfout_d02_2010-06-13_21:00:00.nc"
+        in_file = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2010-06-13_21:00:00.nc"
     end if 
     
     if (.not. isvar("out_file")) then
@@ -48,20 +48,21 @@
     
     ; Do the interpolation routines manually
     ; 3D vertical cross section
-    z  = wrf_user_getvar(input_file, "z", 0)        ; grid point height
-    p  = wrf_user_getvar(input_file, "pressure", 0) ; total pressure
+    z := wrf_user_getvar(input_file, "z", 0)        ; grid point height
+    p := wrf_user_getvar(input_file, "pressure", 0) ; total pressure
     
     dimsz = dimsizes(z)
     pivot = (/ dimsz(2)/2, dimsz(1)/2 /)    ; pivot point is center of domain
     
     ht_cross = wrf_user_intrp3d(z,p,"v",pivot,90.0,False)
+    
     p_cross = wrf_user_intrp3d(p,z,"v",pivot,90.0,False)
     p_cross!0 = "Vertical_p"
     p_cross!1 = "Horizontal_p"
     
     fout->ht_cross = ht_cross
     fout->p_cross = p_cross
-    
+;    
     ; For the new cross section routine
     xopt = True
     xopt@use_pivot = True
@@ -90,17 +91,17 @@
     
     fout->ht_vertcross2 = ht_vertcross2
     
-    ; Note for the default file:
-    ; >>> np.amin(lats)
-    ; 32.77895
-    ; >>> np.amax(lats)
-    ; 40.209736
-    ; >>> np.amin(lons)
-    ; -104.71738
-    ; >>> np.amax(lons)
-    ; -95.153076
-    
-    start_end = (/ -102.0, 35.0, -98.0, 37.0 /)
+    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))
+            
+    start_lon = min(lons) + .25d*(max(lons) - min(lons))
+    end_lon = min(lons) + .75d*(max(lons) - min(lons))
+    
+    start_end = (/ start_lon, start_lat, end_lon, end_lat /)
+    
     ; For the new cross section routine
     xopt := True
     xopt@use_pivot = False
@@ -108,6 +109,7 @@
     xopt@file_handle = input_file
     xopt@timeidx = 0
     xopt@linecoords = True
+    xopt@autolevels = 1000
     
     ht_vertcross3 = wrf_user_vertcross(z, p, start_end, xopt)
     ht_vertcross3!0 = "vertical3"
@@ -118,7 +120,6 @@
     ; 3D horizontal interpolation
     plev = 500.   ; 500 MB
     z_500 = wrf_user_intrp3d(z,p,"h",plev,0.,False)
-    ;z_500_dims = dimsizes(z_500)
     
     fout->z_500 = z_500
     
@@ -221,8 +222,8 @@
     
     
     ; lat/lon to x/y and x/y to lat/lon routines
-    lats = (/-55, -60, -65 /)
-    lons = (/25, 30, 35 /)
+    lats := (/-55, -60, -65 /)
+    lons := (/25, 30, 35 /)
     i_s = (/10, 100, 150 /)
     j_s = (/10, 100, 150 /)
     ij = wrf_user_ll_to_ij(input_file, lons, lats, True) 

test/utests.py

@@ -65,6 +65,10 @@ def make_test(varname, wrf_in, referent, multi=False, repeat=3, pynio=False):
             timeidx = 0
             if not pynio:
                 in_wrfnc = NetCDF(wrf_in)
+                try:
+                    in_wrfnc.set_auto_mask(False)
+                except:
+                    pass
             else:
                 # Note: Python doesn't like it if you reassign an outer scoped
                 # variable (wrf_in in this case)
@@ -77,6 +81,10 @@ def make_test(varname, wrf_in, referent, multi=False, repeat=3, pynio=False):
             timeidx = None
             if not pynio:
                 nc = NetCDF(wrf_in)
+                try:
+                    nc.set_auto_mask(False)
+                except:
+                    pass
                 in_wrfnc = [nc for i in xrange(repeat)]
             else:
                 if not wrf_in.endswith(".nc"):
@@ -88,6 +96,10 @@ def make_test(varname, wrf_in, referent, multi=False, repeat=3, pynio=False):
             
         
         refnc = NetCDF(referent)
+        try:
+            refnc.set_auto_mask(False)
+        except:
+            pass
         
         # These have a left index that defines the product type
         multiproduct = varname in ("uvmet", "uvmet10", "cape_2d", "cape_3d", 
@@ -141,10 +153,15 @@ def make_test(varname, wrf_in, referent, multi=False, repeat=3, pynio=False):
             tol = 1/100.
             atol = .1 # Note:  NCL uses 273.16 as conversion for some reason
             nt.assert_allclose(to_np(my_vals), ref_vals, tol, atol)
+        elif (varname == "height_agl"):
+            # Change the vert_type to height_agl when NCL gets updated.
+            my_vals = getvar(in_wrfnc, "z", timeidx=timeidx, msl=False)
+            tol = 1/100.
+            atol = .1 # Note:  NCL uses 273.16 as conversion for some reason
+            nt.assert_allclose(to_np(my_vals), ref_vals, tol, atol)
         elif (varname == "cfrac"):
             # Change the vert_type to height_agl when NCL gets updated.
-            my_vals = getvar(in_wrfnc, "cfrac", timeidx=timeidx, 
-                             vert_type="pres")
+            my_vals = getvar(in_wrfnc, "cfrac", timeidx=timeidx)
             tol = 1/100.
             atol = .1 # Note:  NCL uses 273.16 as conversion for some reason
             nt.assert_allclose(to_np(my_vals), ref_vals, tol, atol)
@@ -191,7 +208,11 @@ def make_test(varname, wrf_in, referent, multi=False, repeat=3, pynio=False):
             try:
                 nt.assert_allclose(to_np(my_vals), ref_vals, tol, atol)
             except:
-                print (np.amax(np.abs(to_np(my_vals) - ref_vals)))
+                absdiff = np.abs(to_np(my_vals) - ref_vals)
+                maxdiff = np.amax(absdiff)
+                print (maxdiff)
+                print np.argwhere(absdiff == maxdiff)
+                
                 raise
     
     
@@ -204,6 +225,10 @@ def _get_refvals(referent, varname, repeat, multi):
         pass
         
     refnc = NetCDF(referent)
+    try:
+        refnc.set_auto_mask(False)
+    except:
+        pass
     
     if not multi:
         ref_vals = refnc.variables[varname][:]
@@ -244,6 +269,10 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
             timeidx = 0
             if not pynio:
                 in_wrfnc = NetCDF(wrf_in)
+                try:
+                    in_wrfnc.set_auto_mask(False)
+                except:
+                    pass
             else:
                 # Note: Python doesn't like it if you reassign an outer scoped
                 # variable (wrf_in in this case)
@@ -256,6 +285,10 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
             timeidx = None
             if not pynio:
                 nc = NetCDF(wrf_in)
+                try:
+                    nc.set_auto_mask(False)
+                except:
+                    pass
                 in_wrfnc = [nc for i in xrange(repeat)]
             else:
                 if not wrf_in.endswith(".nc"):
@@ -279,22 +312,35 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
         elif (varname == "vertcross"):
             ref_ht_cross = _get_refvals(referent, "ht_cross", repeat, multi)
             ref_p_cross = _get_refvals(referent, "p_cross", repeat, multi)
+            ref_ht_vertcross1 = _get_refvals(referent, "ht_vertcross1", repeat, 
+                                            multi)
+            ref_ht_vertcross2 = _get_refvals(referent, "ht_vertcross2", repeat, 
+                                            multi)
+            ref_ht_vertcross3 = _get_refvals(referent, "ht_vertcross3", repeat, 
+                                            multi)
             
             hts = getvar(in_wrfnc, "z", timeidx=timeidx)
             p = getvar(in_wrfnc, "pressure", timeidx=timeidx)
             
             pivot_point = CoordPair(hts.shape[-1] / 2, hts.shape[-2] / 2) 
-            #ht_cross = vertcross(to_np(hts), p, pivot_point=pivot_point, 
-            #                     angle=90., latlon=True)
+            
+            # Beginning in wrf-python 1.3, users can select number of levels.
+            # Originally, for pressure, dz was 10, so let's back calculate
+            # the number of levels.
+            p_max = np.floor(np.amax(p)/10) * 10     # bottom value
+            p_min = np.ceil(np.amin(p)/10) * 10      # top value
+            
+            p_autolevels = int((p_max - p_min) /10)
+
             # Make sure the numpy versions work first
+            
             ht_cross = vertcross(to_np(hts), to_np(p), 
-                                 pivot_point=pivot_point, angle=90.)
-            ht_cross = vertcross(hts, p, pivot_point=pivot_point, angle=90.)
+                                 pivot_point=pivot_point, angle=90.,
+                                 autolevels=p_autolevels)
+            ht_cross = vertcross(hts, p, pivot_point=pivot_point, angle=90.,
+                                 autolevels=p_autolevels)
             
-            # Note:  Until the bug is fixed in NCL, the wrf-python cross 
-            # sections will contain one extra point
-            nt.assert_allclose(to_np(ht_cross)[...,0:-1], ref_ht_cross, 
-                               rtol=.01)
+            nt.assert_allclose(to_np(ht_cross), ref_ht_cross, rtol=.01)
             
             # Test the manual projection method with lat/lon
             lats = hts.coords["XLAT"]
@@ -304,6 +350,7 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
                                        int(lats.shape[-1]/2)],
                               lon=lons[int(lons.shape[-2]/2), 
                                        int(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"], 
@@ -313,9 +360,10 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
             nt.assert_allclose(to_np(v1), to_np(v2), rtol=.01)
             
             # Test opposite
+            
             p_cross1 = vertcross(p,hts,pivot_point=pivot_point, angle=90.0)
             
-            nt.assert_allclose(to_np(p_cross1)[...,0:-1], 
+            nt.assert_allclose(to_np(p_cross1), 
                                ref_p_cross, 
                                rtol=.01)
             # Test point to point
@@ -328,6 +376,50 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
             
             nt.assert_allclose(to_np(p_cross1), 
                                to_np(p_cross2))
+            
+            # Check the new vertcross routine
+            pivot_point = CoordPair(hts.shape[-1] / 2, hts.shape[-2] / 2) 
+            ht_cross = vertcross(hts, p, 
+                                 pivot_point=pivot_point, angle=90.,
+                                 latlon=True)
+            
+            nt.assert_allclose(to_np(ht_cross), 
+                               to_np(ref_ht_vertcross1), atol=.01)
+            
+            
+            
+            levels = [1000., 850., 700., 500., 250.]
+            ht_cross = vertcross(hts, p, 
+                                 pivot_point=pivot_point, angle=90.,
+                                 levels=levels, latlon=True)
+            
+            nt.assert_allclose(to_np(ht_cross), 
+                               to_np(ref_ht_vertcross2), atol=.01)
+            
+            
+            start_lat = np.amin(lats) + .25*(np.amax(lats) - np.amin(lats))
+            end_lat = np.amin(lats) + .75*(np.amax(lats) - np.amin(lats))
+            
+            start_lon = np.amin(lons) + .25*(np.amax(lons) - np.amin(lons))
+            end_lon = np.amin(lons) + .75*(np.amax(lons) - np.amin(lons))
+            
+            start_point = CoordPair(lat=start_lat, lon=start_lon)
+            end_point = CoordPair(lat=end_lat, lon=end_lon)
+            
+            # ll_point and projection came from above 
+            ht_cross = vertcross(hts, p, 
+                                 start_point=start_point, 
+                                 end_point=end_point, 
+                                 projection=hts.attrs["projection"], 
+                                 ll_point=ll_point, 
+                                 latlon=True,
+                                 autolevels=1000)
+            
+            
+            nt.assert_allclose(to_np(ht_cross), 
+                               to_np(ref_ht_vertcross3),
+                               rtol=.01)
+            
               
         elif (varname == "interpline"):
             
@@ -336,15 +428,12 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
             t2 = getvar(in_wrfnc, "T2", timeidx=timeidx)
             pivot_point = CoordPair(t2.shape[-1] / 2, t2.shape[-2] / 2)
             
-            #t2_line1 = interpline(to_np(t2), pivot_point=pivot_point, 
-            #                      angle=90.0, latlon=True)
             # Make sure the numpy version works
             t2_line1 = interpline(to_np(t2), pivot_point=pivot_point, 
                                   angle=90.0)
             t2_line1 = interpline(t2, pivot_point=pivot_point, angle=90.0)
             
-            # Note:  After NCL is fixed, remove the slice.
-            nt.assert_allclose(to_np(t2_line1)[...,0:-1], ref_t2_line)
+            nt.assert_allclose(to_np(t2_line1), ref_t2_line)
             
             # Test the manual projection method with lat/lon
             lats = t2.coords["XLAT"]
@@ -575,6 +664,10 @@ def make_latlon_test(testid, wrf_in, referent, single, multi=False, repeat=3,
             timeidx = 0
             if not pynio:
                 in_wrfnc = NetCDF(wrf_in)
+                try:
+                    in_wrfnc.set_auto_mask(False)
+                except:
+                    pass
             else:
                 # Note: Python doesn't like it if you reassign an outer scoped
                 # variable (wrf_in in this case)
@@ -587,6 +680,10 @@ def make_latlon_test(testid, wrf_in, referent, single, multi=False, repeat=3,
             timeidx = None
             if not pynio:
                 nc = NetCDF(wrf_in)
+                try:
+                    nc.set_auto_mask(False)
+                except:
+                    pass
                 in_wrfnc = [nc for i in xrange(repeat)]
             else:
                 if not wrf_in.endswith(".nc"):
@@ -597,6 +694,10 @@ def make_latlon_test(testid, wrf_in, referent, single, multi=False, repeat=3,
                 in_wrfnc = [nc for i in xrange(repeat)]
                 
         refnc = NetCDF(referent)
+        try:
+            refnc.set_auto_mask(False)
+        except:
+            pass
                 
         if testid == "xy":
             # Since this domain is not moving, the reference values are the 
@@ -691,7 +792,8 @@ if __name__ == "__main__":
                 "geopt", "helicity", "lat", "lon", "omg", "p", "pressure", 
                 "pvo", "pw", "rh2", "rh", "slp", "ter", "td2", "td", "tc", 
                 "theta", "tk", "tv", "twb", "updraft_helicity", "ua", "va", 
-                "wa", "uvmet10", "uvmet", "z", "cfrac", "zstag", "geopt_stag"]
+                "wa", "uvmet10", "uvmet", "z", "cfrac", "zstag", "geopt_stag",
+                "height_agl"]
     interp_methods = ["interplevel", "vertcross", "interpline", "vinterp"]
     latlon_tests = ["xy", "ll"]