From 4eeeadc9ce0bdcb85b5d4daaf4aed469b5af9dc2 Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Mon, 3 Dec 2018 17:09:30 -0700 Subject: [PATCH] Fix issues with moving nests and line interpolation. Updated unit tests. --- src/wrf/interp.py | 152 +++++++++++++------------------------- src/wrf/metadecorators.py | 67 +++++++++++++++-- test/ncl_get_var.ncl | 29 +++++++- test/utests.py | 74 ++++++++++++++++++- 4 files changed, 208 insertions(+), 114 deletions(-) diff --git a/src/wrf/interp.py b/src/wrf/interp.py index bc1e89c..81e5eb1 100755 --- a/src/wrf/interp.py +++ b/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,33 +262,6 @@ 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 - - if timeidx is None: - if (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 - - try: xy = cache["xy"] var2dz = cache["var2dz"] @@ -345,6 +272,31 @@ def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64), 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)): + + 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, @@ -512,30 +464,6 @@ def interpline(field2d, pivot_point=None, """ - if timeidx is None: - if (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: - _timeidx = 0 - - # If using grid coordinates, then don't care about lat/lon - # coordinates. Just use 0. - else: - _timeidx = 0 - else: - _timeidx = timeidx try: xy = cache["xy"] @@ -544,6 +472,32 @@ def interpline(field2d, pivot_point=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)): + + 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, @@ -566,10 +520,10 @@ def interpline(field2d, pivot_point=None, 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) - + return _interpline(field2d, xy) diff --git a/src/wrf/metadecorators.py b/src/wrf/metadecorators.py index 6957713..c8ff1b0 100644 --- a/src/wrf/metadecorators.py +++ b/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: @@ -1149,10 +1175,36 @@ def _set_line_meta(wrapped, instance, args, kwargs): start_point_xy = None 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 diff --git a/test/ncl_get_var.ncl b/test/ncl_get_var.ncl index bd4d3a6..a7b8aa5 100644 --- a/test/ncl_get_var.ncl +++ b/test/ncl_get_var.ncl @@ -236,26 +236,31 @@ 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 xopt@file_handle = input_file xopt@timeidx = 0 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 diff --git a/test/utests.py b/test/utests.py index 5310a2e..64f54fb 100644 --- a/test/utests.py +++ b/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,7 +207,10 @@ def _get_refvals(referent, varname, multi): else: v = refnc.variables[varname][:] if v.ndim == 2: - ref_vals = v + if varname in interpline: + ref_vals = v[0,:] + else: + ref_vals = v else: ref_vals = v[0,:] else: @@ -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",