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

@@ -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)

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

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")

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
 

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
     

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