diff --git a/src/wrf/interputils.py b/src/wrf/interputils.py index 04807d1..8e04c86 100644 --- a/src/wrf/interputils.py +++ b/src/wrf/interputils.py @@ -156,13 +156,6 @@ def _calc_xy(xdim, ydim, pivot_point=None, angle=None, if x1 >= xdim or y1 >= ydim: raise ValueError("end_point {} is outside of domain " "with shape {}".format(end_point, (xdim, ydim))) - - # From the original NCL code, but the error above will be thrown - # instead. - if ( x1 > xdim-1 ): - x1 = xdim - 1 - if ( y1 > ydim-1): - y1 = ydim - 1 else: raise ValueError("invalid start/end or pivot/angle arguments") diff --git a/src/wrf/util.py b/src/wrf/util.py index be174b4..b9431be 100644 --- a/src/wrf/util.py +++ b/src/wrf/util.py @@ -3880,13 +3880,21 @@ def pairs_to_latlon(pairs): return lats, lons - +def is_latlon_pair(pair): + """Return True if the :class:`wrf.CoordPair` is a lat/lon pair + + Args: + + pair (:class:`wrf.CoordPair`): A single :class:`wrf.CoordPair` object. + + Returns: + + :obj:`bool`: True if the pair is a lat/lon pair. + + """ + return (pair.lat is not None and pair.lon is not None) - - - - diff --git a/test/ncl_get_var.ncl b/test/ncl_get_var.ncl index 86d85fb..bd4d3a6 100644 --- a/test/ncl_get_var.ncl +++ b/test/ncl_get_var.ncl @@ -5,9 +5,9 @@ ;system("printenv") if (.not. isvar("dir")) then - dir = "/Users/ladwig/Documents/wrf_files/wrf_vortex_multi" + dir = "/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/moving_nest" ;in_file = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2010-06-13_21:00:00.nc" - ;in_file = "/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/wrfout_d02_2005-08-28_00:00:00.nc" + ;in_file = "/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/moving_nest/wrfout_d02_2005-08-28_00:00:00.nc" end if if (.not. isvar("pattern")) then @@ -63,10 +63,10 @@ z := wrf_user_getvar(input_file, "z", time) ; grid point height p := wrf_user_getvar(input_file, "pressure", time) ; total pressure + dimsz = dimsizes(z) - pivot = (/ dimsz(2)/2, dimsz(1)/2 /) ; pivot point is center of domain + pivot = (/ dimsz(3)/2, dimsz(2)/2 /) ; pivot point is center of domain - time = 0 ht_cross = wrf_user_intrp3d(z,p,"v",pivot,90.0,False) p_cross = wrf_user_intrp3d(p,z,"v",pivot,90.0,False) @@ -75,7 +75,10 @@ fout->ht_cross = ht_cross fout->p_cross = p_cross -; + + + time = 0 + ; For the new cross section routine xopt = True xopt@use_pivot = True @@ -85,7 +88,6 @@ xopt@linecoords = True ht_vertcross1 = wrf_user_vertcross(z, p, pivot, xopt) - printVarSummary(ht_vertcross1) fout->ht_vertcross1 = ht_vertcross1 @@ -99,8 +101,8 @@ xopt@linecoords = True ht_vertcross2 = wrf_user_vertcross(z, p, pivot, xopt) - ht_vertcross2!0 = "vertical2" - ht_vertcross2!1 = "cross_line_idx2" + ht_vertcross2!1 = "vertical2" + ht_vertcross2!2 = "cross_line_idx2" fout->ht_vertcross2 = ht_vertcross2 @@ -113,10 +115,10 @@ 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 /) @@ -125,20 +127,39 @@ xopt@use_pivot = False xopt@latlon = True xopt@file_handle = input_file - xopt@timeidx = -1 + xopt@timeidx = 0 xopt@linecoords = True xopt@autolevels = 1000 ht_vertcross3 = wrf_user_vertcross(z, p, start_end, xopt) - printVarSummary(ht_vertcross3) - ht_vertcross3!0 = "Time" ht_vertcross3!1 = "vertical3" ht_vertcross3!2 = "cross_line_idx3" fout->ht_vertcross3 = ht_vertcross3 + ; Test the moving nest with lat/lon over time + + times = wrf_user_getvar(input_file, "times", -1) + ntimes = dimsizes(times) + + do i=0,ntimes-1 + xopt@timeidx = i + name = sprinti("ht_vertcross_t%i", i) + p_var := p(i,:,:,:) + z_var := z(i,:,:,:) + + ht_vertcross := wrf_user_vertcross(z_var, p_var, start_end, xopt) + + dim0name = sprinti("vertical_t%i",i) + dim1name = sprinti("cross_line_idx_t%i",i) + ht_vertcross!0 = dim0name + ht_vertcross!1 = dim1name + + fout->$name$ = ht_vertcross + end do + ;;;;;;;;;;;;;;;;;;;;;;;; 3D horizontal interpolation time = -1 @@ -202,7 +223,52 @@ t2_line = wrf_user_intrp2d(t2, pivot, 90.0, False) - fout->t2_line = (/t2_line/) + fout->t2_line = t2_line + + ; For the new interplevel routine + xopt := True + xopt@use_pivot = True + xopt@angle = 90.0 + xopt@latlon = False + xopt@file_handle = input_file + xopt@timeidx = 0 + xopt@linecoords = True + + t2_line2 = wrf_user_interpline(t2, pivot, xopt) + + 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 + xopt@latlon = True + xopt@file_handle = input_file + xopt@timeidx = 0 + xopt@linecoords = True + + t2_line3 = wrf_user_interpline(t2, start_end, xopt) + + times = wrf_user_getvar(input_file, "times", -1) + ntimes = dimsizes(times) + + do i=0,ntimes-1 + xopt@timeidx = i + name = sprinti("t2_line_t%i", i) + dim0name = sprinti("lineidx_t%i",i) + var := t2(i,:,:) + t2_line := wrf_user_interpline(var, start_end, xopt) + t2_line!0 = dim0name + fout->$name$ = t2_line + end do ;;;;;;;;;;;;;;;;;;;;;;; 3D interpolation to new vertical coordinates diff --git a/test/utests.py b/test/utests.py index dc8eeb7..5310a2e 100644 --- a/test/utests.py +++ b/test/utests.py @@ -15,7 +15,7 @@ from wrf.util import is_multi_file NCL_EXE = "/Users/ladwig/miniconda2/envs/ncl_build/bin/ncl" NCARG_ROOT = "/Users/ladwig/miniconda2/envs/ncl_build" #TEST_FILE = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2010-06-13_21:00:00" -DIR = "/Users/ladwig/Documents/wrf_files/wrf_vortex_multi" +DIR = "/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/moving_nest" PATTERN = "wrfout_d02_*" REF_NC_FILE = "/tmp/wrftest.nc" @@ -74,7 +74,7 @@ def make_test(varname, dir, pattern, referent, multi=False, pynio=False): if not pynio: f = NetCDF(fname) try: - f.set_auto_mask(False) + f.set_always_mask(False) except: pass wrfin.append(f) @@ -190,7 +190,8 @@ def _get_refvals(referent, varname, multi): refnc = NetCDF(referent) try: - refnc.set_auto_mask(False) + pass + #refnc.set_auto_mask(False) except: pass @@ -236,7 +237,7 @@ def make_interp_test(varname, dir, pattern, referent, multi=False, if not pynio: f = NetCDF(fname) try: - f.set_auto_mask(False) + f.set_always_mask(False) except: pass wrfin.append(f) @@ -353,33 +354,40 @@ def make_interp_test(varname, dir, pattern, referent, multi=False, 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., autolevels=p_autolevels) + ht_cross = vertcross(hts, p, pivot_point=pivot_point, angle=90., autolevels=p_autolevels) nt.assert_allclose(to_np(ht_cross), ref_ht_cross, rtol=.01) - # Test the manual projection method with lat/lon lats = hts.coords["XLAT"] lons = hts.coords["XLONG"] - ll_point = ll_points(lats, lons) - pivot = CoordPair(lat=lats[int(lats.shape[-2]/2), + + # Test the manual projection method with lat/lon + # Only do this for the non-multi case, since the domain + # might be moving + if not multi: + + ll_point = ll_points(lats, lons) + + pivot = CoordPair(lat=lats[int(lats.shape[-2]/2), int(lats.shape[-1]/2)], - lon=lons[int(lons.shape[-2]/2), + lon=lons[int(lons.shape[-2]/2), int(lons.shape[-1]/2)]) - v1 = vertcross(hts,p,wrfin,pivot_point=pivot_point, + v1 = vertcross(hts,p,wrfin=wrfin,pivot_point=pivot_point, angle=90.0) - v2 = vertcross(hts,p,projection=hts.attrs["projection"], + v2 = vertcross(hts,p,projection=hts.attrs["projection"], ll_point=ll_point, pivot_point=pivot_point, angle=90.) - nt.assert_allclose(to_np(v1), to_np(v2), rtol=.01) + nt.assert_allclose(to_np(v1), to_np(v2), rtol=.01) # Test opposite @@ -390,7 +398,7 @@ def make_interp_test(varname, dir, pattern, referent, multi=False, rtol=.01) # Test point to point start_point = CoordPair(0, hts.shape[-2]/2) - end_point = CoordPair(-1,hts.shape[-2]/2) + end_point = CoordPair(-1, hts.shape[-2]/2) p_cross2 = vertcross(p,hts,start_point=start_point, @@ -418,17 +426,19 @@ def make_interp_test(varname, dir, pattern, referent, multi=False, nt.assert_allclose(to_np(ht_cross), to_np(ref_ht_vertcross2), atol=.01) + idxs = (0, slice(None)) if multi else (slice(None),) - 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_lat = np.amin(lats[idxs]) + .25*(np.amax(lats[idxs]) - np.amin(lats[idxs])) + end_lat = np.amin(lats[idxs]) + .65*(np.amax(lats[idxs]) - np.amin(lats[idxs])) - 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_lon = np.amin(lons[idxs]) + .25*(np.amax(lons[idxs]) - np.amin(lons[idxs])) + end_lon = np.amin(lons[idxs]) + .65*(np.amax(lons[idxs]) - np.amin(lons[idxs])) 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 + ll_point = ll_points(lats[idxs], lons[idxs]) + ht_cross = vertcross(hts, p, start_point=start_point, end_point=end_point, @@ -442,6 +452,27 @@ def make_interp_test(varname, dir, pattern, referent, multi=False, to_np(ref_ht_vertcross3), rtol=.01) + if multi: + ntimes = hts.shape[0] + + for t in range(ntimes): + hts = getvar(wrfin, "z", timeidx=t) + p = getvar(wrfin, "pressure", timeidx=t) + + ht_cross = vertcross(hts, p, + start_point=start_point, + end_point=end_point, + wrfin=wrfin, + timeidx=t, + latlon=True, + autolevels=1000) + + refname = "ht_vertcross_t{}".format(t) + ref_ht_vertcross = _get_refvals(referent, refname, False) + + nt.assert_allclose(to_np(ht_cross), + to_np(ref_ht_vertcross),rtol=.02) + elif (varname == "interpline"): @@ -700,7 +731,7 @@ def make_latlon_test(testid, dir, pattern, referent, single, refnc = NetCDF(referent) try: - refnc.set_auto_mask(False) + refnc.set_always_mask(False) except: pass @@ -813,9 +844,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",