From 5799ccf6d0f84a53009d44a5e1c462738c2967a2 Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Thu, 29 Mar 2018 15:12:18 -0600 Subject: [PATCH 1/2] Fixed bugs related to input failures with dictionaries. Modified routines that use qvapor to make copies so that scipy.io.netcdf works in the newer release. Added input tests. --- doc/source/new.rst | 2 + src/wrf/g_dewpoint.py | 8 +- src/wrf/g_latlon.py | 131 ++++++++++++++++++++-- src/wrf/g_rh.py | 8 +- src/wrf/g_slp.py | 4 +- src/wrf/interp.py | 22 ++-- src/wrf/latlonutils.py | 27 ++--- src/wrf/metadecorators.py | 8 +- src/wrf/util.py | 3 +- test/comp_utest.py | 6 +- test/ipynb/Doc_Examples.ipynb | 2 +- test/ipynb/Test_CTT.ipynb | 2 +- test/ipynb/WRF_python_demo.ipynb | 37 +++++-- test/test_inputs.py | 181 +++++++++++++++++++++++++++++++ test/utests.py | 2 +- 15 files changed, 381 insertions(+), 62 deletions(-) create mode 100644 test/test_inputs.py diff --git a/doc/source/new.rst b/doc/source/new.rst index c07e5a5..6103c43 100644 --- a/doc/source/new.rst +++ b/doc/source/new.rst @@ -21,8 +21,10 @@ v1.1.3 - Added 'th' alias for the theta product. - Fixed a crash issue related to updraft helicity when a dictionary is used as the input. +- Dictionary inputs now work correctly with xy_to_ll and ll_to_xy. - The cape_2d diagnostic can now work with a single column of data, just like cape_3d. + v1.1.2 ^^^^^^^^^^^^^^ diff --git a/src/wrf/g_dewpoint.py b/src/wrf/g_dewpoint.py index ad4a245..e9778d2 100755 --- a/src/wrf/g_dewpoint.py +++ b/src/wrf/g_dewpoint.py @@ -76,7 +76,9 @@ def get_dp(wrfin, timeidx=0, method="cat", squeeze=True, p = ncvars["P"] pb = ncvars["PB"] - qvapor = ncvars["QVAPOR"] + # Copy needed for the mmap nonsense of scipy.io.netcdf, which seems to + # break with every release + qvapor = ncvars["QVAPOR"].copy() # Algorithm requires hPa full_p = .01*(p + pb) @@ -152,7 +154,9 @@ def get_dp_2m(wrfin, timeidx=0, method="cat", squeeze=True, # Algorithm requires hPa psfc = .01*(ncvars["PSFC"]) - q2 = ncvars["Q2"] + # Copy needed for the mmap nonsense of scipy.io.netcdf, which seems to + # break with every release + q2 = ncvars["Q2"].copy() q2[q2 < 0] = 0 td = _td(psfc, q2) diff --git a/src/wrf/g_latlon.py b/src/wrf/g_latlon.py index d25f64a..44363c9 100755 --- a/src/wrf/g_latlon.py +++ b/src/wrf/g_latlon.py @@ -1,9 +1,16 @@ from __future__ import (absolute_import, division, print_function, unicode_literals) -from .util import extract_vars, get_id -from .latlonutils import (_lat_varname, _lon_varname, _ll_to_xy, _xy_to_ll) +from collections import OrderedDict + +import numpy as np +from xarray import DataArray + +from .util import extract_vars, get_id, get_iterable, is_mapping, to_np +from .py3compat import viewkeys +from .latlonutils import _lat_varname, _lon_varname, _ll_to_xy, _xy_to_ll from .metadecorators import set_latlon_metadata +from .config import xarray_enabled def get_lat(wrfin, timeidx=0, method="cat", squeeze=True, @@ -151,7 +158,101 @@ def get_lon(wrfin, timeidx=0, method="cat", squeeze=True, return lon_var[varname] -# TODO: Do we need the user to know about method, squeeze, cache for this? + +def _llxy_mapping(wrfin, x_or_lat, y_or_lon, func, timeidx, stagger, + squeeze, meta, as_int=None): + + keynames = [] + # This might not work once mapping iterators are implemented + numkeys = len(wrfin) + + key_iter = iter(viewkeys(wrfin)) + first_key = next(key_iter) + keynames.append(first_key) + + first_args = [wrfin[first_key], x_or_lat, y_or_lon, timeidx, squeeze, + meta, stagger] + if as_int is not None: + first_args.append(as_int) + + first_array = func(*first_args) + + # Create the output data numpy array based on the first array + outdims = [numkeys] + outdims += first_array.shape + outdata = np.empty(outdims, first_array.dtype) + outdata[0,:] = first_array[:] + + idx = 1 + while True: + try: + key = next(key_iter) + except StopIteration: + break + else: + keynames.append(key) + + args = [wrfin[first_key], x_or_lat, y_or_lon, timeidx, squeeze, + meta, stagger] + if as_int is not None: + args.append(as_int) + + result_array = func(*args) + + if outdata.shape[1:] != result_array.shape: + raise ValueError("data sequences must have the " + "same size for all dictionary keys") + outdata[idx,:] = to_np(result_array)[:] + idx += 1 + + if xarray_enabled() and meta: + outname = str(first_array.name) + # Note: assumes that all entries in dict have same coords + outcoords = OrderedDict(first_array.coords) + + # First find and store all the existing key coord names/values + # This is applicable only if there are nested dictionaries. + key_coordnames = [] + coord_vals = [] + existing_cnt = 0 + while True: + key_coord_name = "key_{}".format(existing_cnt) + + if key_coord_name not in first_array.dims: + break + + key_coordnames.append(key_coord_name) + coord_vals.append(to_np(first_array.coords[key_coord_name])) + + existing_cnt += 1 + + # Now add the key coord name and values for THIS dictionary. + # Put the new key_n name at the bottom, but the new values will + # be at the top to be associated with key_0 (left most). This + # effectively shifts the existing 'key_n' coordinate values to the + # right one dimension so *this* dicionary's key coordinate values + # are at 'key_0'. + key_coordnames.append(key_coord_name) + coord_vals.insert(0, keynames) + + # make it so that key_0 is leftmost + outdims = key_coordnames + list(first_array.dims[existing_cnt:]) + + + # Create the new 'key_n', value pairs + for coordname, coordval in zip(key_coordnames, coord_vals): + outcoords[coordname] = coordval + + + outattrs = OrderedDict(first_array.attrs) + + outarr = DataArray(outdata, name=outname, coords=outcoords, + dims=outdims, attrs=outattrs) + else: + outarr = outdata + + return outarr + @set_latlon_metadata(xy=True) def ll_to_xy(wrfin, latitude, longitude, timeidx=0, @@ -215,9 +316,15 @@ def ll_to_xy(wrfin, latitude, longitude, timeidx=0, be a :class:`numpy.ndarray` object with no metadata. """ + if is_mapping(wrfin): + return _llxy_mapping(wrfin, latitude, longitude, ll_to_xy, + timeidx, stagger, squeeze, meta, as_int) _key = get_id(wrfin) - return _ll_to_xy(latitude, longitude, wrfin, timeidx, stagger, "cat", + _wrfin = get_iterable(wrfin) + + return _ll_to_xy(latitude, longitude, _wrfin, timeidx, stagger, "cat", squeeze, None, _key, as_int, **{}) + @set_latlon_metadata(xy=True) @@ -319,10 +426,10 @@ def ll_to_xy_proj(latitude, longitude, meta=True, squeeze=True, as_int=True, return _ll_to_xy(latitude, longitude, None, 0, True, "cat", squeeze, None, None, as_int, **projparams) - - + + @set_latlon_metadata(xy=False) -def xy_to_ll(wrfin, x, y, timeidx=0, stagger=None, squeeze=True, meta=True): +def xy_to_ll(wrfin, x, y, timeidx=0, squeeze=True, meta=True, stagger=None): """Return the latitude and longitude for specified x,y coordinates. The *x* and *y* arguments can be a single value or a sequence of values. @@ -370,9 +477,6 @@ def xy_to_ll(wrfin, x, y, timeidx=0, stagger=None, squeeze=True, meta=True): - 'v': Use the same staggered grid as the v wind component, which has a staggered south_north (y) dimension. - as_int (:obj:`bool`): Set to True to return the x,y values as - :obj:`int`, otherwise they will be returned as :obj:`float`. - Returns: :class:`xarray.DataArray` or :class:`numpy.ndarray`: The latitude and longitude values whose leftmost dimension is 2 @@ -382,8 +486,13 @@ def xy_to_ll(wrfin, x, y, timeidx=0, stagger=None, squeeze=True, meta=True): be a :class:`numpy.ndarray` object with no metadata. """ + if is_mapping(wrfin): + return _llxy_mapping(wrfin, x, y, xy_to_ll, + timeidx, stagger, squeeze, meta) + _key = get_id(wrfin) - return _xy_to_ll(x, y, wrfin, timeidx, stagger, "cat", True, None, + _wrfin = get_iterable(wrfin) + return _xy_to_ll(x, y, _wrfin, timeidx, stagger, "cat", True, None, _key, **{}) diff --git a/src/wrf/g_rh.py b/src/wrf/g_rh.py index 5af8ae5..6e4bb6b 100755 --- a/src/wrf/g_rh.py +++ b/src/wrf/g_rh.py @@ -71,7 +71,9 @@ def get_rh(wrfin, timeidx=0, method="cat", squeeze=True, cache=None, t = ncvars["T"] p = ncvars["P"] pb = ncvars["PB"] - qvapor = ncvars["QVAPOR"] + # Copy needed for the mmap nonsense of scipy.io.netcdf, which seems to + # break with every release + qvapor = ncvars["QVAPOR"].copy() full_t = t + Constants.T_BASE full_p = p + pb @@ -144,7 +146,9 @@ def get_rh_2m(wrfin, timeidx=0, method="cat", squeeze=True, cache=None, meta=False, _key=_key) t2 = ncvars["T2"] psfc = ncvars["PSFC"] - q2 = ncvars["Q2"] + # Copy needed for the mmap nonsense of scipy.io.netcdf, which seems to + # break with every release + q2 = ncvars["Q2"].copy() q2[q2 < 0] = 0 rh = _rh(q2, psfc, t2) diff --git a/src/wrf/g_slp.py b/src/wrf/g_slp.py index a1f5a96..77bac48 100755 --- a/src/wrf/g_slp.py +++ b/src/wrf/g_slp.py @@ -83,7 +83,9 @@ def get_slp(wrfin, timeidx=0, method="cat", squeeze=True, t = ncvars["T"] p = ncvars["P"] pb = ncvars["PB"] - qvapor = ncvars["QVAPOR"] + # Copy needed for the mmap nonsense of scipy.io.netcdf, which seems to + # break with every release + qvapor = ncvars["QVAPOR"].copy() ph = ncvars["PH"] phb = ncvars["PHB"] diff --git a/src/wrf/interp.py b/src/wrf/interp.py index 1892947..34d7e6f 100755 --- a/src/wrf/interp.py +++ b/src/wrf/interp.py @@ -8,7 +8,7 @@ from .extension import (_interpz3d, _vertcross, _interpline, _smooth2d, _monotonic, _vintrp) from .metadecorators import set_interp_metadata -from .util import extract_vars, is_staggered, get_id, to_np +from .util import extract_vars, is_staggered, get_id, to_np, get_iterable from .py3compat import py3range from .interputils import get_xy, get_xy_z_params, to_xy_coords from .constants import Constants, default_fill, ConversionFactors @@ -548,6 +548,8 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False, """ _key = get_id(wrfin) + _wrfin = get_iterable(wrfin) + # Remove case sensitivity field_type = field_type.lower() if field_type is not None else "none" vert_coord = vert_coord.lower() if vert_coord is not None else "none" @@ -583,7 +585,7 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False, interp_levels = np.asarray(interp_levels, np.float64) # TODO: Check if field is staggered - if is_staggered(wrfin, field): + if is_staggered(_wrfin, field): raise RuntimeError("Please unstagger field in the vertical") # Check for valid coord @@ -606,23 +608,23 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False, # Extract vriables #timeidx = -1 # Should this be an argument? - ncvars = extract_vars(wrfin, timeidx, ("PSFC", "QVAPOR", "F"), + ncvars = extract_vars(_wrfin, timeidx, ("PSFC", "QVAPOR", "F"), method, squeeze, cache, meta=False, _key=_key) sfp = ncvars["PSFC"] * ConversionFactors.PA_TO_HPA qv = ncvars["QVAPOR"] coriolis = ncvars["F"] - terht = get_terrain(wrfin, timeidx, units="m", + terht = get_terrain(_wrfin, timeidx, units="m", method=method, squeeze=squeeze, cache=cache, meta=False, _key=_key) - tk = get_temp(wrfin, timeidx, units="k", + tk = get_temp(_wrfin, timeidx, units="k", method=method, squeeze=squeeze, cache=cache, meta=False, _key=_key) - p = get_pressure(wrfin, timeidx, units="pa", + p = get_pressure(_wrfin, timeidx, units="pa", method=method, squeeze=squeeze, cache=cache, meta=False, _key=_key) - ght = get_height(wrfin, timeidx, msl=True, units="m", + ght = get_height(_wrfin, timeidx, msl=True, units="m", method=method, squeeze=squeeze, cache=cache, meta=False, _key=_key) @@ -639,7 +641,7 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False, vcord_array = np.exp(-ght/sclht) elif vert_coord == "ght_agl": - ht_agl = get_height(wrfin, timeidx, msl=False, units="m", + ht_agl = get_height(_wrfin, timeidx, msl=False, units="m", method=method, squeeze=squeeze, cache=cache, meta=False, _key=_key) @@ -647,7 +649,7 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False, vcord_array = np.exp(-ht_agl/sclht) elif vert_coord in ("theta", "th"): - t = get_theta(wrfin, timeidx, units="k", + t = get_theta(_wrfin, timeidx, units="k", method=method, squeeze=squeeze, cache=cache, meta=False, _key=_key) @@ -671,7 +673,7 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False, idir = 1 delta = 0.01 - eth = get_eth(wrfin, timeidx, method=method, squeeze=squeeze, + eth = get_eth(_wrfin, timeidx, method=method, squeeze=squeeze, cache=cache, meta=False, _key=_key) p_hpa = p * ConversionFactors.PA_TO_HPA diff --git a/src/wrf/latlonutils.py b/src/wrf/latlonutils.py index ccd5221..00940ef 100644 --- a/src/wrf/latlonutils.py +++ b/src/wrf/latlonutils.py @@ -131,8 +131,9 @@ def _get_proj_params(wrfin, timeidx, stagger, method, squeeze, cache, _key): """ - if timeidx < 0: - raise ValueError("'timeidx' must be greater than 0") + if timeidx is not None: + if timeidx < 0: + raise ValueError("'timeidx' must be greater than 0") attrs = extract_global_attrs(wrfin, attrs=("MAP_PROJ", "TRUELAT1", "TRUELAT2", "STAND_LON", @@ -383,12 +384,10 @@ def _ll_to_xy(latitude, longitude, wrfin=None, timeidx=0, if ref_lat.size == 1: outdim = [2, lats.size] - #outdim = [lats.size, 2] extra_dims = [outdim[1]] else: # Moving domain will have moving ref_lats/ref_lons outdim = [2, ref_lat.size, lats.size] - #outdim = [lats.size, ref_lat.size, 2] extra_dims = outdim[1:] result = np.empty(outdim, np.float64) @@ -402,11 +401,11 @@ def _ll_to_xy(latitude, longitude, wrfin=None, timeidx=0, ref_lat_val = ref_lat[0] ref_lon_val = ref_lon[0] else: - ref_lat_val = ref_lat[left_idxs[-1]] - ref_lon_val = ref_lon[left_idxs[-1]] + ref_lat_val = ref_lat[left_idxs[-2]] + ref_lon_val = ref_lon[left_idxs[-2]] - lat = lats[left_idxs[0]] - lon = lons[left_idxs[0]] + lat = lats[left_idxs[-1]] + lon = lons[left_idxs[-1]] xy = _lltoxy(map_proj, truelat1, truelat2, stdlon, ref_lat_val, ref_lon_val, pole_lat, pole_lon, @@ -548,14 +547,10 @@ def _xy_to_ll(x, y, wrfin=None, timeidx=0, stagger=None, raise ValueError("'x' and 'y' must be the same length") if ref_lat.size == 1: - #outdim = [x_arr.size, 2] - #extra_dims = [outdim[0]] outdim = [2, x_arr.size] extra_dims = [outdim[1]] else: # Moving domain will have moving ref_lats/ref_lons - #outdim = [x_arr.size, ref_lat.size, 2] - #extra_dims = outdim[0:2] outdim = [2, ref_lat.size, x_arr.size] extra_dims = outdim[1:] @@ -570,11 +565,11 @@ def _xy_to_ll(x, y, wrfin=None, timeidx=0, stagger=None, ref_lat_val = ref_lat[0] ref_lon_val = ref_lon[0] else: - ref_lat_val = ref_lat[left_idxs[-1]] - ref_lon_val = ref_lon[left_idxs[-1]] + ref_lat_val = ref_lat[left_idxs[-2]] + ref_lon_val = ref_lon[left_idxs[-2]] - x_val = x_arr[left_idxs[0]] - y_val = y_arr[left_idxs[0]] + x_val = x_arr[left_idxs[-1]] + y_val = y_arr[left_idxs[-1]] ll = _xytoll(map_proj, truelat1, truelat2, stdlon, ref_lat_val, ref_lon_val, pole_lat, pole_lon, known_x, known_y, diff --git a/src/wrf/metadecorators.py b/src/wrf/metadecorators.py index 40e8e49..771164b 100644 --- a/src/wrf/metadecorators.py +++ b/src/wrf/metadecorators.py @@ -10,7 +10,7 @@ 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) + from_var, iter_left_indexes, is_mapping) from .coordpair import CoordPair from .py3compat import viewkeys, viewitems, py3range, ucode from .interputils import get_xy_z_params, get_xy, to_xy_coords @@ -586,8 +586,10 @@ def set_latlon_metadata(xy=False): @wrapt.decorator def func_wrapper(wrapped, instance, args, kwargs): - do_meta = from_args(wrapped, ("meta",), *args, **kwargs)["meta"] - + argvars = from_args(wrapped, ("wrfin", "meta"), *args, **kwargs) + # If it's a mapping, then this is handled as a special case in g_latlon + do_meta = (not is_mapping(argvars["wrfin"]) and argvars["meta"]) + if do_meta is None: do_meta = True diff --git a/src/wrf/util.py b/src/wrf/util.py index 7bf9927..2c0fd9e 100644 --- a/src/wrf/util.py +++ b/src/wrf/util.py @@ -3073,7 +3073,8 @@ def get_id(obj, prefix=''): # For sequences, the hashing string will be the list ID and the # path for the first file in the sequence if not is_mapping(obj): - _next = next(iter(obj)) + _obj = get_iterable(obj) + _next = next(iter(_obj)) return get_id(_next, prefix + str(id(obj))) # For each key in the mapping, recursively call get_id until diff --git a/test/comp_utest.py b/test/comp_utest.py index 864bb7a..de46a5f 100644 --- a/test/comp_utest.py +++ b/test/comp_utest.py @@ -642,15 +642,11 @@ def test_cape2d_1d(wrfnc): ter_follow = 1 result = cape_2d(p_hpa, tkel, qv, z, ter, psfc_hpa, ter_follow) - - print ("RESULT", result) ref = getvar(wrfnc, "cape_2d") ref = ref[(slice(None),) + col_idxs[1:]] - print ("REF", ref) - nt.assert_allclose(to_np(result), to_np(ref)) return func @@ -658,7 +654,7 @@ def test_cape2d_1d(wrfnc): 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()-1) + omp_set_num_threads(omp_get_num_procs()//2) omp_set_schedule(OMP_SCHED_STATIC, 0) omp_set_dynamic(False) diff --git a/test/ipynb/Doc_Examples.ipynb b/test/ipynb/Doc_Examples.ipynb index 4b8d794..81ba5bd 100644 --- a/test/ipynb/Doc_Examples.ipynb +++ b/test/ipynb/Doc_Examples.ipynb @@ -1275,7 +1275,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", - "version": "2.7.13" + "version": "2.7.14" } }, "nbformat": 4, diff --git a/test/ipynb/Test_CTT.ipynb b/test/ipynb/Test_CTT.ipynb index 533ea13..c0c1a36 100644 --- a/test/ipynb/Test_CTT.ipynb +++ b/test/ipynb/Test_CTT.ipynb @@ -206,7 +206,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", - "version": "2.7.13" + "version": "2.7.14" } }, "nbformat": 4, diff --git a/test/ipynb/WRF_python_demo.ipynb b/test/ipynb/WRF_python_demo.ipynb index 977aaed..8f073d0 100644 --- a/test/ipynb/WRF_python_demo.ipynb +++ b/test/ipynb/WRF_python_demo.ipynb @@ -251,7 +251,7 @@ " return self\n", " \n", " def next(self):\n", - " if self._i > self._total:\n", + " if self._i >= self._total:\n", " raise StopIteration\n", " else:\n", " val = self.ncfile[self._i]\n", @@ -653,7 +653,7 @@ "metadata": {}, "outputs": [], "source": [ - "from wrf.latlon import xy_to_ll, ll_to_xy \n", + "from wrf import xy_to_ll, ll_to_xy \n", "\n", "a = xy_to_ll(ncfile, 400, 200)\n", "a1 = ll_to_xy(ncfile, a[0], a[1])\n", @@ -1134,7 +1134,7 @@ "metadata": {}, "outputs": [], "source": [ - "from wrf.latlon import xy_to_ll, ll_to_xy \n", + "from wrf import xy_to_ll, ll_to_xy \n", "\n", "a = xy_to_ll(ncfiles, 400, 200)\n", "a = xy_to_ll(ncfiles, [400,105], [200,205])\n", @@ -1196,6 +1196,27 @@ "print (\"\\n\")\n" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, { "cell_type": "code", "execution_count": null, @@ -1206,21 +1227,21 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 2", + "display_name": "Python 3", "language": "python", - "name": "python2" + "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", - "version": 2 + "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython2", - "version": "2.7.13" + "pygments_lexer": "ipython3", + "version": "3.6.4" } }, "nbformat": 4, diff --git a/test/test_inputs.py b/test/test_inputs.py new file mode 100644 index 0000000..1fa986e --- /dev/null +++ b/test/test_inputs.py @@ -0,0 +1,181 @@ +import unittest as ut +import numpy.testing as nt +import numpy as np +import numpy.ma as ma +import os +import sys +import subprocess + +from wrf import (getvar, interplevel, interpline, vertcross, vinterp, + disable_xarray, xarray_enabled, to_np, + xy_to_ll, ll_to_xy, xy_to_ll_proj, ll_to_xy_proj, + extract_global_attrs, viewitems, CoordPair, ll_points) +from wrf.util import is_multi_file + +IN_DIR = "/Users/ladwig/Documents/wrf_files/wrf_vortex_multi" +TEST_FILES = [os.path.join(IN_DIR, "wrfout_d02_2005-08-28_00:00:00"), + os.path.join(IN_DIR, "wrfout_d02_2005-08-28_12:00:00"), + os.path.join(IN_DIR, "wrfout_d02_2005-08-29_00:00:00")] + +def wrfin_gen(wrf_in): + for x in wrf_in: + yield x + + +class wrf_in_iter_class(object): + def __init__(self, wrf_in): + self._wrf_in = wrf_in + self._total = len(wrf_in) + self._i = 0 + + def __iter__(self): + return self + + def next(self): + if self._i >= self._total: + raise StopIteration + else: + result = self._wrf_in[self._i] + self._i += 1 + return result + + # Python 3 + def __next__(self): + return self.next() + + +def make_test(varname, wrf_in): + def test(self): + x = getvar(wrf_in, varname, 0) + xa = getvar(wrf_in, varname, None) + + return test + +def make_interp_test(varname, wrf_in): + def test(self): + + # Only testing vinterp since other interpolation just use variables + if (varname == "vinterp"): + for timeidx in (0, None): + eth = getvar(wrf_in, "eth", timeidx=timeidx) + interp_levels = [850,500,5] + field = vinterp(wrf_in, + field=eth, + vert_coord="pressure", + interp_levels=interp_levels, + extrapolate=True, + field_type="theta-e", + timeidx=timeidx, + log_p=True) + else: + pass + + + return test + + +def make_latlon_test(testid, wrf_in): + def test(self): + if testid == "xy_out": + # Lats/Lons taken from NCL script, just hard-coding for now + lats = [-55, -60, -65] + lons = [25, 30, 35] + xy = ll_to_xy(wrf_in, lats, lons, timeidx=0) + xy = ll_to_xy(wrf_in, lats, lons, timeidx=None) + else: + i_s = np.asarray([10, 100, 150], int) + j_s = np.asarray([10, 100, 150], int) + ll = xy_to_ll(wrf_in, i_s, j_s, timeidx=0) + ll = xy_to_ll(wrf_in, i_s, j_s, timeidx=None) + + return test + + +class WRFVarsTest(ut.TestCase): + longMessage = True + +class WRFInterpTest(ut.TestCase): + longMessage = True + +class WRFLatLonTest(ut.TestCase): + longMessage = True + +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()-1) + 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", + "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"] + interp_methods = ["interplevel", "vertcross", "interpline", "vinterp"] + latlon_tests = ["xy_out", "ll_out"] + + + for nc_lib in ("netcdf4", "pynio", "scipy"): + if nc_lib == "netcdf4": + try: + from netCDF4 import Dataset + except ImportError: + print ("netcdf4-python not installed") + continue + else: + test_in = [Dataset(x) for x in TEST_FILES] + elif nc_lib == "pynio": + try: + from Nio import open_file + except ImportError: + print ("PyNIO not installed") + continue + else: + test_in = [open_file(x +".nc", "r") for x in TEST_FILES] + elif nc_lib == "scipy": + try: + from scipy.io.netcdf import netcdf_file + except ImportError: + print ("scipy.io.netcdf not installed") + else: + test_in = [netcdf_file(x, mmap=False) for x in TEST_FILES] + + input0 = test_in[0] + input1 = test_in + input2 = tuple(input1) + input3 = wrfin_gen(test_in) + input4 = wrf_in_iter_class(test_in) + input5 = {"A" : input1, + "B" : input2} + input6 = {"A" : {"AA" : input1}, + "B" : {"BB" : input2}} + + for i,input in enumerate((input0, input1, input2, + input3, input4, input5, input6)): + for var in wrf_vars: + if var in ignore_vars: + continue + test_func1 = make_test(var, input) + + setattr(WRFVarsTest, "test_{0}_input{1}_{2}".format(nc_lib, + i, var), + test_func1) + + + for method in interp_methods: + test_interp_func1 = make_interp_test(method, input) + setattr(WRFInterpTest, + "test_{0}_input{1}_{2}".format(nc_lib, + i, method), + test_interp_func1) + + for testid in latlon_tests: + test_ll_func = make_latlon_test(testid, input) + test_name = "test_{0}_input{1}_{2}".format(nc_lib, i, testid) + setattr(WRFLatLonTest, test_name, test_ll_func) + + ut.main() + + \ No newline at end of file diff --git a/test/utests.py b/test/utests.py index 4962ef5..ed1d452 100644 --- a/test/utests.py +++ b/test/utests.py @@ -682,7 +682,7 @@ 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()-1) + omp_set_num_threads(omp_get_num_procs()//2) omp_set_schedule(OMP_SCHED_STATIC, 0) omp_set_dynamic(False) From 2c35edd7d7e42d5388615ae2b876345ae8b6502e Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Thu, 29 Mar 2018 16:00:28 -0600 Subject: [PATCH 2/2] Added ctt test --- test/ctt_test.py | 50 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 50 insertions(+) create mode 100644 test/ctt_test.py diff --git a/test/ctt_test.py b/test/ctt_test.py new file mode 100644 index 0000000..6e97909 --- /dev/null +++ b/test/ctt_test.py @@ -0,0 +1,50 @@ +from netCDF4 import Dataset +import matplotlib.pyplot as plt +from matplotlib.cm import get_cmap +import cartopy.crs as crs +from cartopy.feature import NaturalEarthFeature + +from wrf import to_np, getvar, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords + +# Open the NetCDF file +ncfile = Dataset("/Users/ladwig/Documents/wrf_files/problem_files/cfrac_bug/wrfout_d02_1987-10-01_00:00:00") + +# Get the sea level pressure +ctt = getvar(ncfile, "ctt") + +# Get the latitude and longitude points +lats, lons = latlon_coords(ctt) + +# Get the cartopy mapping object +cart_proj = get_cartopy(ctt) + +# Create a figure +fig = plt.figure(figsize=(12,9)) +# Set the GeoAxes to the projection used by WRF +ax = plt.axes(projection=cart_proj) + +# Download and add the states and coastlines +states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none', + name='admin_1_states_provinces_shp') +ax.add_feature(states, linewidth=.5) +ax.coastlines('50m', linewidth=0.8) + +# Make the contour outlines and filled contours for the smoothed sea level pressure. +plt.contour(to_np(lons), to_np(lats), to_np(ctt), 10, colors="black", + transform=crs.PlateCarree()) +plt.contourf(to_np(lons), to_np(lats), to_np(ctt), 10, transform=crs.PlateCarree(), + cmap=get_cmap("jet")) + +# Add a color bar +plt.colorbar(ax=ax, shrink=.62) + +# Set the map limits. Not really necessary, but used for demonstration. +ax.set_xlim(cartopy_xlim(ctt)) +ax.set_ylim(cartopy_ylim(ctt)) + +# Add the gridlines +ax.gridlines(color="black", linestyle="dotted") + +plt.title("Cloud Top Temperature") + +plt.show()