From a3880b1e10886fb7b95db73b1f3a781e666792e4 Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Fri, 26 Oct 2018 16:30:20 -0600 Subject: [PATCH] Added autolevels parameter for vertcross Users can now specify how many autogenerated levels they want Updated unit tests with NCL for the new routines. --- src/wrf/interp.py | 9 ++- src/wrf/interputils.py | 18 +++--- src/wrf/metadecorators.py | 5 +- test/ncl_get_var.ncl | 37 ++++++------ test/utests.py | 121 +++++++++++++++++++++++++++++++++----- 5 files changed, 146 insertions(+), 44 deletions(-) diff --git a/src/wrf/interp.py b/src/wrf/interp.py index ec6e9d2..424587a 100755 --- a/src/wrf/interp.py +++ b/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) diff --git a/src/wrf/interputils.py b/src/wrf/interputils.py index 44a6ab7..04807d1 100644 --- a/src/wrf/interputils.py +++ b/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) diff --git a/src/wrf/metadecorators.py b/src/wrf/metadecorators.py index 63e1195..e406b6e 100644 --- a/src/wrf/metadecorators.py +++ b/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: diff --git a/test/ncl_get_var.ncl b/test/ncl_get_var.ncl index 602e393..5db6a14 100644 --- a/test/ncl_get_var.ncl +++ b/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) diff --git a/test/utests.py b/test/utests.py index ed1d452..41d6201 100644 --- a/test/utests.py +++ b/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", @@ -204,6 +216,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 +260,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 +276,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 +303,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 +341,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 +351,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 +367,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 +419,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 +655,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 +671,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 +685,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 +783,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"]