From d67586c09b57c03a044c960c9270831c59e5b7af Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Thu, 11 Oct 2018 14:51:46 -0600 Subject: [PATCH] Fixed the smoothing algorithm. It now works like a typical smoothing kernel using convolution. Users can now specify the center weight. The vintrp algorithm needed to be modified to use the new function signature. Updated the documentation to give a better description of how it works. Fixes #67. --- doc/source/conf.py | 7 +++++-- fortran/wrf_user.f90 | 37 +++++++++++++++++++++---------------- src/wrf/computation.py | 27 ++++++++++++++++++++++++--- src/wrf/extension.py | 7 ++++--- src/wrf/interp.py | 2 +- src/wrf/metadecorators.py | 4 +++- 6 files changed, 58 insertions(+), 26 deletions(-) diff --git a/doc/source/conf.py b/doc/source/conf.py index b248434..7d339dc 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -75,15 +75,18 @@ extensions = [ 'sphinx.ext.napoleon', 'sphinx.ext.autosummary', 'sphinx.ext.intersphinx', + 'sphinx.ext.mathjax' ] +#mathjax_path = "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/MathJax.js?config=TeX-MML-AM_CHTML" + intersphinx_mapping = { 'matplotlib': ('http://matplotlib.org/', None), 'basemap' : ('http://matplotlib.org/basemap/', None), - 'cartopy' : ( - 'http://scitools.org.uk/cartopy/docs/latest/index.html', None), 'python': ('http://docs.python.org/3/', None), 'numpy': ('http://docs.scipy.org/doc/numpy/', None), + 'scipy' : + ('https://docs.scipy.org/doc/scipy/reference/', None), 'xarray': ('http://xarray.pydata.org/en/stable/', None), 'wrapt': ('http://wrapt.readthedocs.io/en/latest/', None) diff --git a/fortran/wrf_user.f90 b/fortran/wrf_user.f90 index e0bfac8..689ac6b 100644 --- a/fortran/wrf_user.f90 +++ b/fortran/wrf_user.f90 @@ -479,7 +479,7 @@ END SUBROUTINE DCOMPUTESEAPRS ! must make the same change below to filter2d. ! NCLFORTSTART -SUBROUTINE DFILTER2D(a, b, nx, ny, it, missing) +SUBROUTINE DFILTER2D(a, b, nx, ny, it, missing, cenweight) IMPLICIT NONE @@ -488,14 +488,16 @@ SUBROUTINE DFILTER2D(a, b, nx, ny, it, missing) INTEGER, INTENT(IN) :: nx, ny, it REAL(KIND=8), DIMENSION(nx, ny), INTENT(INOUT) :: a - REAL(KIND=8), INTENT(IN) :: missing + REAL(KIND=8), INTENT(IN) :: missing, cenweight REAL(KIND=8), DIMENSION(nx, ny), INTENT(INOUT) :: b ! NCLEND - REAL(KIND=8), PARAMETER :: COEF=0.25D0 - INTEGER :: i, j, iter + REAL(KIND=8) :: cenmult, coef + + cenmult = (cenweight) / 2. + coef = 1.0 / (4. + cenweight) DO iter=1,it !$OMP PARALLEL DO COLLAPSE(2) SCHEDULE(runtime) @@ -508,25 +510,25 @@ SUBROUTINE DFILTER2D(a, b, nx, ny, it, missing) !$OMP PARALLEL DO COLLAPSE(2) SCHEDULE(runtime) DO j=2,ny-1 - DO i=1,nx + DO i=2,nx-1 IF (b(i,j-1) .EQ. missing .OR. b(i,j) .EQ. missing .OR. & b(i,j+1) .EQ. missing) THEN a(i,j) = a(i,j) ELSE - a(i,j) = a(i,j) + COEF*(b(i,j-1) - 2*b(i,j) + b(i,j+1)) + a(i,j) = coef*(b(i,j-1) + cenmult*b(i,j) + b(i,j+1)) END IF END DO END DO !$OMP END PARALLEL DO !$OMP PARALLEL DO COLLAPSE(2) SCHEDULE(runtime) - DO j=1,ny + DO j=2,ny-1 DO i=2,nx-1 IF (b(i-1,j) .EQ. missing .OR. b(i,j) .EQ. missing .OR. & b(i+1,j) .EQ. missing) THEN a(i,j) = a(i,j) ELSE - a(i,j) = a(i,j) + COEF*(b(i-1,j) - 2*b(i,j) + b(i+1,j)) + a(i,j) = a(i,j) + coef*(b(i-1,j) + cenmult*b(i,j) + b(i+1,j)) END IF END DO END DO @@ -556,7 +558,7 @@ END SUBROUTINE DFILTER2D ! must make the same change below to dfilter2d. ! NCLFORTSTART -SUBROUTINE FILTER2D(a, b, nx, ny, it, missing) +SUBROUTINE FILTER2D(a, b, nx, ny, it, missing, cenweight) IMPLICIT NONE !f2py threadsafe @@ -564,14 +566,17 @@ SUBROUTINE FILTER2D(a, b, nx, ny, it, missing) INTEGER, INTENT(IN) :: nx, ny, it REAL(KIND=4), DIMENSION(nx, ny), INTENT(INOUT) :: a - REAL(KIND=4), INTENT(IN) :: missing + REAL(KIND=4), INTENT(IN) :: missing, cenweight REAL(KIND=4), DIMENSION(nx, ny), INTENT(INOUT) :: b ! NCLEND - REAL(KIND=4), PARAMETER :: COEF=0.25 - + !REAL(KIND=8), PARAMETER :: COEF = .125 INTEGER :: i, j, iter + REAL(KIND=8) :: cenmult, coef + + cenmult = (cenweight) / 2. + coef = 1.0 / (4. + cenweight) !$OMP PARALLEL @@ -586,25 +591,25 @@ SUBROUTINE FILTER2D(a, b, nx, ny, it, missing) !$OMP DO COLLAPSE(2) SCHEDULE(runtime) DO j=2,ny-1 - DO i=1,nx + DO i=2,nx-1 IF (b(i,j-1) .EQ. missing .OR. b(i,j) .EQ. missing .OR. & b(i,j+1) .EQ. missing) THEN a(i,j) = a(i,j) ELSE - a(i,j) = a(i,j) + COEF*(b(i,j-1) - 2*b(i,j) + b(i,j+1)) + a(i,j) = coef*(b(i,j-1) + cenmult*b(i,j) + b(i,j+1)) END IF END DO END DO !$OMP END DO !$OMP DO COLLAPSE(2) SCHEDULE(runtime) - DO j=1,ny + DO j=2,ny-1 DO i=2,nx-1 IF (b(i-1,j) .EQ. missing .OR. b(i,j) .EQ. missing .OR. & b(i+1,j) .EQ. missing) THEN a(i,j) = a(i,j) ELSE - a(i,j) = a(i,j) + COEF*(b(i-1,j) - 2*b(i,j)+b(i+1,j)) + a(i,j) = a(i,j) + coef*(b(i-1,j) + cenmult*b(i,j) + b(i+1,j)) END IF END DO END DO diff --git a/src/wrf/computation.py b/src/wrf/computation.py index 68486b2..68ee1cf 100644 --- a/src/wrf/computation.py +++ b/src/wrf/computation.py @@ -677,10 +677,24 @@ def uvmet(u, v, lat, lon, cen_long, cone, meta=True, units="m s-1"): @set_smooth_metdata() -def smooth2d(field, passes, meta=True): +def smooth2d(field, passes, cenweight=2.0, meta=True): """Return the field smoothed. - This routine does not modify the original data. + The smoothing kernel applied is: + + .. math:: + + \\frac{1}{4 + cenweight} * \\begin{bmatrix} + 0 & 1 & 0 \\\\ + 1 & cenweight & 1 \\\\ + 0 & 1 & 0 + \\end{bmatrix} + + Data values along the borders are left unchanged. This routine does not + modify the original data supplied by the *field* parameter.. + + If you need more general purpose multidimensional filtering tools, + try the :meth:`scipy.ndimage.convolve` method. Args: @@ -692,6 +706,9 @@ def smooth2d(field, passes, meta=True): passes (:obj:`int`): The number of smoothing passes. + cenweight (:obj:`float`, optional): The weight to apply to the + center of the smoothing kernel. Default is 2.0. + meta (:obj:`bool`): Set to False to disable metadata and return :class:`numpy.ndarray` instead of :class:`xarray.DataArray`. Default is True. @@ -705,9 +722,13 @@ def smooth2d(field, passes, meta=True): be either a :class:`numpy.ndarray` or a :class:`numpy.ma.MaskedArray` depending on the type for *field*. + See Also: + + :meth:`scipy.ndimage.convolve` + """ - return _smooth2d(field, passes) + return _smooth2d(field, passes, cenweight) @set_cape_alg_metadata(is2d=True, copyarg="pres_hpa") diff --git a/src/wrf/extension.py b/src/wrf/extension.py index d7636c5..0c4355d 100755 --- a/src/wrf/extension.py +++ b/src/wrf/extension.py @@ -799,10 +799,10 @@ def _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, fill_nocloud, @check_args(0, 2, (2,)) -@left_iteration(2, 2, ref_var_idx=0, ignore_args=(1,)) +@left_iteration(2, 2, ref_var_idx=0, ignore_args=(1,2)) @cast_type(arg_idxs=(0,)) @extract_and_transpose() -def _smooth2d(field, passes, outview=None): +def _smooth2d(field, passes, cenweight, outview=None): """Wrapper for dfilter2d. Located in wrf_user.f90. @@ -826,7 +826,8 @@ def _smooth2d(field, passes, outview=None): dfilter2d(outview, field_tmp, passes, - missing) + missing, + cenweight) return outview diff --git a/src/wrf/interp.py b/src/wrf/interp.py index 978f854..ec6e9d2 100755 --- a/src/wrf/interp.py +++ b/src/wrf/interp.py @@ -656,7 +656,7 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False, method=method, squeeze=squeeze, cache=cache, meta=False, _key=_key) - smsfp = _smooth2d(sfp, 3) + smsfp = _smooth2d(sfp, 3, 2.0) vcor = 0 diff --git a/src/wrf/metadecorators.py b/src/wrf/metadecorators.py index 573c7fc..63e1195 100644 --- a/src/wrf/metadecorators.py +++ b/src/wrf/metadecorators.py @@ -1839,11 +1839,12 @@ def set_smooth_metdata(): if not xarray_enabled() or not do_meta: return wrapped(*args, **kwargs) - argvars = from_args(wrapped, ("field", "passes"), + argvars = from_args(wrapped, ("field", "passes", "cenweight"), *args, **kwargs) field = argvars["field"] passes = argvars["passes"] + cenweight = argvars["cenweight"] result = wrapped(*args, **kwargs) @@ -1858,6 +1859,7 @@ def set_smooth_metdata(): outcoords.update(field.coords) outattrs.update(field.attrs) outattrs["passes"] = passes + outattrs["cenweight"] = cenweight if isinstance(result, ma.MaskedArray): outattrs["_FillValue"] = result.fill_value