Browse Source

Merge branch 'feature/contrib' into develop

lon0
Bill Ladwig 6 years ago
parent
commit
b40d1d6e2d
  1. 63
      doc/source/contrib.rst
  2. 10
      doc/source/index.rst
  3. 13
      src/wrf/__init__.py
  4. 1
      src/wrf/api.py
  5. 5
      src/wrf/cache.py
  6. 13
      src/wrf/computation.py
  7. 4
      src/wrf/config.py
  8. 11
      src/wrf/constants.py
  9. 31
      src/wrf/coordpair.py
  10. 40
      src/wrf/decorators.py
  11. 2
      src/wrf/destag.py
  12. 41
      src/wrf/extension.py
  13. 5
      src/wrf/g_cape.py
  14. 20
      src/wrf/g_cloudfrac.py
  15. 1
      src/wrf/g_ctt.py
  16. 2
      src/wrf/g_dbz.py
  17. 3
      src/wrf/g_dewpoint.py
  18. 5
      src/wrf/g_geoht.py
  19. 7
      src/wrf/g_helicity.py
  20. 28
      src/wrf/g_latlon.py
  21. 1
      src/wrf/g_omega.py
  22. 2
      src/wrf/g_precip.py
  23. 5
      src/wrf/g_pressure.py
  24. 5
      src/wrf/g_pw.py
  25. 2
      src/wrf/g_rh.py
  26. 2
      src/wrf/g_slp.py
  27. 3
      src/wrf/g_temp.py
  28. 3
      src/wrf/g_terrain.py
  29. 1
      src/wrf/g_times.py
  30. 15
      src/wrf/g_uvmet.py
  31. 1
      src/wrf/g_vorticity.py
  32. 7
      src/wrf/g_wind.py
  33. 5
      src/wrf/geobnds.py
  34. 65
      src/wrf/interp.py
  35. 8
      src/wrf/interputils.py
  36. 23
      src/wrf/latlonutils.py
  37. 99
      src/wrf/metadecorators.py
  38. 74
      src/wrf/projection.py
  39. 1
      src/wrf/projutils.py
  40. 4
      src/wrf/py3compat.py
  41. 9
      src/wrf/routines.py
  42. 15
      src/wrf/specialdec.py
  43. 14
      src/wrf/units.py
  44. 137
      src/wrf/util.py
  45. 1
      src/wrf/version.py

63
doc/source/contrib.rst

@ -0,0 +1,63 @@
.. _contrib_guide:
Contributor Guide
=================================
.. note::
This contributor guide is written for wrf-python v1.3.x. In the
not-too-distant future, wrf-python will undergo a significant refactoring
to remove the wrapt decorators (which don't serialize for dask), but the
concepts will remain the same as described below.
Ways to Contribute
-----------------------------
Users are encouraged to contribute various ways. This includes:
- Submitting a bug report
- Submitting bug fixes
- Submitting new Fortran computational routines
- Submitting new Python computational routines
- Submitting fully wrapped computational routines
Getting the source code
------------------------------
The source code is available on GitHub:
https://github.com/NCAR/wrf-python
To checkout the code::
git clone https://github.com/NCAR/wrf-python
Git Flow
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
This project follows the GitFlow Workflow, which you can read about here if it
is new to you:
https://leanpub.com/git-flow/read
When you first clone the repository, by default you will be on the 'develop'
branch, which is what you should use for your development.
Pull Requests
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
In order to submit changes, you must use GitHub to issue a pull request.
Overview of WRF-Python Internals
----------------------------------
WRF-Python is a collection of diagnostic and interpolation routines for WRF-ARW
data. The API consists of a handful of functions

10
doc/source/index.rst

@ -13,6 +13,14 @@
university corporation for atmospheric research, university corporation for atmospheric research,
pynio, pyngl, interpolation pynio, pyngl, interpolation
.. .. image:: _static/images/nsf.png
.. :scale: 100%
.. :align: right
.. |
.. |
wrf-python wrf-python
=========== ===========
@ -57,6 +65,7 @@ Indices and tables
* :ref:`modindex` * :ref:`modindex`
* :ref:`search` * :ref:`search`
-------------------- --------------------
*The National Center for Atmospheric Research is sponsored by the National *The National Center for Atmospheric Research is sponsored by the National
@ -64,4 +73,3 @@ Science Foundation. Any opinions, findings and conclusions or recommendations
expressed in this material do not necessarily reflect the views of the expressed in this material do not necessarily reflect the views of the
National Science Foundation.* National Science Foundation.*

13
src/wrf/__init__.py

@ -2,19 +2,24 @@ from __future__ import (absolute_import, division, print_function)
import os import os
import pkg_resources import pkg_resources
# For gfortran+msvc combination, extra shared libraries may exist (stored by numpy.distutils)
if os.name == "nt":
try: try:
from . import api
from .api import *
except ImportError:
# For gfortran+msvc combination, extra shared libraries may exist
# (stored by numpy.distutils)
if os.name == "nt":
req = pkg_resources.Requirement.parse("wrf-python") req = pkg_resources.Requirement.parse("wrf-python")
extra_dll_dir = pkg_resources.resource_filename(req, extra_dll_dir = pkg_resources.resource_filename(req,
"wrf-python/.libs") "wrf-python/.libs")
if os.path.isdir(extra_dll_dir): if os.path.isdir(extra_dll_dir):
os.environ["PATH"] += os.pathsep + extra_dll_dir os.environ["PATH"] += os.pathsep + extra_dll_dir
except ImportError:
pass
from . import api from . import api
from .api import * from .api import *
else:
raise
__all__ = [] __all__ = []
__all__.extend(api.__all__) __all__.extend(api.__all__)

1
src/wrf/api.py

@ -110,4 +110,3 @@ __all__ += ["CoordPair"]
__all__ += ["to_xy_coords"] __all__ += ["to_xy_coords"]
__all__ += ["cache_item", "get_cached_item"] __all__ += ["cache_item", "get_cached_item"]
__all__ += ["__version__"] __all__ += ["__version__"]

5
src/wrf/cache.py

@ -161,8 +161,3 @@ def _get_cache():
_shrink_cache() _shrink_cache()
return getattr(_local_storage, "cache", None) return getattr(_local_storage, "cache", None)

13
src/wrf/computation.py

@ -15,6 +15,7 @@ from .metadecorators import (set_alg_metadata, set_uvmet_alg_metadata,
set_smooth_metdata) set_smooth_metdata)
from .interputils import get_xy from .interputils import get_xy
@set_interp_metadata("xy") @set_interp_metadata("xy")
def xy(field, pivot_point=None, angle=None, start_point=None, end_point=None, def xy(field, pivot_point=None, angle=None, start_point=None, end_point=None,
meta=True): meta=True):
@ -601,8 +602,8 @@ def uvmet(u, v, lat, lon, cen_long, cone, meta=True, units="m s-1"):
and *v*, but with rightmost dimensions south_north x and *v*, but with rightmost dimensions south_north x
west_east and the same leftmost dimensions as *u* and *v* west_east and the same leftmost dimensions as *u* and *v*
- multi-dimensional with one fewer dimensions as *u* and *v*, - multi-dimensional with one fewer dimensions as *u* and *v*,
with rightmost dimensions south_north x west_east and the same with rightmost dimensions south_north x west_east and the
leftmost dimensions as *u* and *v*, minus the same leftmost dimensions as *u* and *v*, minus the
third-from-the-right dimension of *u* and *v*. third-from-the-right dimension of *u* and *v*.
Note: Note:
@ -622,8 +623,8 @@ def uvmet(u, v, lat, lon, cen_long, cone, meta=True, units="m s-1"):
and *v*, but with rightmost dimensions south_north x and *v*, but with rightmost dimensions south_north x
west_east and the same leftmost dimensions as *u* and *v* west_east and the same leftmost dimensions as *u* and *v*
- multi-dimensional with one fewer dimensions as *u* and *v*, - multi-dimensional with one fewer dimensions as *u* and *v*,
with rightmost dimensions south_north x west_east and the same with rightmost dimensions south_north x west_east and the
leftmost dimensions as *u* and *v*, minus the same leftmost dimensions as *u* and *v*, minus the
third-from-the-right dimension of *u* and *v*. third-from-the-right dimension of *u* and *v*.
@ -1201,7 +1202,6 @@ def ctt(pres_hpa, tkel, qv, qcld, height, terrain, qice=None,
return ma.masked_values(ctt, missing) return ma.masked_values(ctt, missing)
@set_alg_metadata(3, "pres", units="dBZ", @set_alg_metadata(3, "pres", units="dBZ",
description="radar reflectivity") description="radar reflectivity")
def dbz(pres, tkel, qv, qr, qs=None, qg=None, use_varint=False, def dbz(pres, tkel, qv, qr, qs=None, qg=None, use_varint=False,
@ -1955,6 +1955,3 @@ def pw(pres, tkel, qv, height, meta=True):
tv = _tv(tkel, qv) tv = _tv(tkel, qv)
return _pw(pres, tv, qv, height) return _pw(pres, tv, qv, height)

4
src/wrf/config.py

@ -9,6 +9,7 @@ from ._wrffortran import (fomp_enabled, fomp_set_num_threads,
_local_config = local() _local_config = local()
def _init_local(): def _init_local():
global _local_config global _local_config
@ -191,7 +192,8 @@ def set_cache_size(size):
@init_local() @init_local()
def get_cache_size(): def get_cache_size():
"""Return the maximum number of items that the threadlocal cache can retain. """Return the maximum number of items that the threadlocal cache can
retain.
Returns: Returns:

11
src/wrf/constants.py

@ -10,9 +10,11 @@ from ._wrffortran import wrf_constants, omp_constants
#: Indicates that all times should be used in a diagnostic routine. #: Indicates that all times should be used in a diagnostic routine.
ALL_TIMES = None ALL_TIMES = None
class Constants(object): class Constants(object):
pass pass
for key, val in viewitems(wrf_constants.__dict__): for key, val in viewitems(wrf_constants.__dict__):
setattr(Constants, key.upper(), np.asscalar(val)) setattr(Constants, key.upper(), np.asscalar(val))
@ -44,10 +46,11 @@ class ProjectionTypes(object):
MERCATOR = 3 MERCATOR = 3
LAT_LON = 6 LAT_LON = 6
# Create the default fill mapping based on type. # Create the default fill mapping based on type.
_DEFAULT_FILL_MAP = {None: Constants.DEFAULT_FILL, _DEFAULT_FILL_MAP = {None: Constants.DEFAULT_FILL,
np.dtype(np.bool_): False, np.dtype(np.bool_): False,
np.dtype(np.intc) : Constants.DEFAULT_FILL_INT32, # Usually true np.dtype(np.intc): Constants.DEFAULT_FILL_INT32,
np.dtype(np.int8): Constants.DEFAULT_FILL_INT8, np.dtype(np.int8): Constants.DEFAULT_FILL_INT8,
np.dtype(np.uint8): 255, np.dtype(np.uint8): 255,
np.dtype(np.int16): Constants.DEFAULT_FILL_INT16, np.dtype(np.int16): Constants.DEFAULT_FILL_INT16,
@ -76,9 +79,3 @@ else:
def default_fill(dtype=None): def default_fill(dtype=None):
dt = np.dtype(dtype) if dtype is not None else None dt = np.dtype(dtype) if dtype is not None else None
return _DEFAULT_FILL_MAP.get(dt, Constants.DEFAULT_FILL) return _DEFAULT_FILL_MAP.get(dt, Constants.DEFAULT_FILL)

31
src/wrf/coordpair.py

@ -35,13 +35,13 @@ def _binary_operator(operator):
""" """
if isinstance(other, CoordPair): if isinstance(other, CoordPair):
args = [ args = [None if getattr(self, attr) is None or
None if getattr(self, attr) is None or getattr(other, attr) is None getattr(other, attr) is None else
else getattr(getattr(self, attr), operator)(getattr(other, attr)) getattr(getattr(self, attr), operator)(getattr(other,
attr))
for attr in ("x", "y", "lat", "lon")] for attr in ("x", "y", "lat", "lon")]
else: else:
args = [ args = [None if getattr(self, attr) is None
None if getattr(self, attr) is None
else getattr(getattr(self, attr), operator)(other) else getattr(getattr(self, attr), operator)(other)
for attr in ("x", "y", "lat", "lon")] for attr in ("x", "y", "lat", "lon")]
@ -151,7 +151,6 @@ class CoordPair(object):
self.lat = lat self.lat = lat
self.lon = lon self.lon = lon
def __repr__(self): def __repr__(self):
args = [] args = []
if self.x is not None: if self.x is not None:
@ -166,11 +165,9 @@ class CoordPair(object):
return "{}({})".format(self.__class__.__name__, argstr) return "{}({})".format(self.__class__.__name__, argstr)
def __str__(self): def __str__(self):
return self.__repr__() return self.__repr__()
def xy_str(self, fmt="{:.4f}, {:.4f}"): def xy_str(self, fmt="{:.4f}, {:.4f}"):
"""Return a :obj:`str` for the (x,y) coordinate pair. """Return a :obj:`str` for the (x,y) coordinate pair.
@ -188,7 +185,6 @@ class CoordPair(object):
return fmt.format(self.x, self.y) return fmt.format(self.x, self.y)
def latlon_str(self, fmt="{:.4f}, {:.4f}"): def latlon_str(self, fmt="{:.4f}, {:.4f}"):
"""Return a :obj:`str` for the (latitude, longitude) coordinate pair. """Return a :obj:`str` for the (latitude, longitude) coordinate pair.
@ -206,7 +202,6 @@ class CoordPair(object):
return fmt.format(self.lat, self.lon) return fmt.format(self.lat, self.lon)
def __round__(self, ndigits=None): def __round__(self, ndigits=None):
"""Return a new :class:`CoordPair` object with all coordinate values """Return a new :class:`CoordPair` object with all coordinate values
rounded to the nearest integer. rounded to the nearest integer.
@ -226,23 +221,20 @@ class CoordPair(object):
return CoordPair(*args) return CoordPair(*args)
def __pow__(self, other, modulo=None): def __pow__(self, other, modulo=None):
if isinstance(other, CoordPair): if isinstance(other, CoordPair):
args = [ args = [None if getattr(self, attr) is None or
None if getattr(self, attr) is None or getattr(other, attr) is None getattr(other, attr) is None
else getattr(getattr(self, attr), "__pow__")(getattr(other, attr), else getattr(getattr(self, attr), "__pow__")(
modulo) getattr(other, attr), modulo)
for attr in ("x", "y", "lat", "lon")] for attr in ("x", "y", "lat", "lon")]
else: else:
args = [ args = [None if getattr(self, attr) is None
None if getattr(self, attr) is None
else getattr(getattr(self, attr), "__pow__")(other, modulo) else getattr(getattr(self, attr), "__pow__")(other, modulo)
for attr in ("x", "y", "lat", "lon")] for attr in ("x", "y", "lat", "lon")]
return CoordPair(*args) return CoordPair(*args)
def __rpow__(self, other): def __rpow__(self, other):
return self.__pow__(other) return self.__pow__(other)
@ -260,6 +252,3 @@ for operator in ("__neg__", "__pos__", "__abs__", "__invert__"):
for operator in ("__lt__", "__le__", "__eq__", "__ne__", "__gt__", "__ge__"): for operator in ("__lt__", "__le__", "__eq__", "__ne__", "__gt__", "__ge__"):
setattr(CoordPair, operator, _cmp_operator(operator)) setattr(CoordPair, operator, _cmp_operator(operator))

40
src/wrf/decorators.py

@ -15,6 +15,7 @@ from .constants import default_fill
if xarray_enabled(): if xarray_enabled():
from xarray import DataArray from xarray import DataArray
def convert_units(unit_type, alg_unit): def convert_units(unit_type, alg_unit):
"""A decorator that converts the units from the wrapped function's output. """A decorator that converts the units from the wrapped function's output.
@ -47,17 +48,6 @@ def convert_units(unit_type, alg_unit):
return func_wrapper return func_wrapper
#def _calc_out_dims(outvar, left_dims):
# """
#
# """
# #left_dims = [x for x in left_dims]
# #right_dims = [x for x in outvar.shape]
# #return left_dims + right_dims
#
# return left_dims + outvar.shape
def left_iteration(ref_var_expected_dims, def left_iteration(ref_var_expected_dims,
ref_var_right_ndims, ref_var_right_ndims,
insert_dims=None, insert_dims=None,
@ -171,7 +161,6 @@ def left_iteration(ref_var_expected_dims,
# the right (e.g. [1,1,:]) # the right (e.g. [1,1,:])
left_and_slice_idxs = left_idxs + (slice(None), ) left_and_slice_idxs = left_idxs + (slice(None), )
# Slice the args if applicable # Slice the args if applicable
new_args = [arg[left_and_slice_idxs] new_args = [arg[left_and_slice_idxs]
if i not in _ignore_args else arg if i not in _ignore_args else arg
@ -222,14 +211,11 @@ def left_iteration(ref_var_expected_dims,
outview.__array_interface__["data"][0]): outview.__array_interface__["data"][0]):
raise RuntimeError("output array was copied") raise RuntimeError("output array was copied")
if len(outd) == 1: if len(outd) == 1:
output = next(iter(viewvalues(outd))) output = next(iter(viewvalues(outd)))
else: else:
output = tuple(arr for arr in viewvalues(outd)) output = tuple(arr for arr in viewvalues(outd))
if cast_output: if cast_output:
if isinstance(output, np.ndarray): if isinstance(output, np.ndarray):
output = output.astype(ref_var_dtype) output = output.astype(ref_var_dtype)
@ -262,8 +248,8 @@ def cast_type(ref_idx=0, arg_idxs=None, karg_names=None,
positional arguments to be used as the reference variable for positional arguments to be used as the reference variable for
determining the :class:`numpy.dtype` to return. Default is 0. determining the :class:`numpy.dtype` to return. Default is 0.
arg_idxs (sequence of :obj:`int`, optional): A sequence of indexes in the arg_idxs (sequence of :obj:`int`, optional): A sequence of indexes in
wrapped function's positional arguments that indicate which the wrapped function's positional arguments that indicate which
arguments to cast. Must be specified if *karg_names* is None. arguments to cast. Must be specified if *karg_names* is None.
Default is None. Default is None.
@ -272,8 +258,8 @@ def cast_type(ref_idx=0, arg_idxs=None, karg_names=None,
arguments to cast. Must be specified if *arg_idxs* is None. arguments to cast. Must be specified if *arg_idxs* is None.
Default is None. Default is None.
alg_dtype (:class:`numpy.dtype` or :obj:`str`): The numpy data type used alg_dtype (:class:`numpy.dtype` or :obj:`str`): The numpy data type
in the wrapped function. used in the wrapped function.
outviews (:obj:`str` or a sequence): A single key or sequence of keys outviews (:obj:`str` or a sequence): A single key or sequence of keys
that indicate the wrapped function's keyword argument to use that indicate the wrapped function's keyword argument to use
@ -300,7 +286,6 @@ def cast_type(ref_idx=0, arg_idxs=None, karg_names=None,
if _outview is not None: if _outview is not None:
has_outview = True has_outview = True
orig_type = args[ref_idx].dtype orig_type = args[ref_idx].dtype
new_args = [arg.astype(alg_dtype) new_args = [arg.astype(alg_dtype)
@ -498,7 +483,8 @@ def check_args(refvaridx, refvarndim, rightdims, stagger=None,
# Check that the number of dims is correct # Check that the number of dims is correct
if (var.ndim - extra_dims != right_var_ndims): if (var.ndim - extra_dims != right_var_ndims):
raise ValueError("invalid number of dimensions for argument " raise ValueError("invalid number of dimensions for argument "
"{} (got {}, expected {}).".format(i, "{} (got {}, expected {}).".format(
i,
var.ndim, var.ndim,
right_var_ndims + extra_dims)) right_var_ndims + extra_dims))
@ -517,19 +503,11 @@ def check_args(refvaridx, refvarndim, rightdims, stagger=None,
ref_right_sizes[-right_var_ndims:]): ref_right_sizes[-right_var_ndims:]):
raise ValueError("invalid shape for argument " raise ValueError("invalid shape for argument "
"{} (got {}, expected {})".format(i, "{} (got {}, expected {})".format(
i,
var.shape[-right_var_ndims:], var.shape[-right_var_ndims:],
ref_right_sizes[-right_var_ndims:])) ref_right_sizes[-right_var_ndims:]))
return wrapped(*args, **kwargs) return wrapped(*args, **kwargs)
return func_wrapper return func_wrapper

2
src/wrf/destag.py

@ -60,5 +60,3 @@ def destagger(var, stagger_dim, meta=False):
result = .5*(var[tuple(dim_ranges_1)] + var[tuple(dim_ranges_2)]) result = .5*(var[tuple(dim_ranges_1)] + var[tuple(dim_ranges_2)])
return result return result

41
src/wrf/extension.py

@ -5,20 +5,21 @@ import numpy as np
from .constants import Constants, default_fill from .constants import Constants, default_fill
from wrf._wrffortran import (dcomputetk, dinterp3dz, dinterp2dxy, dinterp1d, from wrf._wrffortran import (dcomputetk, dinterp3dz, dinterp2dxy, dinterp1d,
dcomputeseaprs, dfilter2d, dcomputerh, dcomputeuvmet, dcomputeseaprs, dfilter2d, dcomputerh,
dcomputetd, dcapecalc2d, dcapecalc3d, dcloudfrac2, dcomputeuvmet, dcomputetd, dcapecalc2d,
wrfcttcalc, calcdbz, dcalrelhl, dcalcuh, dcomputepv, dcapecalc3d, dcloudfrac2, wrfcttcalc, calcdbz,
dcomputeabsvort, dlltoij, dijtoll, deqthecalc, dcalrelhl, dcalcuh, dcomputepv, dcomputeabsvort,
omgcalc, virtual_temp, wetbulbcalc, dcomputepw, dlltoij, dijtoll, deqthecalc, omgcalc,
virtual_temp, wetbulbcalc, dcomputepw,
wrf_monotonic, wrf_vintrp, dcomputewspd, wrf_monotonic, wrf_vintrp, dcomputewspd,
dcomputewdir, dinterp3dz_2dlev, dcomputewdir, dinterp3dz_2dlev,
fomp_set_num_threads, fomp_get_num_threads, fomp_set_num_threads, fomp_get_num_threads,
fomp_get_max_threads, fomp_get_thread_num, fomp_get_max_threads, fomp_get_thread_num,
fomp_get_num_procs, fomp_in_parallel, fomp_get_num_procs, fomp_in_parallel,
fomp_set_dynamic, fomp_get_dynamic, fomp_set_nested, fomp_set_dynamic, fomp_get_dynamic,
fomp_get_nested, fomp_set_schedule, fomp_set_nested, fomp_get_nested,
fomp_get_schedule, fomp_get_thread_limit, fomp_set_schedule, fomp_get_schedule,
fomp_set_max_active_levels, fomp_get_thread_limit, fomp_set_max_active_levels,
fomp_get_max_active_levels, fomp_get_level, fomp_get_max_active_levels, fomp_get_level,
fomp_get_ancestor_thread_num, fomp_get_team_size, fomp_get_ancestor_thread_num, fomp_get_team_size,
fomp_get_active_level, fomp_in_final, fomp_get_active_level, fomp_in_final,
@ -37,6 +38,7 @@ from .specialdec import (uvmet_left_iter, cape_left_iter,
cloudfrac_left_iter, check_cape_args, cloudfrac_left_iter, check_cape_args,
interplevel_left_iter, check_interplevel_args) interplevel_left_iter, check_interplevel_args)
class DiagnosticError(Exception): class DiagnosticError(Exception):
"""Raised when an error occurs in a diagnostic routine.""" """Raised when an error occurs in a diagnostic routine."""
def __init__(self, message=None): def __init__(self, message=None):
@ -66,6 +68,7 @@ class DiagnosticError(Exception):
""" """
raise self.__class__(message) raise self.__class__(message)
# The routines below are thin wrappers around the Fortran functions. These # The routines below are thin wrappers around the Fortran functions. These
# are not meant to be called by end users. Use the public API instead for # are not meant to be called by end users. Use the public API instead for
# that purpose. # that purpose.
@ -73,7 +76,6 @@ class DiagnosticError(Exception):
# IMPORTANT! Unless otherwise noted, all variables used in the routines # IMPORTANT! Unless otherwise noted, all variables used in the routines
# below assume that Fortran-ordered views are being used. This allows # below assume that Fortran-ordered views are being used. This allows
# f2py to pass the array pointers directly to the Fortran routine. # f2py to pass the array pointers directly to the Fortran routine.
@check_interplevel_args(is2dlev=False) @check_interplevel_args(is2dlev=False)
@interplevel_left_iter(is2dlev=False) @interplevel_left_iter(is2dlev=False)
@cast_type(arg_idxs=(0, 1, 2)) @cast_type(arg_idxs=(0, 1, 2))
@ -348,9 +350,8 @@ def _avo(u, v, msfu, msfv, msfm, cor, dx, dy, outview=None):
return result return result
@check_args(0, 3, (3,3,3,3,2,2,2,2), stagger=(-1,-2,None,None,-1,-2,None, @check_args(0, 3, (3, 3, 3, 3, 2, 2, 2, 2),
None), stagger=(-1, -2, None, None, -1, -2, None, None), refstagdim=-1)
refstagdim=-1)
@left_iteration(3, 3, ref_var_idx=2, ignore_args=(8, 9)) @left_iteration(3, 3, ref_var_idx=2, ignore_args=(8, 9))
@cast_type(arg_idxs=(0, 1, 2, 3, 4, 5, 6, 7)) @cast_type(arg_idxs=(0, 1, 2, 3, 4, 5, 6, 7))
@extract_and_transpose() @extract_and_transpose()
@ -403,8 +404,10 @@ def _eth(qv, tk, p, outview=None):
@cast_type(arg_idxs=(0, 1, 2, 3)) @cast_type(arg_idxs=(0, 1, 2, 3))
@extract_and_transpose() @extract_and_transpose()
def _uvmet(u, v, lat, lon, cen_long, cone, isstag=0, has_missing=False, def _uvmet(u, v, lat, lon, cen_long, cone, isstag=0, has_missing=False,
umissing=default_fill(np.float64), vmissing=default_fill(np.float64), umissing=default_fill(np.float64),
uvmetmissing=default_fill(np.float64), outview=None): vmissing=default_fill(np.float64),
uvmetmissing=default_fill(np.float64),
outview=None):
"""Wrapper for dcomputeuvmet. """Wrapper for dcomputeuvmet.
Located in wrf_user.f90. Located in wrf_user.f90.
@ -674,12 +677,14 @@ def _cape(p_hpa, tk, qv, ht, ter, sfp, missing, i3dflag, ter_follow,
return result return result
@check_args(0, 3, (3, 3)) @check_args(0, 3, (3, 3))
@cloudfrac_left_iter() @cloudfrac_left_iter()
@cast_type(arg_idxs=(0, 1), outviews=("lowview", "midview", "highview")) @cast_type(arg_idxs=(0, 1), outviews=("lowview", "midview", "highview"))
@extract_and_transpose(outviews=("lowview", "midview", "highview")) @extract_and_transpose(outviews=("lowview", "midview", "highview"))
def _cloudfrac(vert, rh, vert_inc_w_height, low_thresh, mid_thresh, def _cloudfrac(vert, rh, vert_inc_w_height, low_thresh, mid_thresh,
high_thresh, missing, lowview=None, midview=None, highview=None): high_thresh, missing, lowview=None, midview=None,
highview=None):
"""Wrapper for dcloudfrac2. """Wrapper for dcloudfrac2.
Located in wrf_cloud_fracf.f90. Located in wrf_cloud_fracf.f90.
@ -933,6 +938,7 @@ def _vintrp(field, pres, tk, qvp, ght, terrain, sfp, smsfp,
return result return result
@check_args(0, 2, (2, 2)) @check_args(0, 2, (2, 2))
@left_iteration(2, 2, ref_var_idx=0) @left_iteration(2, 2, ref_var_idx=0)
@cast_type(arg_idxs=(0, 1)) @cast_type(arg_idxs=(0, 1))
@ -1651,6 +1657,3 @@ def omp_get_wtick():
""" """
return fomp_get_wtick() return fomp_get_wtick()

5
src/wrf/g_cape.py

@ -9,6 +9,7 @@ from .constants import default_fill, Constants, ConversionFactors
from .util import extract_vars from .util import extract_vars
from .metadecorators import set_cape_metadata from .metadecorators import set_cape_metadata
@set_cape_metadata(is2d=True) @set_cape_metadata(is2d=True)
def get_2dcape(wrfin, timeidx=0, method="cat", squeeze=True, cache=None, def get_2dcape(wrfin, timeidx=0, method="cat", squeeze=True, cache=None,
meta=True, _key=None, missing=default_fill(np.float64)): meta=True, _key=None, missing=default_fill(np.float64)):
@ -221,7 +222,6 @@ def get_3dcape(wrfin, timeidx=0, method="cat",
cape_cin = _cape(p_hpa, tk, qv, z, ter, psfc_hpa, missing, i3dflag, cape_cin = _cape(p_hpa, tk, qv, z, ter, psfc_hpa, missing, i3dflag,
ter_follow) ter_follow)
return ma.masked_values(cape_cin, missing) return ma.masked_values(cape_cin, missing)
@ -629,6 +629,3 @@ def get_3dcin_only(wrfin, timeidx=0, method="cat",
result.attrs["units"] = "J kg-1" result.attrs["units"] = "J kg-1"
return result return result

20
src/wrf/g_cloudfrac.py

@ -164,8 +164,9 @@ def get_cloudfrac(wrfin, timeidx=0, method="cat", squeeze=True,
def get_low_cloudfrac(wrfin, timeidx=0, method="cat", squeeze=True, def get_low_cloudfrac(wrfin, timeidx=0, method="cat", squeeze=True,
cache=None, meta=True, _key=None, cache=None, meta=True, _key=None,
vert_type="height_agl", low_thresh=None, mid_thresh=None, vert_type="height_agl", low_thresh=None,
high_thresh=None, missing=default_fill(np.float64)): mid_thresh=None, high_thresh=None,
missing=default_fill(np.float64)):
"""Return the cloud fraction for the low level clouds. """Return the cloud fraction for the low level clouds.
If the vertical coordinate type is 'height_agl' or 'height_msl', the If the vertical coordinate type is 'height_agl' or 'height_msl', the
@ -264,9 +265,8 @@ def get_low_cloudfrac(wrfin, timeidx=0, method="cat", squeeze=True,
""" """
result = get_cloudfrac(wrfin, timeidx, method, squeeze, result = get_cloudfrac(wrfin, timeidx, method, squeeze,
cache, meta, _key, cache, meta, _key, vert_type, low_thresh,
vert_type, low_thresh, mid_thresh, mid_thresh, high_thresh, missing)[0, :]
high_thresh, missing)[0,:]
if meta: if meta:
result.attrs["description"] = "low clouds" result.attrs["description"] = "low clouds"
@ -276,8 +276,9 @@ def get_low_cloudfrac(wrfin, timeidx=0, method="cat", squeeze=True,
def get_mid_cloudfrac(wrfin, timeidx=0, method="cat", squeeze=True, def get_mid_cloudfrac(wrfin, timeidx=0, method="cat", squeeze=True,
cache=None, meta=True, _key=None, cache=None, meta=True, _key=None,
vert_type="height_agl", low_thresh=None, mid_thresh=None, vert_type="height_agl", low_thresh=None,
high_thresh=None, missing=default_fill(np.float64)): mid_thresh=None, high_thresh=None,
missing=default_fill(np.float64)):
"""Return the cloud fraction for the mid level clouds. """Return the cloud fraction for the mid level clouds.
If the vertical coordinate type is 'height_agl' or 'height_msl', the If the vertical coordinate type is 'height_agl' or 'height_msl', the
@ -388,8 +389,9 @@ def get_mid_cloudfrac(wrfin, timeidx=0, method="cat", squeeze=True,
def get_high_cloudfrac(wrfin, timeidx=0, method="cat", squeeze=True, def get_high_cloudfrac(wrfin, timeidx=0, method="cat", squeeze=True,
cache=None, meta=True, _key=None, cache=None, meta=True, _key=None,
vert_type="height_agl", low_thresh=None, mid_thresh=None, vert_type="height_agl", low_thresh=None,
high_thresh=None, missing=default_fill(np.float64)): mid_thresh=None, high_thresh=None,
missing=default_fill(np.float64)):
"""Return the cloud fraction for the high level clouds. """Return the cloud fraction for the high level clouds.
If the vertical coordinate type is 'height_agl' or 'height_msl', the If the vertical coordinate type is 'height_agl' or 'height_msl', the

1
src/wrf/g_ctt.py

@ -3,7 +3,6 @@ from __future__ import (absolute_import, division, print_function)
import numpy as np import numpy as np
import numpy.ma as ma import numpy.ma as ma
#from .extension import computectt, computetk
from .extension import _ctt, _tk from .extension import _ctt, _tk
from .constants import Constants, ConversionFactors, default_fill from .constants import Constants, ConversionFactors, default_fill
from .destag import destagger from .destag import destagger

2
src/wrf/g_dbz.py

@ -2,7 +2,6 @@ from __future__ import (absolute_import, division, print_function)
import numpy as np import numpy as np
#from .extension import computedbz,computetk
from .extension import _dbz, _tk from .extension import _dbz, _tk
from .constants import Constants from .constants import Constants
from .util import extract_vars, to_np from .util import extract_vars, to_np
@ -196,4 +195,3 @@ def get_max_dbz(wrfin, timeidx=0, method="cat",
return np.amax(to_np(get_dbz(wrfin, timeidx, method, squeeze, cache, meta, return np.amax(to_np(get_dbz(wrfin, timeidx, method, squeeze, cache, meta,
_key, use_varint, use_liqskin)), _key, use_varint, use_liqskin)),
axis=-3) axis=-3)

3
src/wrf/g_dewpoint.py

@ -1,6 +1,5 @@
from __future__ import (absolute_import, division, print_function) from __future__ import (absolute_import, division, print_function)
#from .extension import computetd
from .extension import _td from .extension import _td
from .decorators import convert_units from .decorators import convert_units
from .metadecorators import copy_and_set_metadata from .metadecorators import copy_and_set_metadata
@ -86,6 +85,7 @@ def get_dp(wrfin, timeidx=0, method="cat", squeeze=True,
td = _td(full_p, qvapor) td = _td(full_p, qvapor)
return td return td
@copy_and_set_metadata(copy_varname="Q2", name="td2", @copy_and_set_metadata(copy_varname="Q2", name="td2",
description="2m dew point temperature") description="2m dew point temperature")
@convert_units("temp", "c") @convert_units("temp", "c")
@ -161,4 +161,3 @@ def get_dp_2m(wrfin, timeidx=0, method="cat", squeeze=True,
td = _td(psfc, q2) td = _td(psfc, q2)
return td return td

5
src/wrf/g_geoht.py

@ -8,6 +8,7 @@ from .decorators import convert_units
from .metadecorators import set_height_metadata from .metadecorators import set_height_metadata
from .util import extract_vars, either from .util import extract_vars, either
def _get_geoht(wrfin, timeidx, method="cat", squeeze=True, def _get_geoht(wrfin, timeidx, method="cat", squeeze=True,
cache=None, meta=True, _key=None, cache=None, meta=True, _key=None,
height=True, msl=True, stag=False): height=True, msl=True, stag=False):
@ -321,8 +322,7 @@ def get_stag_geopt(wrfin, timeidx=0, method="cat", squeeze=True, cache=None,
@set_height_metadata(geopt=False, stag=True) @set_height_metadata(geopt=False, stag=True)
@convert_units("height", "m") @convert_units("height", "m")
def get_stag_height(wrfin, timeidx=0, method="cat", squeeze=True, def get_stag_height(wrfin, timeidx=0, method="cat", squeeze=True,
cache=None, meta=True, _key=None, cache=None, meta=True, _key=None, msl=True, units="m"):
msl=True, units="m"):
"""Return the geopotential height for the vertically staggered grid. """Return the geopotential height for the vertically staggered grid.
If *msl* is True, then geopotential height is returned as Mean Sea Level If *msl* is True, then geopotential height is returned as Mean Sea Level
@ -391,4 +391,3 @@ def get_stag_height(wrfin, timeidx=0, method="cat", squeeze=True,
return _get_geoht(wrfin, timeidx, method, squeeze, cache, meta, _key, return _get_geoht(wrfin, timeidx, method, squeeze, cache, meta, _key,
True, msl, stag=True) True, msl, stag=True)

7
src/wrf/g_helicity.py

@ -9,6 +9,7 @@ from .util import extract_vars, extract_global_attrs, either
from .metadecorators import copy_and_set_metadata from .metadecorators import copy_and_set_metadata
from .g_latlon import get_lat from .g_latlon import get_lat
@copy_and_set_metadata(copy_varname="HGT", name="srh", @copy_and_set_metadata(copy_varname="HGT", name="srh",
description="storm relative helicity", description="storm relative helicity",
units="m2 s-2") units="m2 s-2")
@ -109,6 +110,7 @@ def get_srh(wrfin, timeidx=0, method="cat", squeeze=True,
return srh return srh
@copy_and_set_metadata(copy_varname="MAPFAC_M", name="updraft_helicity", @copy_and_set_metadata(copy_varname="MAPFAC_M", name="updraft_helicity",
description="updraft helicity", description="updraft helicity",
units="m2 s-2") units="m2 s-2")
@ -206,8 +208,3 @@ def get_uh(wrfin, timeidx=0, method="cat", squeeze=True,
uh = _udhel(zp, mapfct, u, v, wstag, dx, dy, bottom, top) uh = _udhel(zp, mapfct, u, v, wstag, dx, dy, bottom, top)
return uh return uh

28
src/wrf/g_latlon.py

@ -301,12 +301,10 @@ def _llxy_mapping(wrfin, x_or_lat, y_or_lon, func, timeidx, stagger,
# make it so that key_0 is leftmost # make it so that key_0 is leftmost
outdims = key_coordnames + list(first_array.dims[existing_cnt:]) outdims = key_coordnames + list(first_array.dims[existing_cnt:])
# Create the new 'key_n', value pairs # Create the new 'key_n', value pairs
for coordname, coordval in zip(key_coordnames, coord_vals): for coordname, coordval in zip(key_coordnames, coord_vals):
outcoords[coordname] = coordval outcoords[coordname] = coordval
outattrs = OrderedDict(first_array.attrs) outattrs = OrderedDict(first_array.attrs)
outarr = DataArray(outdata, name=outname, coords=outcoords, outarr = DataArray(outdata, name=outname, coords=outcoords,
@ -605,10 +603,11 @@ def ll_to_xy_proj(latitude, longitude, meta=True, squeeze=True, as_int=True,
""" """
loc = locals() loc = locals()
_projparams = {name: loc[name] for name in ("map_proj", "truelat1", _projparams = {name: loc[name] for name in ("map_proj", "truelat1",
"truelat2", "stand_lon", "ref_lat", "truelat2", "stand_lon",
"ref_lon", "pole_lat", "pole_lon", "ref_lat", "ref_lon",
"known_x", "known_y", "dx", "dy", "pole_lat", "pole_lon",
"latinc", "loninc")} "known_x", "known_y", "dx",
"dy", "latinc", "loninc")}
projparams = _set_defaults(_projparams) projparams = _set_defaults(_projparams)
@ -773,19 +772,12 @@ def xy_to_ll_proj(x, y, meta=True, squeeze=True, map_proj=None, truelat1=None,
""" """
loc = locals() loc = locals()
_projparams = {name: loc[name] for name in ("map_proj", "truelat1", _projparams = {name: loc[name] for name in ("map_proj", "truelat1",
"truelat2", "stand_lon", "ref_lat", "truelat2", "stand_lon",
"ref_lon", "pole_lat", "pole_lon", "ref_lat", "ref_lon",
"known_x", "known_y", "dx", "dy", "pole_lat", "pole_lon",
"latinc", "loninc")} "known_x", "known_y", "dx",
"dy", "latinc", "loninc")}
projparams = _set_defaults(_projparams) projparams = _set_defaults(_projparams)
return _xy_to_ll(x, y, None, 0, None, "cat", squeeze, None, None, return _xy_to_ll(x, y, None, 0, None, "cat", squeeze, None, None,
**projparams) **projparams)

1
src/wrf/g_omega.py

@ -81,4 +81,3 @@ def get_omega(wrfin, timeidx=0, method="cat", squeeze=True, cache=None,
omega = _omega(qv, tk, wa, full_p) omega = _omega(qv, tk, wa, full_p)
return omega return omega

2
src/wrf/g_precip.py

@ -2,7 +2,6 @@ from __future__ import (absolute_import, division, print_function)
from .util import extract_vars from .util import extract_vars
__all__ = ["get_accum_precip", "get_precip_diff"]
def get_accum_precip(wrfin, timeidx=0): def get_accum_precip(wrfin, timeidx=0):
ncvars = extract_vars(wrfin, timeidx, varnames=("RAINC", "RAINNC")) ncvars = extract_vars(wrfin, timeidx, varnames=("RAINC", "RAINNC"))
@ -13,6 +12,7 @@ def get_accum_precip(wrfin, timeidx=0):
return rainsum return rainsum
def get_precip_diff(wrfin1, wrfin2, timeidx=0): def get_precip_diff(wrfin1, wrfin2, timeidx=0):
vars1 = extract_vars(wrfin1, timeidx, varnames=("RAINC", "RAINNC")) vars1 = extract_vars(wrfin1, timeidx, varnames=("RAINC", "RAINNC"))
vars2 = extract_vars(wrfin2, timeidx, varnames=("RAINC", "RAINNC")) vars2 = extract_vars(wrfin2, timeidx, varnames=("RAINC", "RAINNC"))

5
src/wrf/g_pressure.py

@ -85,6 +85,7 @@ def get_pressure(wrfin, timeidx=0, method="cat", squeeze=True,
return pres return pres
def get_pressure_hpa(wrfin, timeidx=0, method="cat", squeeze=True, def get_pressure_hpa(wrfin, timeidx=0, method="cat", squeeze=True,
cache=None, meta=True, _key=None): cache=None, meta=True, _key=None):
"""Return the pressure in [hPa]. """Return the pressure in [hPa].
@ -142,7 +143,3 @@ def get_pressure_hpa(wrfin, timeidx=0, method="cat", squeeze=True,
""" """
return get_pressure(wrfin, timeidx, method, squeeze, cache, meta, _key, return get_pressure(wrfin, timeidx, method, squeeze, cache, meta, _key,
units="hPa") units="hPa")

5
src/wrf/g_pw.py

@ -1,6 +1,5 @@
from __future__ import (absolute_import, division, print_function) from __future__ import (absolute_import, division, print_function)
#from .extension import computepw,computetv,computetk
from .extension import _pw, _tv, _tk from .extension import _pw, _tv, _tk
from .constants import Constants from .constants import Constants
from .util import extract_vars from .util import extract_vars
@ -85,7 +84,3 @@ def get_pw(wrfin, timeidx=0, method="cat", squeeze=True, cache=None,
tv = _tv(tk, qv) tv = _tv(tk, qv)
return _pw(full_p, tv, qv, ht) return _pw(full_p, tv, qv, ht)

2
src/wrf/g_rh.py

@ -1,7 +1,6 @@
from __future__ import (absolute_import, division, print_function) from __future__ import (absolute_import, division, print_function)
from .constants import Constants from .constants import Constants
#from .extension import computerh, computetk
from .extension import _rh, _tk from .extension import _rh, _tk
from .util import extract_vars from .util import extract_vars
from .metadecorators import copy_and_set_metadata from .metadecorators import copy_and_set_metadata
@ -153,4 +152,3 @@ def get_rh_2m(wrfin, timeidx=0, method="cat", squeeze=True, cache=None,
rh = _rh(q2, psfc, t2) rh = _rh(q2, psfc, t2)
return rh return rh

2
src/wrf/g_slp.py

@ -1,6 +1,5 @@
from __future__ import (absolute_import, division, print_function) from __future__ import (absolute_import, division, print_function)
#from .extension import computeslp, computetk
from .extension import _slp, _tk from .extension import _slp, _tk
from .constants import Constants from .constants import Constants
from .destag import destagger from .destag import destagger
@ -100,4 +99,3 @@ def get_slp(wrfin, timeidx=0, method="cat", squeeze=True,
slp = _slp(destag_ph, tk, full_p, qvapor) slp = _slp(destag_ph, tk, full_p, qvapor)
return slp return slp

3
src/wrf/g_temp.py

@ -511,6 +511,3 @@ def get_tc(wrfin, timeidx=0, method="cat", squeeze=True, cache=None,
""" """
return get_temp(wrfin, timeidx, method, squeeze, cache, meta, _key, return get_temp(wrfin, timeidx, method, squeeze, cache, meta, _key,
units="degC") units="degC")

3
src/wrf/g_terrain.py

@ -73,6 +73,3 @@ def get_terrain(wrfin, timeidx=0, method="cat", squeeze=True,
return extract_vars(wrfin, timeidx, varname, return extract_vars(wrfin, timeidx, varname,
method, squeeze, cache, meta=False, method, squeeze, cache, meta=False,
_key=_key)[varname] _key=_key)[varname]

1
src/wrf/g_times.py

@ -118,4 +118,3 @@ def get_xtimes(wrfin, timeidx=0, method="cat", squeeze=True, cache=None,
""" """
return extract_times(wrfin, timeidx, method, squeeze, cache, return extract_times(wrfin, timeidx, method, squeeze, cache,
meta=meta, do_xtime=True) meta=meta, do_xtime=True)

15
src/wrf/g_uvmet.py

@ -4,7 +4,6 @@ from math import fabs, log, tan, sin, cos
import numpy as np import numpy as np
#from .extension import computeuvmet
from .extension import _uvmet from .extension import _uvmet
from .destag import destagger from .destag import destagger
from .constants import Constants from .constants import Constants
@ -165,7 +164,6 @@ def _get_uvmet(wrfin, timeidx=0, method="cat", squeeze=True,
else: else:
cen_lon = lon_attrs["STAND_LON"] cen_lon = lon_attrs["STAND_LON"]
varname = either("XLAT_M", "XLAT")(wrfin) varname = either("XLAT_M", "XLAT")(wrfin)
xlat_var = extract_vars(wrfin, timeidx, varname, xlat_var = extract_vars(wrfin, timeidx, varname,
method, squeeze, cache, meta=False, method, squeeze, cache, meta=False,
@ -181,11 +179,12 @@ def _get_uvmet(wrfin, timeidx=0, method="cat", squeeze=True,
if map_proj == 1: if map_proj == 1:
if((fabs(true_lat1 - true_lat2) > 0.1) and if((fabs(true_lat1 - true_lat2) > 0.1) and
(fabs(true_lat2 - 90.) > 0.1)): (fabs(true_lat2 - 90.) > 0.1)):
cone = (log(cos(true_lat1*radians_per_degree)) cone = (log(cos(true_lat1*radians_per_degree)) -
- log(cos(true_lat2*radians_per_degree))) log(cos(true_lat2*radians_per_degree)))
cone = (cone / cone = (cone /
(log(tan((45.-fabs(true_lat1/2.))*radians_per_degree)) (log(tan((45.-fabs(true_lat1/2.))*radians_per_degree))
- log(tan((45.-fabs(true_lat2/2.))*radians_per_degree)))) - log(tan((45.-fabs(true_lat2/2.)) *
radians_per_degree))))
else: else:
cone = sin(fabs(true_lat1)*radians_per_degree) cone = sin(fabs(true_lat1)*radians_per_degree)
else: else:
@ -517,7 +516,8 @@ def get_uvmet10_wspd_wdir(wrfin, timeidx=0, method="cat", squeeze=True,
uvmet10 = _get_uvmet(wrfin, timeidx, method, squeeze, cache, meta, _key, uvmet10 = _get_uvmet(wrfin, timeidx, method, squeeze, cache, meta, _key,
True, units="m s-1") True, units="m s-1")
return _calc_wspd_wdir(uvmet10[0,...,:,:], uvmet10[1,...,:,:], True, units) return _calc_wspd_wdir(uvmet10[0, ..., :, :], uvmet10[1, ..., :, :],
True, units)
def get_uvmet_wspd(wrfin, timeidx=0, method="cat", squeeze=True, def get_uvmet_wspd(wrfin, timeidx=0, method="cat", squeeze=True,
@ -801,6 +801,3 @@ def get_uvmet10_wdir(wrfin, timeidx=0, method="cat", squeeze=True,
result.attrs["description"] = "10m earth rotated wdir" result.attrs["description"] = "10m earth rotated wdir"
return result return result

1
src/wrf/g_vorticity.py

@ -164,4 +164,3 @@ def get_pvo(wrfin, timeidx=0, method="cat", squeeze=True, cache=None,
full_p = p + pb full_p = p + pb
return _pvo(u, v, full_t, full_p, msfu, msfv, msfm, cor, dx, dy) return _pvo(u, v, full_t, full_p, msfu, msfv, msfm, cor, dx, dy)

7
src/wrf/g_wind.py

@ -93,12 +93,10 @@ def _calc_wspd_wdir(u, v, two_d, units):
result = np.zeros(outdims, wspd.dtype) result = np.zeros(outdims, wspd.dtype)
idxs0 = ((0, Ellipsis, slice(None), slice(None), slice(None)) idxs0 = ((0, Ellipsis, slice(None), slice(None), slice(None))
if not two_d else if not two_d else (0, Ellipsis, slice(None), slice(None)))
(0, Ellipsis, slice(None), slice(None)))
idxs1 = ((1, Ellipsis, slice(None), slice(None), slice(None)) idxs1 = ((1, Ellipsis, slice(None), slice(None), slice(None))
if not two_d else if not two_d else (1, Ellipsis, slice(None), slice(None)))
(1, Ellipsis, slice(None), slice(None)))
result[idxs0] = wspd[:] result[idxs0] = wspd[:]
result[idxs1] = wdir[:] result[idxs1] = wdir[:]
@ -784,4 +782,3 @@ def get_destag_wdir10(wrfin, timeidx=0, method="cat",
result.attrs["description"] = "10m wdir in projection space" result.attrs["description"] = "10m wdir in projection space"
return result return result

5
src/wrf/geobnds.py

@ -2,6 +2,7 @@ from __future__ import (absolute_import, division, print_function)
from .coordpair import CoordPair from .coordpair import CoordPair
class GeoBounds(object): class GeoBounds(object):
"""A class that stores the geographic boundaries. """A class that stores the geographic boundaries.
@ -84,7 +85,3 @@ class NullGeoBounds(GeoBounds):
def __repr__(self): def __repr__(self):
return "{}()".format(self.__class__.__name__) return "{}()".format(self.__class__.__name__)

65
src/wrf/interp.py

@ -76,11 +76,24 @@ def interplevel(field3d, vert, desiredlev, missing=default_fill(np.float64),
p = getvar(wrfin, "pressure") p = getvar(wrfin, "pressure")
ht = getvar(wrfin, "z", units="dm") ht = getvar(wrfin, "z", units="dm")
pblh = getvar(wrfin, "PBLH")
ht_500 = interplevel(ht, p, 500.0) ht_500 = interplevel(ht, p, 500.0)
Interpolate Relative Humidity to Boundary Layer Heights
.. code-block:: python
from netCDF4 import Dataset
from wrf import getvar, interplevel
wrfin = Dataset("wrfout_d02_2010-06-13_21:00:00")
rh = getvar(wrfin, "rh")
z = getvar(wrfin, "z")
pblh = getvar(wrfin, "PBLH")
rh_pblh = interplevel(rh, p, pblh)
""" """
@ -162,9 +175,9 @@ def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64),
timeidx (:obj:`int`, optional): The timeidx (:obj:`int`, optional): The
desired time index when obtaining map boundary information desired time index when obtaining map boundary information
from moving nests. This value can be a positive or negative integer. from moving nests. This value can be a positive or negative
Only required when *wrfin* is specified and the nest is moving. integer. Only required when *wrfin* is specified and the nest is
Currently, :data:`wrf.ALL_TIMES` is not supported. moving. Currently, :data:`wrf.ALL_TIMES` is not supported.
Default is 0. Default is 0.
stagger (:obj:`str`): If using latitude, longitude coordinate pairs stagger (:obj:`str`): If using latitude, longitude coordinate pairs
@ -286,12 +299,12 @@ def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64),
# domain could move outside of the line, which causes # domain could move outside of the line, which causes
# crashes or different line lengths. # crashes or different line lengths.
if is_moving: if is_moving:
raise ValueError("Requesting all times with a moving nest " raise ValueError("Requesting all times with a moving "
"is not supported when using lat/lon " "nest is not supported when using "
"cross sections because the domain could " "lat/lon cross sections because the "
"move outside of the cross section. " "domain could move outside of the "
"You must request each time " "cross section. You must request "
"individually.") "each time individually.")
else: else:
# Domain not moving, just use 0 # Domain not moving, just use 0
_timeidx = 0 _timeidx = 0
@ -379,9 +392,9 @@ def interpline(field2d, wrfin=None, timeidx=0, stagger=None, projection=None,
timeidx (:obj:`int`, optional): The timeidx (:obj:`int`, optional): The
desired time index when obtaining map boundary information desired time index when obtaining map boundary information
from moving nests. This value can be a positive or negative integer. from moving nests. This value can be a positive or negative
Only required when *wrfin* is specified and the nest is moving. integer. Only required when *wrfin* is specified and the nest is
Currently, :data:`wrf.ALL_TIMES` is not supported. moving. Currently, :data:`wrf.ALL_TIMES` is not supported.
Default is 0. Default is 0.
stagger (:obj:`str`): If using latitude, longitude coordinate pairs stagger (:obj:`str`): If using latitude, longitude coordinate pairs
@ -496,10 +509,10 @@ def interpline(field2d, wrfin=None, timeidx=0, stagger=None, projection=None,
# domain could move outside of the line, which causes # domain could move outside of the line, which causes
# crashes or different line lengths. # crashes or different line lengths.
if is_moving: if is_moving:
raise ValueError("Requesting all times with a moving nest " raise ValueError("Requesting all times with a moving "
"is not supported when using a lat/lon " "nest is not supported when using a "
"line because the domain could " "lat/lon line because the domain "
"move outside of line. " "could move outside of line. "
"You must request each time " "You must request each time "
"individually.") "individually.")
else: else:
@ -572,8 +585,8 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False,
* 'ght_msl': grid point height msl [km] * 'ght_msl': grid point height msl [km]
* 'ght_agl': grid point height agl [km] * 'ght_agl': grid point height agl [km]
* 'theta', 'th': potential temperature [K] * 'theta', 'th': potential temperature [K]
* 'theta-e', 'thetae', 'eth': equivalent potential temperature \ * 'theta-e', 'thetae', 'eth': equivalent potential \
[K] temperature [K]
interp_levels (sequence): A 1D sequence of vertical levels to interp_levels (sequence): A 1D sequence of vertical levels to
interpolate to. Values must be in the same units as specified interpolate to. Values must be in the same units as specified
@ -686,7 +699,6 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False,
"pressure_hpa": 1.0/ConversionFactors.PA_TO_HPA, "pressure_hpa": 1.0/ConversionFactors.PA_TO_HPA,
"z_km": 1.0/ConversionFactors.M_TO_KM, "z_km": 1.0/ConversionFactors.M_TO_KM,
"ght_km": 1.0/ConversionFactors.M_TO_KM, "ght_km": 1.0/ConversionFactors.M_TO_KM,
} }
out_unitmap = {"p_hpa": ConversionFactors.PA_TO_HPA, out_unitmap = {"p_hpa": ConversionFactors.PA_TO_HPA,
@ -694,7 +706,6 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False,
"pressure_hpa": ConversionFactors.PA_TO_HPA, "pressure_hpa": ConversionFactors.PA_TO_HPA,
"z_km": ConversionFactors.M_TO_KM, "z_km": ConversionFactors.M_TO_KM,
"ght_km": ConversionFactors.M_TO_KM, "ght_km": ConversionFactors.M_TO_KM,
} }
# These constants match what's in the fortran code. # These constants match what's in the fortran code.
@ -715,12 +726,12 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False,
# Check for valid coord # Check for valid coord
if vert_coord not in valid_coords: if vert_coord not in valid_coords:
raise ValueError("'%s' is not a valid vertical " raise ValueError("'{}' is not a valid vertical "
"coordinate type" % vert_coord) "coordinate type".format(vert_coord))
# Check for valid field type # Check for valid field type
if field_type not in valid_field_types: if field_type not in valid_field_types:
raise ValueError("'%s' is not a valid field type" % field_type) raise ValueError("'{}' is not a valid field type".format(field_type))
log_p_int = 1 if log_p else 0 log_p_int = 1 if log_p else 0
@ -848,9 +859,3 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False,
res_ = res res_ = res
return ma.masked_values(res_, missing) return ma.masked_values(res_, missing)

8
src/wrf/interputils.py

@ -85,7 +85,6 @@ def _calc_xy(xdim, ydim, pivot_point=None, angle=None,
if (angle > 315.0 or angle < 45.0 if (angle > 315.0 or angle < 45.0
or ((angle > 135.0) and (angle < 225.0))): or ((angle > 135.0) and (angle < 225.0))):
#x = y*slope + intercept
slope = -(360.-angle)/45. slope = -(360.-angle)/45.
if(angle < 45.): if(angle < 45.):
slope = angle/45. slope = angle/45.
@ -121,7 +120,7 @@ def _calc_xy(xdim, ydim, pivot_point=None, angle=None,
slope = (270.-angle)/45. slope = (270.-angle)/45.
intercept = yp - xp*slope intercept = yp - xp*slope
#find intersections with domain boundaries # Find intersections with domain boundaries
x0 = 0. x0 = 0.
y0 = x0*slope + intercept y0 = x0*slope + intercept
@ -323,6 +322,7 @@ def get_xy(var, pivot_point=None, angle=None,
return xy return xy
def to_xy_coords(pairs, wrfin=None, timeidx=0, stagger=None, projection=None, def to_xy_coords(pairs, wrfin=None, timeidx=0, stagger=None, projection=None,
ll_point=None): ll_point=None):
"""Return the coordinate pairs in grid space. """Return the coordinate pairs in grid space.
@ -389,7 +389,8 @@ def to_xy_coords(pairs, wrfin=None, timeidx=0, stagger=None, projection=None,
if wrfin is not None: if wrfin is not None:
xy_vals = _ll_to_xy(lat, lon, wrfin=wrfin, timeidx=timeidx, xy_vals = _ll_to_xy(lat, lon, wrfin=wrfin, timeidx=timeidx,
squeeze=True, meta=False, stagger=stagger, as_int=True) squeeze=True, meta=False, stagger=stagger,
as_int=True)
else: else:
map_proj = projection.map_proj map_proj = projection.map_proj
@ -427,7 +428,6 @@ def to_xy_coords(pairs, wrfin=None, timeidx=0, stagger=None, projection=None,
xy_vals = xy_vals.squeeze() xy_vals = xy_vals.squeeze()
if xy_vals.ndim == 1: if xy_vals.ndim == 1:
return CoordPair(x=xy_vals[0], y=xy_vals[1]) return CoordPair(x=xy_vals[0], y=xy_vals[1])
else: else:

23
src/wrf/latlonutils.py

@ -7,8 +7,8 @@ import numpy as np
from .constants import Constants, ProjectionTypes from .constants import Constants, ProjectionTypes
from .extension import _lltoxy, _xytoll from .extension import _lltoxy, _xytoll
from .util import (extract_vars, extract_global_attrs, from .util import (extract_vars, extract_global_attrs,
either, is_moving_domain, is_multi_time_req, either, is_moving_domain, iter_left_indexes,
iter_left_indexes, is_mapping, is_multi_file) is_mapping, is_multi_file)
from .py3compat import viewkeys, viewitems from .py3compat import viewkeys, viewitems
from .projutils import dict_keys_to_upper from .projutils import dict_keys_to_upper
@ -46,6 +46,7 @@ def _lat_varname(wrfin, stagger):
return varname return varname
def _lon_varname(wrfin, stagger): def _lon_varname(wrfin, stagger):
"""Return the longitude variable name for the specified stagger type. """Return the longitude variable name for the specified stagger type.
@ -79,6 +80,7 @@ def _lon_varname(wrfin, stagger):
return varname return varname
def _get_proj_params(wrfin, timeidx, stagger, method, squeeze, cache, _key): def _get_proj_params(wrfin, timeidx, stagger, method, squeeze, cache, _key):
"""Return the map projection parameters. """Return the map projection parameters.
@ -391,7 +393,6 @@ def _ll_to_xy(latitude, longitude, wrfin=None, timeidx=0,
result = np.empty(outdim, np.float64) result = np.empty(outdim, np.float64)
for left_idxs in iter_left_indexes(extra_dims): for left_idxs in iter_left_indexes(extra_dims):
#left_and_slice_idxs = left_idxs + (slice(None), )
# Left indexes is a misnomer, since these will be on the right # Left indexes is a misnomer, since these will be on the right
x_idxs = (0,) + left_idxs x_idxs = (0,) + left_idxs
y_idxs = (1,) + left_idxs y_idxs = (1,) + left_idxs
@ -426,7 +427,6 @@ def _ll_to_xy(latitude, longitude, wrfin=None, timeidx=0,
result[0] = fort_out[1] result[0] = fort_out[1]
result[1] = fort_out[0] result[1] = fort_out[0]
# Make indexes 0-based # Make indexes 0-based
result = result - 1 result = result - 1
@ -435,6 +435,7 @@ def _ll_to_xy(latitude, longitude, wrfin=None, timeidx=0,
return result return result
# X and Y should be 0-based # X and Y should be 0-based
def _xy_to_ll(x, y, wrfin=None, timeidx=0, stagger=None, def _xy_to_ll(x, y, wrfin=None, timeidx=0, stagger=None,
method="cat", squeeze=True, cache=None, _key=None, method="cat", squeeze=True, cache=None, _key=None,
@ -528,7 +529,6 @@ def _xy_to_ll(x, y, wrfin=None, timeidx=0, stagger=None,
pole_lat, pole_lon, known_x, known_y, dx, dy, latinc, pole_lat, pole_lon, known_x, known_y, dx, dy, latinc,
loninc) = _kwarg_proj_params(**projparams) loninc) = _kwarg_proj_params(**projparams)
if isinstance(x, Iterable): if isinstance(x, Iterable):
x_arr = np.asarray(x) x_arr = np.asarray(x)
y_arr = np.asarray(y) y_arr = np.asarray(y)
@ -555,7 +555,6 @@ def _xy_to_ll(x, y, wrfin=None, timeidx=0, stagger=None,
result = np.empty(outdim, np.float64) result = np.empty(outdim, np.float64)
for left_idxs in iter_left_indexes(extra_dims): for left_idxs in iter_left_indexes(extra_dims):
#left_and_slice_idxs = left_idxs + (slice(None), )
lat_idxs = (0,) + left_idxs lat_idxs = (0,) + left_idxs
lon_idxs = (1,) + left_idxs lon_idxs = (1,) + left_idxs
@ -573,7 +572,6 @@ def _xy_to_ll(x, y, wrfin=None, timeidx=0, stagger=None,
ref_lon_val, pole_lat, pole_lon, known_x, known_y, ref_lon_val, pole_lat, pole_lon, known_x, known_y,
dx, dy, latinc, loninc, x_val, y_val) dx, dy, latinc, loninc, x_val, y_val)
#result[left_and_slice_idxs] = ll[:]
result[lat_idxs] = ll[0] result[lat_idxs] = ll[0]
result[lon_idxs] = ll[1] result[lon_idxs] = ll[1]
@ -582,13 +580,8 @@ def _xy_to_ll(x, y, wrfin=None, timeidx=0, stagger=None,
x_val = x + 1 x_val = x + 1
y_val = y + 1 y_val = y + 1
result = _xytoll(map_proj, truelat1, truelat2, stdlon, ref_lat, ref_lon, result = _xytoll(map_proj, truelat1, truelat2, stdlon, ref_lat,
pole_lat, pole_lon, known_x, known_y, dx, dy, latinc, ref_lon, pole_lat, pole_lon, known_x, known_y,
loninc, x_val, y_val) dx, dy, latinc, loninc, x_val, y_val)
return result return result

99
src/wrf/metadecorators.py

@ -112,7 +112,6 @@ def copy_and_set_metadata(copy_varname=None, delete_attrs=None, name=None,
var_to_copy = None if cache is None else cache.get(_copy_varname, var_to_copy = None if cache is None else cache.get(_copy_varname,
None) None)
if var_to_copy is None: if var_to_copy is None:
var_to_copy = extract_vars(wrfin, timeidx, (_copy_varname,), var_to_copy = extract_vars(wrfin, timeidx, (_copy_varname,),
method, squeeze, cache, method, squeeze, cache,
@ -155,7 +154,6 @@ def copy_and_set_metadata(copy_varname=None, delete_attrs=None, name=None,
except KeyError: except KeyError:
pass pass
if name is not None: if name is not None:
outname = name outname = name
@ -277,10 +275,8 @@ def set_wind_metadata(copy_varname, name, description,
outattrs = OrderedDict() outattrs = OrderedDict()
outdimnames = list(copy_var.dims) outdimnames = list(copy_var.dims)
#outcoords.update(copy_var.coords)
outattrs.update(copy_var.attrs) outattrs.update(copy_var.attrs)
if wind_ncvar: if wind_ncvar:
outcoords.update(copy_var.coords) outcoords.update(copy_var.coords)
elif not wspd_wdir: elif not wspd_wdir:
@ -414,7 +410,6 @@ def set_cape_metadata(is2d):
outattrs["_FillValue"] = missing outattrs["_FillValue"] = missing
outattrs["missing_value"] = missing outattrs["missing_value"] = missing
# xarray doesn't line up coordinate dimensions based on # xarray doesn't line up coordinate dimensions based on
# names, it just remembers the index it originally mapped to. # names, it just remembers the index it originally mapped to.
# So, need to rebuild the XLAT, XLONG, coordinates again since the # So, need to rebuild the XLAT, XLONG, coordinates again since the
@ -433,7 +428,6 @@ def set_cape_metadata(is2d):
else: else:
outcoords["cape_cin"] = ["cape", "cin"] outcoords["cape_cin"] = ["cape", "cin"]
return DataArray(result, name=outname, coords=outcoords, return DataArray(result, name=outname, coords=outcoords,
dims=outdimnames, attrs=outattrs) dims=outdimnames, attrs=outattrs)
@ -632,8 +626,8 @@ def set_latlon_metadata(xy=False):
for x in zip(arr1, arr2)]) for x in zip(arr1, arr2)])
coords[dimnames[0]] = ["lat", "lon"] coords[dimnames[0]] = ["lat", "lon"]
else: else:
coords["latlon_coord"] = (dimnames[-1], [CoordPair(lat=x[0], coords["latlon_coord"] = (dimnames[-1],
lon=x[1]) [CoordPair(lat=x[0], lon=x[1])
for x in zip(arr1, arr2)]) for x in zip(arr1, arr2)])
coords[dimnames[0]] = ["x", "y"] coords[dimnames[0]] = ["x", "y"]
@ -750,11 +744,12 @@ def set_height_metadata(geopt=False, stag=False):
"(mass grid)".format(height_type)) "(mass grid)".format(height_type))
else: else:
outattrs["description"] = ("model height - [{}] (vertically " outattrs["description"] = ("model height - [{}] (vertically "
"staggered grid)".format(height_type)) "staggered grid)".format(
height_type))
return DataArray(result, name=outname, dims=outdimnames,
coords=outcoords, attrs=outattrs)
return DataArray(result, name=outname,
dims=outdimnames, coords=outcoords, attrs=outattrs)
return func_wrapper return func_wrapper
@ -819,7 +814,6 @@ def _set_horiz_meta(wrapped, instance, args, kwargs):
if isinstance(z, DataArray): if isinstance(z, DataArray):
vert_units = z.attrs.get("units", None) vert_units = z.attrs.get("units", None)
if isinstance(field3d, DataArray): if isinstance(field3d, DataArray):
outcoords = OrderedDict() outcoords = OrderedDict()
outdimnames = list(field3d.dims) outdimnames = list(field3d.dims)
@ -993,7 +987,8 @@ def _set_cross_meta(wrapped, instance, args, kwargs):
end_point_xy = (end_point.x, end_point.y) end_point_xy = (end_point.x, end_point.y)
xy, var2dz, z_var2d = get_xy_z_params(to_np(z), pivot_point_xy, angle, xy, var2dz, z_var2d = get_xy_z_params(to_np(z), pivot_point_xy, angle,
start_point_xy, end_point_xy, levels, autolevels) start_point_xy, end_point_xy,
levels, autolevels)
# Make a copy so we don't modify a user supplied cache # Make a copy so we don't modify a user supplied cache
if cache is not None: if cache is not None:
@ -1039,7 +1034,6 @@ def _set_cross_meta(wrapped, instance, args, kwargs):
except KeyError: except KeyError:
pass # Xarray 0.9 pass # Xarray 0.9
# Delete any lat,lon coords # Delete any lat,lon coords
delkeys = [key for key in viewkeys(outcoords) if is_coordvar(key)] delkeys = [key for key in viewkeys(outcoords) if is_coordvar(key)]
for key in delkeys: for key in delkeys:
@ -1069,12 +1063,10 @@ def _set_cross_meta(wrapped, instance, args, kwargs):
lats = _interpline(latcoord, xy) lats = _interpline(latcoord, xy)
lons = _interpline(loncoord, xy) lons = _interpline(loncoord, xy)
outcoords["xy_loc"] = ("cross_line_idx", outcoords["xy_loc"] = ("cross_line_idx", np.asarray(tuple(
np.asarray(tuple(
CoordPair(x=xy[i, 0], y=xy[i, 1], CoordPair(x=xy[i, 0], y=xy[i, 1],
lat=lats[i], lon=lons[i]) lat=lats[i], lon=lons[i])
for i in py3range(xy.shape[-2]))) for i in py3range(xy.shape[-2]))))
)
# Moving domain # Moving domain
else: else:
extra_dims = latcoord.shape[0:-2] extra_dims = latcoord.shape[0:-2]
@ -1086,12 +1078,11 @@ def _set_cross_meta(wrapped, instance, args, kwargs):
lats = _interpline(latcoord[idxs], xy) lats = _interpline(latcoord[idxs], xy)
lons = _interpline(loncoord[idxs], xy) lons = _interpline(loncoord[idxs], xy)
latlon_loc[idxs] = np.asarray(tuple( latlon_loc[idxs] = np.asarray(
CoordPair(x=xy[i,0], y=xy[i,1], tuple(CoordPair(
x=xy[i, 0], y=xy[i, 1],
lat=lats[i], lon=lons[i]) lat=lats[i], lon=lons[i])
for i in py3range(xy.shape[-2])) for i in py3range(xy.shape[-2])))[:]
)[:]
extra_dimnames = latcoord.dims[0:-2] extra_dimnames = latcoord.dims[0:-2]
loc_dimnames = extra_dimnames + ("cross_line_idx",) loc_dimnames = extra_dimnames + ("cross_line_idx",)
@ -1100,12 +1091,14 @@ def _set_cross_meta(wrapped, instance, args, kwargs):
else: else:
warnings.warn("'latlon' is set to True, but 'field3d' " warnings.warn("'latlon' is set to True, but 'field3d' "
"contains no coordinate information") "contains no coordinate information")
outcoords["xy_loc"] = ("cross_line_idx", np.asarray(tuple( outcoords["xy_loc"] = ("cross_line_idx",
np.asarray(tuple(
CoordPair(xy[i, 0], xy[i, 1]) CoordPair(xy[i, 0], xy[i, 1])
for i in py3range(xy.shape[-2])))) for i in py3range(xy.shape[-2]))))
else: else:
outcoords["xy_loc"] = ("cross_line_idx", np.asarray(tuple( outcoords["xy_loc"] = ("cross_line_idx",
np.asarray(tuple(
CoordPair(xy[i, 0], xy[i, 1]) CoordPair(xy[i, 0], xy[i, 1])
for i in py3range(xy.shape[-2])))) for i in py3range(xy.shape[-2]))))
@ -1160,11 +1153,10 @@ def _set_line_meta(wrapped, instance, args, kwargs):
:mod:`wrapt` :mod:`wrapt`
""" """
argvars = from_args(wrapped, ("field2d", argvars = from_args(wrapped, ("field2d", "wrfin", "timeidx", "stagger",
"wrfin", "timeidx", "stagger", "projection", "projection", "ll_point", "pivot_point",
"ll_point", "pivot_point", "angle", "angle", "start_point", "end_point",
"start_point", "end_point", "latlon", "latlon", "cache"),
"cache"),
*args, **kwargs) *args, **kwargs)
field2d = argvars["field2d"] field2d = argvars["field2d"]
@ -1223,7 +1215,6 @@ def _set_line_meta(wrapped, instance, args, kwargs):
# to avoid problems downstream # to avoid problems downstream
_timeidx = 0 _timeidx = 0
if pivot_point is not None: if pivot_point is not None:
if pivot_point.lat is not None and pivot_point.lon 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,
@ -1232,7 +1223,6 @@ def _set_line_meta(wrapped, instance, args, kwargs):
else: else:
pivot_point_xy = (pivot_point.x, pivot_point.y) pivot_point_xy = (pivot_point.x, pivot_point.y)
if start_point is not None and end_point is not None: 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: 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,
@ -1248,7 +1238,6 @@ def _set_line_meta(wrapped, instance, args, kwargs):
else: else:
end_point_xy = (end_point.x, end_point.y) end_point_xy = (end_point.x, end_point.y)
xy = get_xy(field2d, pivot_point_xy, angle, start_point_xy, end_point_xy) 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 # Make a copy so we don't modify a user supplied cache
@ -1276,8 +1265,7 @@ def _set_line_meta(wrapped, instance, args, kwargs):
cross_str = "({0}, {1}) to ({2}, {3})".format(st_x, st_y, ed_x, ed_y) cross_str = "({0}, {1}) to ({2}, {3})".format(st_x, st_y, ed_x, ed_y)
if angle is not None: if angle is not None:
cross_str += " ; center={0} ; angle={1}".format(pivot_point, cross_str += " ; center={0} ; angle={1}".format(pivot_point, angle)
angle)
if isinstance(field2d, DataArray): if isinstance(field2d, DataArray):
outcoords = OrderedDict() outcoords = OrderedDict()
@ -1318,12 +1306,10 @@ def _set_line_meta(wrapped, instance, args, kwargs):
lats = _interpline(latcoord, xy) lats = _interpline(latcoord, xy)
lons = _interpline(loncoord, xy) lons = _interpline(loncoord, xy)
outcoords["xy_loc"] = ("line_idx", outcoords["xy_loc"] = ("line_idx", np.asarray(tuple(
np.asarray(tuple(
CoordPair(x=xy[i, 0], y=xy[i, 1], CoordPair(x=xy[i, 0], y=xy[i, 1],
lat=lats[i], lon=lons[i]) lat=lats[i], lon=lons[i])
for i in py3range(xy.shape[-2]))) for i in py3range(xy.shape[-2]))))
)
# Moving domain # Moving domain
else: else:
@ -1339,9 +1325,7 @@ def _set_line_meta(wrapped, instance, args, kwargs):
latlon_loc[idxs] = np.asarray(tuple( latlon_loc[idxs] = np.asarray(tuple(
CoordPair(x=xy[i, 0], y=xy[i, 1], CoordPair(x=xy[i, 0], y=xy[i, 1],
lat=lats[i], lon=lons[i]) lat=lats[i], lon=lons[i])
for i in py3range(xy.shape[-2])) for i in py3range(xy.shape[-2])))[:]
)[:]
extra_dimnames = latcoord.dims[0:-2] extra_dimnames = latcoord.dims[0:-2]
loc_dimnames = extra_dimnames + ("line_idx",) loc_dimnames = extra_dimnames + ("line_idx",)
@ -1350,15 +1334,13 @@ def _set_line_meta(wrapped, instance, args, kwargs):
else: else:
warnings.warn("'latlon' is set to True, but 'field2d' " warnings.warn("'latlon' is set to True, but 'field2d' "
"contains no coordinate information") "contains no coordinate information")
outcoords["xy_loc"] = ("line_idx", np.asarray(tuple( outcoords["xy_loc"] = ("line_idx", np.asarray(
CoordPair(xy[i,0], xy[i,1]) tuple(CoordPair(xy[i, 0], xy[i, 1])
for i in py3range(xy.shape[-2])))) for i in py3range(xy.shape[-2]))))
else: else:
outcoords["xy_loc"] = ("line_idx", np.asarray(tuple( outcoords["xy_loc"] = ("line_idx", np.asarray(
CoordPair(xy[i,0], xy[i,1]) tuple(CoordPair(xy[i, 0], xy[i, 1])
for i in py3range(xy.shape[-2])))) for i in py3range(xy.shape[-2]))))
else: else:
if inc_latlon: if inc_latlon:
warnings.warn("'latlon' is set to True, but 'field2d' is " warnings.warn("'latlon' is set to True, but 'field2d' is "
@ -1427,7 +1409,6 @@ def _set_vinterp_meta(wrapped, instance, args, kwargs):
outcoords = None outcoords = None
outattrs = OrderedDict() outattrs = OrderedDict()
if isinstance(field, DataArray): if isinstance(field, DataArray):
outcoords = OrderedDict() outcoords = OrderedDict()
outdimnames = list(field.dims) outdimnames = list(field.dims)
@ -1443,7 +1424,6 @@ def _set_vinterp_meta(wrapped, instance, args, kwargs):
outcoords["interp_level"] = interp_levels outcoords["interp_level"] = interp_levels
outattrs.update(field.attrs) outattrs.update(field.attrs)
outname = field.name outname = field.name
else: else:
@ -1503,8 +1483,7 @@ def _set_2dxy_meta(wrapped, instance, args, kwargs):
ed_x = xy[-1, 0] ed_x = xy[-1, 0]
ed_y = xy[-1, 1] ed_y = xy[-1, 1]
cross_str = "({0},{1}) to ({2},{3})".format(st_x, st_y, cross_str = "({0},{1}) to ({2},{3})".format(st_x, st_y, ed_x, ed_y)
ed_x, ed_y)
outname = None outname = None
outdimnames = None outdimnames = None
@ -1532,7 +1511,6 @@ def _set_2dxy_meta(wrapped, instance, args, kwargs):
del outcoords[key] del outcoords[key]
outdimnames.append("line_idx") outdimnames.append("line_idx")
#outattrs.update(field3d.attrs)
desc = field3d.attrs.get("description", None) desc = field3d.attrs.get("description", None)
if desc is not None: if desc is not None:
@ -1544,7 +1522,8 @@ def _set_2dxy_meta(wrapped, instance, args, kwargs):
outname = "{0}_2dxy".format(field3d.name) outname = "{0}_2dxy".format(field3d.name)
outcoords["xy_loc"] = ("line_idx", [CoordPair(xy[i,0], xy[i,1]) outcoords["xy_loc"] = ("line_idx",
[CoordPair(xy[i, 0], xy[i, 1])
for i in py3range(xy.shape[-2])]) for i in py3range(xy.shape[-2])])
for key in ("MemoryOrder",): for key in ("MemoryOrder",):
@ -1824,7 +1803,6 @@ def set_alg_metadata(alg_ndims, refvarname,
if not xarray_enabled() or not do_meta: if not xarray_enabled() or not do_meta:
return wrapped(*args, **kwargs) return wrapped(*args, **kwargs)
result = wrapped(*args, **kwargs) result = wrapped(*args, **kwargs)
outname = wrapped.__name__ outname = wrapped.__name__
@ -1858,7 +1836,6 @@ def set_alg_metadata(alg_ndims, refvarname,
if _units is not None: if _units is not None:
outattrs["units"] = _units outattrs["units"] = _units
if description is not None: if description is not None:
if isinstance(description, from_var): if isinstance(description, from_var):
desc = description(wrapped, *args, **kwargs) desc = description(wrapped, *args, **kwargs)
@ -1867,7 +1844,6 @@ def set_alg_metadata(alg_ndims, refvarname,
else: else:
outattrs["description"] = description outattrs["description"] = description
# Copy the dimnames from the reference variable, otherwise, use # Copy the dimnames from the reference variable, otherwise, use
# the supplied dimnames # the supplied dimnames
if refvarname is not None: if refvarname is not None:
@ -2094,7 +2070,6 @@ def set_cape_alg_metadata(is2d, copyarg="pres_hpa"):
outattrs = OrderedDict() outattrs = OrderedDict()
if is2d: if is2d:
if is1d: if is1d:
outname = "cape_2d" outname = "cape_2d"
@ -2114,7 +2089,6 @@ def set_cape_alg_metadata(is2d, copyarg="pres_hpa"):
outattrs["description"] = "cape; cin" outattrs["description"] = "cape; cin"
outattrs["units"] = "J kg-1 ; J kg-1" outattrs["units"] = "J kg-1 ; J kg-1"
if isinstance(p, DataArray): if isinstance(p, DataArray):
if is2d: if is2d:
if not is1d: if not is1d:
@ -2131,7 +2105,6 @@ def set_cape_alg_metadata(is2d, copyarg="pres_hpa"):
else: else:
outdims[1] = p.dims[0] outdims[1] = p.dims[0]
outcoords = {} outcoords = {}
# Left-most is always cape_cin or cape_cin_lcl_lfc # Left-most is always cape_cin or cape_cin_lcl_lfc
if is2d: if is2d:
@ -2214,7 +2187,6 @@ def set_cloudfrac_alg_metadata(copyarg="vert"):
# Left dims # Left dims
outdims[1:-2] = cp.dims[0:-3] outdims[1:-2] = cp.dims[0:-3]
outcoords = {} outcoords = {}
# Left-most is always low_mid_high # Left-most is always low_mid_high
outdims[0] = "low_mid_high" outdims[0] = "low_mid_high"
@ -2294,8 +2266,3 @@ def set_destag_metadata():
return out return out
return func_wrapper return func_wrapper

74
src/wrf/projection.py

@ -61,7 +61,8 @@ if cartopy_enabled():
super(crs.Mercator, self).__init__(proj4_params, globe=globe) super(crs.Mercator, self).__init__(proj4_params, globe=globe)
# Calculate limits. # Calculate limits.
limits = self.transform_points(crs.Geodetic(), limits = self.transform_points(
crs.Geodetic(),
np.array([-180, 180]) + central_longitude, np.array([-180, 180]) + central_longitude,
np.array([min_latitude, max_latitude])) np.array([min_latitude, max_latitude]))
@ -195,7 +196,6 @@ class WrfProj(object):
if self.stand_lon is None: if self.stand_lon is None:
self.stand_lon = self._cen_lon self.stand_lon = self._cen_lon
@staticmethod @staticmethod
def _context_equal(x, y, ctx): def _context_equal(x, y, ctx):
"""Return True if both objects are equal based on the provided context. """Return True if both objects are equal based on the provided context.
@ -230,7 +230,6 @@ class WrfProj(object):
return True return True
def __eq__(self, other): def __eq__(self, other):
"""Return True if this projection object is the same as *other*. """Return True if this projection object is the same as *other*.
@ -262,7 +261,6 @@ class WrfProj(object):
WrfProj._context_equal(self.dx, other.dx, ctx) and WrfProj._context_equal(self.dx, other.dx, ctx) and
WrfProj._context_equal(self.dy, other.dy, ctx)) WrfProj._context_equal(self.dy, other.dy, ctx))
def _basemap(self, geobounds, **kwargs): def _basemap(self, geobounds, **kwargs):
return None return None
@ -275,17 +273,15 @@ class WrfProj(object):
def _calc_extents(self, geobounds): def _calc_extents(self, geobounds):
# Need to modify the extents for the new projection # Need to modify the extents for the new projection
pc = crs.PlateCarree() pc = crs.PlateCarree()
xs, ys, _ = self._cartopy().transform_points(pc, xs, ys, _ = self._cartopy().transform_points(
np.array([geobounds.bottom_left.lon, pc,
geobounds.top_right.lon]), np.array([geobounds.bottom_left.lon, geobounds.top_right.lon]),
np.array([geobounds.bottom_left.lat, np.array([geobounds.bottom_left.lat, geobounds.top_right.lat])).T
geobounds.top_right.lat])).T
_xlimits = xs.tolist() _xlimits = xs.tolist()
_ylimits = ys.tolist() _ylimits = ys.tolist()
return (_xlimits, _ylimits) return (_xlimits, _ylimits)
def _cart_extents(self, geobounds): def _cart_extents(self, geobounds):
try: try:
_ = len(geobounds) _ = len(geobounds)
@ -529,10 +525,9 @@ class LambertConformal(WrfProj):
if self.truelat2 is not None: if self.truelat2 is not None:
self._std_parallels.append(self.truelat2) self._std_parallels.append(self.truelat2)
def _cf_params(self): def _cf_params(self):
_cf_params = {} _cf_params = {}
_cf_params["grid_mapping_name"] = "lambert_conformal_conic"; _cf_params["grid_mapping_name"] = "lambert_conformal_conic"
_cf_params["standard_parallel"] = self._std_parallels _cf_params["standard_parallel"] = self._std_parallels
_cf_params["longitude_of_central_meridian"] = self.stand_lon _cf_params["longitude_of_central_meridian"] = self.stand_lon
_cf_params["latitude_of_projection_origin"] = self.moad_cen_lat _cf_params["latitude_of_projection_origin"] = self.moad_cen_lat
@ -540,13 +535,11 @@ class LambertConformal(WrfProj):
return _cf_params return _cf_params
def _pyngl(self, geobounds, **kwargs): def _pyngl(self, geobounds, **kwargs):
if not pyngl_enabled(): if not pyngl_enabled():
return None return None
truelat2 = (self.truelat1 truelat2 = (self.truelat1 if _ismissing(self.truelat2)
if _ismissing(self.truelat2)
else self.truelat2) else self.truelat2)
_pyngl = Resources() _pyngl = Resources()
@ -567,7 +560,6 @@ class LambertConformal(WrfProj):
return _pyngl return _pyngl
def _basemap(self, geobounds, **kwargs): def _basemap(self, geobounds, **kwargs):
if not basemap_enabled(): if not basemap_enabled():
return None return None
@ -650,16 +642,14 @@ class Mercator(WrfProj):
""" """
super(Mercator, self).__init__(**proj_params) super(Mercator, self).__init__(**proj_params)
self._lat_ts = (None self._lat_ts = (
if self.truelat1 == 0. or _ismissing(self.truelat1) None if self.truelat1 == 0. or _ismissing(self.truelat1)
else self.truelat1) else self.truelat1)
self._stand_lon = (0. if _ismissing(self.stand_lon, islat=False) self._stand_lon = (0. if _ismissing(self.stand_lon, islat=False)
else self.stand_lon) else self.stand_lon)
def _cf_params(self): def _cf_params(self):
_cf_params = {} _cf_params = {}
_cf_params["grid_mapping_name"] = "mercator" _cf_params["grid_mapping_name"] = "mercator"
_cf_params["longitude_of_projection_origin"] = self.stand_lon _cf_params["longitude_of_projection_origin"] = self.stand_lon
@ -667,7 +657,6 @@ class Mercator(WrfProj):
return _cf_params return _cf_params
def _pyngl(self, geobounds, **kwargs): def _pyngl(self, geobounds, **kwargs):
if not pyngl_enabled(): if not pyngl_enabled():
return None return None
@ -689,7 +678,6 @@ class Mercator(WrfProj):
return _pyngl return _pyngl
def _basemap(self, geobounds, **kwargs): def _basemap(self, geobounds, **kwargs):
if not basemap_enabled(): if not basemap_enabled():
return None return None
@ -710,25 +698,20 @@ class Mercator(WrfProj):
return _basemap return _basemap
def _cartopy(self): def _cartopy(self):
if not cartopy_enabled(): if not cartopy_enabled():
return None return None
if self._lat_ts == 0.0: if self._lat_ts == 0.0:
_cartopy = crs.Mercator( _cartopy = crs.Mercator(central_longitude=self._stand_lon,
central_longitude = self._stand_lon,
globe=self._globe()) globe=self._globe())
else: else:
_cartopy = MercatorWithLatTS( _cartopy = MercatorWithLatTS(central_longitude=self._stand_lon,
central_longitude = self._stand_lon,
latitude_true_scale=self._lat_ts, latitude_true_scale=self._lat_ts,
globe=self._globe()) globe=self._globe())
return _cartopy return _cartopy
def _proj4(self): def _proj4(self):
_proj4 = ("+proj=merc +units=meters +a={} +b={} " _proj4 = ("+proj=merc +units=meters +a={} +b={} "
@ -740,6 +723,7 @@ class Mercator(WrfProj):
return _proj4 return _proj4
class PolarStereographic(WrfProj): class PolarStereographic(WrfProj):
"""A :class:`wrf.WrfProj` subclass for Polar Stereographic projections. """A :class:`wrf.WrfProj` subclass for Polar Stereographic projections.
@ -750,7 +734,6 @@ class PolarStereographic(WrfProj):
:class:`Mercator`, :class:`LambertConformal` :class:`Mercator`, :class:`LambertConformal`
""" """
def __init__(self, **proj_params): def __init__(self, **proj_params):
"""Initialize a :class:`wrf.PolarStereographic` object. """Initialize a :class:`wrf.PolarStereographic` object.
@ -770,10 +753,7 @@ class PolarStereographic(WrfProj):
""" """
super(PolarStereographic, self).__init__(**proj_params) super(PolarStereographic, self).__init__(**proj_params)
self._hemi = -90. if self.truelat1 < 0 else 90. self._hemi = -90. if self.truelat1 < 0 else 90.
self._lat_ts = (None self._lat_ts = (None if _ismissing(self.truelat1) else self.truelat1)
if _ismissing(self.truelat1)
else self.truelat1)
def _cf_params(self): def _cf_params(self):
_cf_params = {} _cf_params = {}
@ -785,7 +765,6 @@ class PolarStereographic(WrfProj):
return _cf_params return _cf_params
def _pyngl(self, geobounds, **kwargs): def _pyngl(self, geobounds, **kwargs):
if not pyngl_enabled(): if not pyngl_enabled():
return None return None
@ -811,7 +790,6 @@ class PolarStereographic(WrfProj):
return _pyngl return _pyngl
def _basemap(self, geobounds, **kwargs): def _basemap(self, geobounds, **kwargs):
if not basemap_enabled(): if not basemap_enabled():
return None return None
@ -832,7 +810,6 @@ class PolarStereographic(WrfProj):
return _basemap return _basemap
def _cartopy(self): def _cartopy(self):
if not cartopy_enabled(): if not cartopy_enabled():
return None return None
@ -843,7 +820,6 @@ class PolarStereographic(WrfProj):
globe=self._globe()) globe=self._globe())
return _cartopy return _cartopy
def _proj4(self): def _proj4(self):
_proj4 = ("+proj=stere +units=meters +a={} +b={} " _proj4 = ("+proj=stere +units=meters +a={} +b={} "
"+lat0={} +lon_0={} +lat_ts={} +nadgrids=@null".format( "+lat0={} +lon_0={} +lat_ts={} +nadgrids=@null".format(
@ -856,7 +832,6 @@ class PolarStereographic(WrfProj):
return _proj4 return _proj4
class LatLon(WrfProj): class LatLon(WrfProj):
"""A :class:`wrf.WrfProj` subclass for Lat Lon projections. """A :class:`wrf.WrfProj` subclass for Lat Lon projections.
@ -886,13 +861,11 @@ class LatLon(WrfProj):
""" """
super(LatLon, self).__init__(**proj_params) super(LatLon, self).__init__(**proj_params)
def _cf_params(self): def _cf_params(self):
_cf_params = {} _cf_params = {}
_cf_params["grid_mapping_name"] = "latitude_longitude" _cf_params["grid_mapping_name"] = "latitude_longitude"
return _cf_params return _cf_params
def _pyngl(self, geobounds, **kwargs): def _pyngl(self, geobounds, **kwargs):
if not pyngl_enabled(): if not pyngl_enabled():
return None return None
@ -914,7 +887,6 @@ class LatLon(WrfProj):
return _pyngl return _pyngl
def _basemap(self, geobounds, **kwargs): def _basemap(self, geobounds, **kwargs):
if not basemap_enabled(): if not basemap_enabled():
return None return None
@ -934,7 +906,6 @@ class LatLon(WrfProj):
return _basemap return _basemap
def _cartopy(self): def _cartopy(self):
if not cartopy_enabled(): if not cartopy_enabled():
return None return None
@ -944,19 +915,19 @@ class LatLon(WrfProj):
return _cartopy return _cartopy
def _cart_extents(self, geobounds): def _cart_extents(self, geobounds):
return ([geobounds.bottom_left.lon, geobounds.top_right.lon], return ([geobounds.bottom_left.lon, geobounds.top_right.lon],
[geobounds.bottom_left.lat, geobounds.top_right.lat]) [geobounds.bottom_left.lat, geobounds.top_right.lat])
def _proj4(self): def _proj4(self):
_proj4 = ("+proj=eqc +units=meters +a={} +b={} " _proj4 = ("+proj=eqc +units=meters +a={} +b={} "
"+lon_0={} +nadgrids=@null".format(Constants.WRF_EARTH_RADIUS, "+lon_0={} +nadgrids=@null".format(
Constants.WRF_EARTH_RADIUS,
Constants.WRF_EARTH_RADIUS, Constants.WRF_EARTH_RADIUS,
self.stand_lon)) self.stand_lon))
return _proj4 return _proj4
# Notes (may not be correct since this projection confuses me): # Notes (may not be correct since this projection confuses me):
# Each projection system handles this differently. # Each projection system handles this differently.
# 1) In WRF, if following the WPS instructions, POLE_LON is mainly used to # 1) In WRF, if following the WPS instructions, POLE_LON is mainly used to
@ -1059,7 +1030,6 @@ class RotatedLatLon(WrfProj):
self._cart_pole_lon = (-self.stand_lon - 180.0 if self._north self._cart_pole_lon = (-self.stand_lon - 180.0 if self._north
else -self.stand_lon) else -self.stand_lon)
def _cf_params(self): def _cf_params(self):
_cf_params = {} _cf_params = {}
# Assuming this follows the same guidelines as cartopy # Assuming this follows the same guidelines as cartopy
@ -1070,7 +1040,6 @@ class RotatedLatLon(WrfProj):
return _cf_params return _cf_params
def _pyngl(self, geobounds, **kwargs): def _pyngl(self, geobounds, **kwargs):
if not pyngl_enabled(): if not pyngl_enabled():
return None return None
@ -1092,7 +1061,6 @@ class RotatedLatLon(WrfProj):
return _pyngl return _pyngl
def _basemap(self, geobounds, **kwargs): def _basemap(self, geobounds, **kwargs):
if not basemap_enabled(): if not basemap_enabled():
return None return None
@ -1113,13 +1081,11 @@ class RotatedLatLon(WrfProj):
return _basemap return _basemap
def _cartopy(self): def _cartopy(self):
if not cartopy_enabled(): if not cartopy_enabled():
return None return None
_cartopy = crs.RotatedPole( _cartopy = crs.RotatedPole(pole_longitude=self._cart_pole_lon,
pole_longitude=self._cart_pole_lon,
pole_latitude=self._bm_cart_pole_lat, pole_latitude=self._bm_cart_pole_lat,
central_rotated_longitude=( central_rotated_longitude=(
180.0 - self.pole_lon), # Probably 180.0 - self.pole_lon), # Probably
@ -1138,6 +1104,7 @@ class RotatedLatLon(WrfProj):
return _proj4 return _proj4
def getproj(**proj_params): def getproj(**proj_params):
"""Return a :class:`wrf.WrfProj` subclass. """Return a :class:`wrf.WrfProj` subclass.
@ -1164,7 +1131,6 @@ def getproj(**proj_params):
specified map projection parameters. specified map projection parameters.
""" """
up_proj_params = dict_keys_to_upper(proj_params) up_proj_params = dict_keys_to_upper(proj_params)
proj_type = up_proj_params.get("MAP_PROJ", 0) proj_type = up_proj_params.get("MAP_PROJ", 0)
@ -1184,5 +1150,3 @@ def getproj(**proj_params):
else: else:
# Unknown projection # Unknown projection
return WrfProj(**proj_params) return WrfProj(**proj_params)

1
src/wrf/projutils.py

@ -2,6 +2,7 @@ from __future__ import (absolute_import, division, print_function)
from .py3compat import viewitems from .py3compat import viewitems
def dict_keys_to_upper(d): def dict_keys_to_upper(d):
"""Return a dictionary with the keys changed to uppercase. """Return a dictionary with the keys changed to uppercase.

4
src/wrf/py3compat.py

@ -3,6 +3,7 @@ from __future__ import (absolute_import, division, print_function)
from sys import version_info from sys import version_info
from math import floor, copysign from math import floor, copysign
# Dictionary python 2-3 compatibility stuff # Dictionary python 2-3 compatibility stuff
def viewitems(d): def viewitems(d):
"""Return either the items or viewitems method for a dictionary. """Return either the items or viewitems method for a dictionary.
@ -57,6 +58,7 @@ def viewvalues(d):
func = d.values func = d.values
return func() return func()
def isstr(s): def isstr(s):
"""Return True if the object is a string type. """Return True if the object is a string type.
@ -165,5 +167,3 @@ def ucode(*args, **kwargs):
return str(*args, **kwargs) return str(*args, **kwargs)
return unicode(*args, **kwargs) return unicode(*args, **kwargs)

9
src/wrf/routines.py

@ -16,7 +16,8 @@ from .g_pressure import get_pressure, get_pressure_hpa
from .g_pw import get_pw from .g_pw import get_pw
from .g_rh import get_rh, get_rh_2m from .g_rh import get_rh, get_rh_2m
from .g_slp import get_slp from .g_slp import get_slp
from .g_temp import get_tc, get_eth, get_temp, get_theta, get_tk, get_tv, get_tw from .g_temp import (get_tc, get_eth, get_temp, get_theta, get_tk, get_tv,
get_tw)
from .g_terrain import get_terrain from .g_terrain import get_terrain
from .g_uvmet import (get_uvmet, get_uvmet10, get_uvmet10_wspd_wdir, from .g_uvmet import (get_uvmet, get_uvmet10, get_uvmet10_wspd_wdir,
get_uvmet_wspd_wdir, get_uvmet_wspd, get_uvmet_wdir, get_uvmet_wspd_wdir, get_uvmet_wspd, get_uvmet_wdir,
@ -194,6 +195,7 @@ _ALIASES = {"cape_2d" : "cape2d",
"wdir_uvmet10": "uvmet10_wdir", "wdir_uvmet10": "uvmet10_wdir",
} }
class ArgumentError(Exception): class ArgumentError(Exception):
def __init__(self, msg): def __init__(self, msg):
self.msg = msg self.msg = msg
@ -201,6 +203,7 @@ class ArgumentError(Exception):
def __str__(self): def __str__(self):
return self.msg return self.msg
def _undo_alias(alias): def _undo_alias(alias):
actual = _ALIASES.get(alias, None) actual = _ALIASES.get(alias, None)
if actual is None: if actual is None:
@ -208,6 +211,7 @@ def _undo_alias(alias):
else: else:
return actual return actual
def _check_kargs(var, kargs): def _check_kargs(var, kargs):
for arg in viewkeys(kargs): for arg in viewkeys(kargs):
if arg not in _VALID_KARGS[var]: if arg not in _VALID_KARGS[var]:
@ -340,10 +344,9 @@ def getvar(wrfin, varname, timeidx=0,
actual_var = _undo_alias(varname) actual_var = _undo_alias(varname)
if actual_var not in _VALID_KARGS: if actual_var not in _VALID_KARGS:
raise ValueError("'%s' is not a valid variable name" % (varname)) raise ValueError("'{}' is not a valid variable name".format(varname))
_check_kargs(actual_var, kwargs) _check_kargs(actual_var, kwargs)
return _FUNC_MAP[actual_var](wrfin, timeidx, method, squeeze, cache, return _FUNC_MAP[actual_var](wrfin, timeidx, method, squeeze, cache,
meta, _key, **kwargs) meta, _key, **kwargs)

15
src/wrf/specialdec.py

@ -99,8 +99,6 @@ def uvmet_left_iter(alg_dtype=np.float64):
if (u.shape[-1] == lat.shape[-1] or u.shape[-2] == lat.shape[-2]): if (u.shape[-1] == lat.shape[-1] or u.shape[-2] == lat.shape[-2]):
raise ValueError("v is staggered but u is not") raise ValueError("v is staggered but u is not")
# No special left side iteration, return the function result # No special left side iteration, return the function result
if (num_left_dims_u == 0): if (num_left_dims_u == 0):
return wrapped(u, v, lat, lon, cen_long, cone, isstag=is_stag, return wrapped(u, v, lat, lon, cen_long, cone, isstag=is_stag,
@ -132,15 +130,14 @@ def uvmet_left_iter(alg_dtype=np.float64):
lat_left_and_slice = left_and_slice_idxs lat_left_and_slice = left_and_slice_idxs
elif mode == 2: elif mode == 2:
# Only need the left-most # Only need the left-most
lat_left_and_slice = tuple(left_idx lat_left_and_slice = tuple(
for left_idx in left_idxs[0:num_left_dims_lat]) left_idx for left_idx in left_idxs[0:num_left_dims_lat])
u_output_idxs = (0,) + left_idxs + (slice(None),) u_output_idxs = (0,) + left_idxs + (slice(None),)
v_output_idxs = (1,) + left_idxs + (slice(None),) v_output_idxs = (1,) + left_idxs + (slice(None),)
u_view_idxs = left_idxs + (0, slice(None)) u_view_idxs = left_idxs + (0, slice(None))
v_view_idxs = left_idxs + (1, slice(None)) v_view_idxs = left_idxs + (1, slice(None))
new_u = u[left_and_slice_idxs] new_u = u[left_and_slice_idxs]
new_v = v[left_and_slice_idxs] new_v = v[left_and_slice_idxs]
new_lat = lat[lat_left_and_slice] new_lat = lat[lat_left_and_slice]
@ -186,7 +183,6 @@ def uvmet_left_iter(alg_dtype=np.float64):
return func_wrapper return func_wrapper
def cape_left_iter(alg_dtype=np.float64): def cape_left_iter(alg_dtype=np.float64):
"""A decorator to handle iterating over the leftmost dimensions for the """A decorator to handle iterating over the leftmost dimensions for the
cape diagnostic. cape diagnostic.
@ -374,7 +370,6 @@ def cape_left_iter(alg_dtype=np.float64):
capeview.__array_interface__["data"][0]): capeview.__array_interface__["data"][0]):
raise RuntimeError("output array was copied") raise RuntimeError("output array was copied")
if flip and not is2d: if flip and not is2d:
output[cape_output_idxs] = ( output[cape_output_idxs] = (
outview_array[view_cape_reverse_idxs].astype(orig_dtype)) outview_array[view_cape_reverse_idxs].astype(orig_dtype))
@ -441,8 +436,8 @@ def cloudfrac_left_iter(alg_dtype=np.float64):
return output return output
# Initial output is ...,low_mid_high,nz,ny,nx to create contiguous
# Initial output is ...,low_mid_high,nz,ny,nx to create contiguous views # views
outdims = vert.shape[0:num_left_dims] outdims = vert.shape[0:num_left_dims]
extra_dims = tuple(outdims) # Copy the left-most dims for iteration extra_dims = tuple(outdims) # Copy the left-most dims for iteration
@ -579,7 +574,6 @@ def interplevel_left_iter(is2dlev, alg_dtype=np.float64):
else: else:
z_slice_idxs = left_idxs + (slice(None),) z_slice_idxs = left_idxs + (slice(None),)
new_args[0] = field3d[field_out_slice_idxs] new_args[0] = field3d[field_out_slice_idxs]
new_args[1] = z[z_slice_idxs] new_args[1] = z[z_slice_idxs]
@ -695,4 +689,3 @@ def check_interplevel_args(is2dlev):
return wrapped(*args, **kwargs) return wrapped(*args, **kwargs)
return func_wrapper return func_wrapper

14
src/wrf/units.py

@ -152,7 +152,6 @@ _UNIT_ALIASES = {"mps" : "m s-1",
"fps": "ft s-1", "fps": "ft s-1",
"fs-1": "ft s-1", "fs-1": "ft s-1",
"feet_per_second": "ft s-1", "feet_per_second": "ft s-1",
"pascal": "pa", "pascal": "pa",
"pascals": "pa", "pascals": "pa",
"hecto_pascal": "hpa", "hecto_pascal": "hpa",
@ -160,7 +159,6 @@ _UNIT_ALIASES = {"mps" : "m s-1",
"millibar": "mb", "millibar": "mb",
"millibars": "mb", "millibars": "mb",
"mbar": "mb", "mbar": "mb",
"kelvin": "k", "kelvin": "k",
"degree_kelvin": "k", "degree_kelvin": "k",
"degrees_kelvin": "k", "degrees_kelvin": "k",
@ -174,7 +172,6 @@ _UNIT_ALIASES = {"mps" : "m s-1",
"degs_k": "k", "degs_k": "k",
"deg k": "k", "deg k": "k",
"degs k": "k", "degs k": "k",
"celsius": "c", "celsius": "c",
"degree_celsius": "c", "degree_celsius": "c",
"degrees_celsius": "c", "degrees_celsius": "c",
@ -188,7 +185,6 @@ _UNIT_ALIASES = {"mps" : "m s-1",
"degs_c": "c", "degs_c": "c",
"deg c": "c", "deg c": "c",
"degs c": "c", "degs c": "c",
"fahrenheit": "f", "fahrenheit": "f",
"degree_fahrenheit": "f", "degree_fahrenheit": "f",
"degrees_fahrenheit": "f", "degrees_fahrenheit": "f",
@ -202,7 +198,6 @@ _UNIT_ALIASES = {"mps" : "m s-1",
"degs_f": "f", "degs_f": "f",
"deg f": "f", "deg f": "f",
"degs f": "f", "degs f": "f",
"meter": "m", "meter": "m",
"meters": "m", "meters": "m",
"metre": "m", "metre": "m",
@ -223,7 +218,6 @@ _UNIT_ALIASES = {"mps" : "m s-1",
"foot": "ft", "foot": "ft",
"mile": "mi", "mile": "mi",
"miles": "mi" "miles": "mi"
} }
# A mapping of unit types to the avaible units # A mapping of unit types to the avaible units
@ -298,6 +292,7 @@ _TEMP_CONV_METHODS = {"c" : _k_to_c,
"f": _k_to_f "f": _k_to_f
} }
def dealias_and_clean_unit(unit): def dealias_and_clean_unit(unit):
"""Return the properly cleaned and dealiased unit name. """Return the properly cleaned and dealiased unit name.
@ -336,7 +331,7 @@ def check_units(unit, unit_type):
""" """
u_cleaned = dealias_and_clean_unit(unit) u_cleaned = dealias_and_clean_unit(unit)
if u_cleaned not in _VALID_UNITS[unit_type]: if u_cleaned not in _VALID_UNITS[unit_type]:
raise ValueError("invalid unit type '%s'" % unit) raise ValueError("invalid unit type '{}'".format(unit))
def do_conversion(var, vartype, var_unit, dest_unit): def do_conversion(var, vartype, var_unit, dest_unit):
@ -365,8 +360,3 @@ def do_conversion(var, vartype, var_unit, dest_unit):
return _apply_conv_fact(var, vartype, var_unit.lower(), u_cleaned) return _apply_conv_fact(var, vartype, var_unit.lower(), u_cleaned)
else: else:
return _apply_temp_conv(var, var_unit.lower(), u_cleaned) return _apply_temp_conv(var, var_unit.lower(), u_cleaned)

137
src/wrf/util.py

@ -214,7 +214,7 @@ def _generator_copy(gen):
res = module.get(funcname)(**argd) res = module.get(funcname)(**argd)
except AttributeError: except AttributeError:
res = getattr(module, funcname)(**argd) res = getattr(module, funcname)(**argd)
except: except Exception:
# This is the old way it used to work, but it looks like this was # This is the old way it used to work, but it looks like this was
# fixed by Python. # fixed by Python.
try: try:
@ -228,7 +228,7 @@ def _generator_copy(gen):
try: try:
argd = {key: argvals.locals[key] for key in argvals.args} argd = {key: argvals.locals[key] for key in argvals.args}
res = getattr(__main__, funcname)(**argd) res = getattr(__main__, funcname)(**argd)
except: except Exception:
# This was the old way it used to work, but appears to have # This was the old way it used to work, but appears to have
# been fixed by Python. # been fixed by Python.
res = getattr(__main__, funcname)(**argvals.locals) res = getattr(__main__, funcname)(**argvals.locals)
@ -669,10 +669,10 @@ def _corners_moved(wrfnc, first_ll_corner, first_ur_corner, latvar, lonvar):
end_idxs[-3] = i end_idxs[-3] = i
end_idxs = tuple(end_idxs) end_idxs = tuple(end_idxs)
if (first_ll_corner[0] != lats[start_idxs] or if (first_ll_corner[0] != lats[start_idxs]
first_ll_corner[1] != lons[start_idxs] or or first_ll_corner[1] != lons[start_idxs]
first_ur_corner[0] != lats[end_idxs] or or first_ur_corner[0] != lats[end_idxs]
first_ur_corner[1] != lons[end_idxs]): or first_ur_corner[1] != lons[end_idxs]):
return True return True
return False return False
@ -757,8 +757,14 @@ def is_moving_domain(wrfin, varname=None, latvar=either("XLAT", "XLAT_M"),
lon_coord = lonvar lon_coord = lonvar
lat_coord = latvar lat_coord = latvar
else: else:
lon_coord = coord_names[0] for name in coord_names:
lat_coord = coord_names[1] if name in _LAT_COORDS:
lat_coord = name
continue
if name in _LON_COORDS:
lon_coord = name
continue
else: else:
lon_coord = lonvar lon_coord = lonvar
lat_coord = latvar lat_coord = latvar
@ -907,9 +913,10 @@ def extract_dim(wrfin, dim):
# figure it out # figure it out
try: try:
s = wrfin.variables["P"].shape s = wrfin.variables["P"].shape
return s[-4] except Exception:
except:
raise ValueError("unsupported NetCDF reader") raise ValueError("unsupported NetCDF reader")
else:
return s[-4]
return d # PyNIO return d # PyNIO
@ -1038,12 +1045,10 @@ def _combine_dict(wrfdict, varname, timeidx, method, meta, _key):
# make it so that key_0 is leftmost # make it so that key_0 is leftmost
outdims = key_coordnames + list(first_array.dims[existing_cnt:]) outdims = key_coordnames + list(first_array.dims[existing_cnt:])
# Create the new 'key_n', value pairs # Create the new 'key_n', value pairs
for coordname, coordval in zip(key_coordnames, coord_vals): for coordname, coordval in zip(key_coordnames, coord_vals):
outcoords[coordname] = coordval outcoords[coordname] = coordval
outattrs = OrderedDict(first_array.attrs) outattrs = OrderedDict(first_array.attrs)
outarr = DataArray(outdata, name=outname, coords=outcoords, outarr = DataArray(outdata, name=outname, coords=outcoords,
@ -1181,15 +1186,22 @@ def _get_coord_names(wrfin, varname):
coord_names = coord_attr.split() coord_names = coord_attr.split()
else: else:
coord_names = coord_attr.decode().split() coord_names = coord_attr.decode().split()
lon_coord = coord_names[0]
lat_coord = coord_names[1]
try: for name in coord_names:
time_coord = coord_names[2] if name in _LAT_COORDS:
except IndexError: lat_coord = name
time_coord = None continue
else:
# Make sure they time variable wasn't removed if name in _LON_COORDS:
lon_coord = name
continue
if name in _TIME_COORD_VARS:
time_coord = name
continue
if time_coord is not None:
# Make sure the time variable wasn't removed
try: try:
_ = wrfnc.variables[time_coord] _ = wrfnc.variables[time_coord]
except KeyError: except KeyError:
@ -1279,7 +1291,6 @@ def _build_data_array(wrfnc, varname, timeidx, is_moving_domain, is_multifile,
except IndexError: except IndexError:
pass pass
coords = OrderedDict() coords = OrderedDict()
# Handle lat/lon coordinates and projection information if available # Handle lat/lon coordinates and projection information if available
@ -1605,8 +1616,6 @@ def _cat_files(wrfseq, varname, timeidx, is_moving, squeeze, meta, _key):
return _find_arr_for_time(wrfseq, varname, timeidx, is_moving, meta, return _find_arr_for_time(wrfseq, varname, timeidx, is_moving, meta,
_key) _key)
#time_idx_or_slice = timeidx if not multitime else slice(None)
# If all times are requested, need to build a new array and cat together # If all times are requested, need to build a new array and cat together
# all of the arrays in the sequence # all of the arrays in the sequence
wrf_iter = iter(wrfseq) wrf_iter = iter(wrfseq)
@ -1652,7 +1661,8 @@ def _cat_files(wrfseq, varname, timeidx, is_moving, squeeze, meta, _key):
outxtimes = get_cached_item(_key, timekey) outxtimes = get_cached_item(_key, timekey)
if outxtimes is None: if outxtimes is None:
outxtimes = np.empty(outdims[0]) outxtimes = np.empty(outdims[0])
outxtimes[startidx:endidx] = to_np(first_var.coords[timename][:]) outxtimes[startidx:endidx] = to_np(
first_var.coords[timename][:])
else: else:
timecached = True timecached = True
@ -1664,7 +1674,8 @@ def _cat_files(wrfseq, varname, timeidx, is_moving, squeeze, meta, _key):
outlats = get_cached_item(_key, latkey) outlats = get_cached_item(_key, latkey)
if outlats is None: if outlats is None:
outlats = np.empty(outcoorddims, first_var.dtype) outlats = np.empty(outcoorddims, first_var.dtype)
outlats[startidx:endidx, :] = to_np(first_var.coords[latname][:]) outlats[startidx:endidx, :] = to_np(
first_var.coords[latname][:])
else: else:
latcached = True latcached = True
@ -1672,11 +1683,11 @@ def _cat_files(wrfseq, varname, timeidx, is_moving, squeeze, meta, _key):
outlons = get_cached_item(_key, lonkey) outlons = get_cached_item(_key, lonkey)
if outlons is None: if outlons is None:
outlons = np.empty(outcoorddims, first_var.dtype) outlons = np.empty(outcoorddims, first_var.dtype)
outlons[startidx:endidx, :] = to_np(first_var.coords[lonname][:]) outlons[startidx:endidx, :] = to_np(
first_var.coords[lonname][:])
else: else:
loncached = True loncached = True
startidx = endidx startidx = endidx
while True: while True:
try: try:
@ -1951,8 +1962,8 @@ def _join_files(wrfseq, varname, timeidx, is_moving, meta, _key):
# For join, the times are a function of fileidx # For join, the times are a function of fileidx
file_times = extract_times(wrfnc, ALL_TIMES, meta=False, file_times = extract_times(wrfnc, ALL_TIMES, meta=False,
do_xtime=False) do_xtime=False)
time_coord[file_idx, 0:numtimes] = np.asarray(file_times, time_coord[file_idx, 0:numtimes] = np.asarray(
"datetime64[ns]")[:] file_times, "datetime64[ns]")[:]
if timename is not None and not timecached: if timename is not None and not timecached:
xtimedata = wrfnc.variables[timename][:] xtimedata = wrfnc.variables[timename][:]
@ -2501,10 +2512,8 @@ def extract_times(wrfin, timeidx, method="cat", squeeze=True, cache=None,
outname = "XTIME" outname = "XTIME"
outarr = DataArray(time_arr, name=outname, coords=outcoords, outarr = DataArray(time_arr, name=outname, coords=outcoords,
dims=outdimnames, attrs=outattrs) dims=outdimnames, attrs=outattrs)
else: else:
outarr = time_arr outarr = time_arr
@ -2657,7 +2666,7 @@ def get_right_slices(var, right_ndims, fixed_val=0):
[slice(None)]*right_ndims) [slice(None)]*right_ndims)
def get_proj_params(wrfin):#, timeidx=0, varname=None): def get_proj_params(wrfin):
"""Return a tuple of latitude, longitude, and projection parameters from """Return a tuple of latitude, longitude, and projection parameters from
a WRF output file object or a sequence of WRF output file objects. a WRF output file object or a sequence of WRF output file objects.
@ -2684,7 +2693,8 @@ def get_proj_params(wrfin):#, timeidx=0, varname=None):
longitude coordinate, and global projection attributes. longitude coordinate, and global projection attributes.
""" """
proj_params = extract_global_attrs(wrfin, attrs=("MAP_PROJ", proj_params = extract_global_attrs(wrfin,
attrs=("MAP_PROJ",
"CEN_LAT", "CEN_LON", "CEN_LAT", "CEN_LON",
"TRUELAT1", "TRUELAT2", "TRUELAT1", "TRUELAT2",
"MOAD_CEN_LAT", "STAND_LON", "MOAD_CEN_LAT", "STAND_LON",
@ -2697,7 +2707,7 @@ def get_proj_params(wrfin):#, timeidx=0, varname=None):
def from_args(func, argnames, *args, **kwargs): def from_args(func, argnames, *args, **kwargs):
"""Return a mapping of argument name to value for the called function. """Return a mapping of argument name to value for the called function.
This function parses the function \*args and \*\*kwargs to obtain the \ This function parses the function args and kwargs to obtain the
desired argument value. If the argument has not been passed in, the value desired argument value. If the argument has not been passed in, the value
is taken from the default keyword argument value. is taken from the default keyword argument value.
@ -2706,7 +2716,7 @@ def from_args(func, argnames, *args, **kwargs):
Note: Note:
This function currently does not work with functions that contain This function currently does not work with functions that contain
\*args or \*\*kwargs arguments. variable length args or kwargs arguments.
Args: Args:
@ -2750,7 +2760,7 @@ def _args_to_list2(func, args, kwargs):
Note: Note:
This function currently does not work with functions that contain This function currently does not work with functions that contain
*args or **kwargs arguments. variable length args or kwargs arguments.
Args: Args:
@ -2837,7 +2847,7 @@ def _args_to_list3(func, args, kwargs):
Note: Note:
This function currently does not work with functions that contain This function currently does not work with functions that contain
*args or **kwargs arguments. variable length args or kwargs arguments.
Args: Args:
@ -2872,7 +2882,7 @@ def args_to_list(func, args, kwargs):
Note: Note:
This function currently does not work with functions that contain This function currently does not work with functions that contain
\*args or \*\*kwargs arguments. variable length args or kwargs arguments.
Args: Args:
@ -3026,13 +3036,27 @@ def psafilepath():
def get_filepath(obj): def get_filepath(obj):
"""Return the file path for the specified object.
This is used to return the file path for a netcdf object. If the
particular object does not have the appropriate file path information,
then one is created based on the timestep in the file.
Args:
obj: An object.
Returns:
:obj:`str`: A string for a file path.
"""
try: try:
path = obj.filepath() path = obj.filepath()
except AttributeError: except AttributeError:
try: try:
path = obj.file.path path = obj.file.path
except: except AttributeError:
# Let's make up a filename from the first file time # Let's make up a filename from the first file time
found = False found = False
times = extract_times(obj, None, meta=False, do_xtime=False) times = extract_times(obj, None, meta=False, do_xtime=False)
@ -3046,6 +3070,7 @@ def get_filepath(obj):
return path return path
def get_id(obj, prefix=''): def get_id(obj, prefix=''):
"""Return the cache id. """Return the cache id.
@ -3221,6 +3246,7 @@ def geo_bounds(var=None, wrfin=None, varname=None, timeidx=0, method="cat",
# Non-moving domains # Non-moving domains
return GeoBounds(lats=lats, lons=lons) return GeoBounds(lats=lats, lons=lons)
def _get_wrf_proj_geobnds(var, wrfin, varname, timeidx, method, squeeze, def _get_wrf_proj_geobnds(var, wrfin, varname, timeidx, method, squeeze,
cache): cache):
"""Return the :class:`wrf.WrfProj` subclass and :class:`wrf.GeoBounds`. """Return the :class:`wrf.WrfProj` subclass and :class:`wrf.GeoBounds`.
@ -3724,7 +3750,8 @@ def cartopy_xlim(var=None, geobounds=None, wrfin=None, varname=None, timeidx=0,
""" """
wrf_proj, native_geobnds = _get_wrf_proj_geobnds(var, wrfin, varname, wrf_proj, native_geobnds = _get_wrf_proj_geobnds(var, wrfin, varname,
timeidx, method, squeeze, cache) timeidx, method, squeeze,
cache)
if geobounds is not None: if geobounds is not None:
return wrf_proj.cartopy_xlim(geobounds) return wrf_proj.cartopy_xlim(geobounds)
@ -3810,7 +3837,8 @@ def cartopy_ylim(var=None, geobounds=None, wrfin=None, varname=None, timeidx=0,
""" """
wrf_proj, native_geobnds = _get_wrf_proj_geobnds(var, wrfin, varname, wrf_proj, native_geobnds = _get_wrf_proj_geobnds(var, wrfin, varname,
timeidx, method, squeeze, cache) timeidx, method, squeeze,
cache)
if geobounds is not None: if geobounds is not None:
return wrf_proj.cartopy_ylim(geobounds) return wrf_proj.cartopy_ylim(geobounds)
@ -3897,30 +3925,3 @@ def is_latlon_pair(pair):
return (pair.lat is not None and pair.lon is not None) return (pair.lat is not None and pair.lon is not None)
else: else:
return False return False

1
src/wrf/version.py

@ -1,2 +1 @@
__version__ = "1.3.1" __version__ = "1.3.1"

Loading…
Cancel
Save