SmoothOperator (#64)

* Fix citation typo for rivers

* Add smoothOperator

Project.toml

@@ -1,7 +1,7 @@
 name = "AIBECS"
 uuid = "ace601d6-714c-11e9-04e5-89b7fad23838"
 authors = ["Benoit Pasquier <[email protected]>"]
-version = "0.8.5"
+version = "0.8.6"
 
 [deps]
 BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
@@ -15,6 +15,7 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 FieldMetadata = "bf96fef3-21d2-5d20-8afa-0e7d4c32a885"
 Flatten = "4c728ea3-d9ee-5c9a-9642-b6f7d7dc04fa"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MetadataArrays = "49441bc9-da82-574f-b07c-a0d10dd4ac13"
@@ -42,6 +43,7 @@ Distributions = "0.23, 0.24"
 FieldMetadata = "0.3"
 Flatten = "0.3, 0.4"
 ForwardDiff = "0.10"
+ImageFiltering = "0.6"
 Interpolations = "0.12, 0.13"
 MetadataArrays = "0.1"
 NCDatasets = "0.10, 0.11"

docs/lit/tutorials/5_river_discharge.jl

@@ -123,3 +123,32 @@ p = RadioRiversParameters(τ = 50.0u"yr")
 prob = SteadyStateProblem(F, ∇ₓF, x, p)
 s_τ50 = solve(prob, CTKAlg()).u * u"mol/m^3"
 plothorizontalslice(s_τ50, grd, zunit=u"μmol/m^3", depth=500, color=cmap, clim=(0,1))
+
+# Point-wise sources like the rivers used here can be problematic numerically because these
+# can generate large gradients and numerical noise with the given spatial discretization:
+
+cmap = :RdYlBu_4
+plothorizontalslice(s_τ50, grd, zunit=u"μmol/m^3", depth=0, color=cmap, clim=(-10,10))
+
+# Note, for an example of numerical noise, the negative values appearing in red.
+# Hence, it might be wise to smooth the 3D field of the river sources.
+# We can do this using the `smooth_operator` function
+# (`smooth_operator` creates an operator that can be applied to smooth the AIBECS vector
+# along the stencil of the given transport matrix and in a volume-conservative way).
+# And we can check that the negative values disappear after smoothing of the sources:
+
+S = smooth_operator(grd, T_OCIM2)
+s_0 = S * (S * s_0) # smooth river sources twice
+s_smooth = solve(prob, CTKAlg()).u * u"mol/m^3"
+plothorizontalslice(s_smooth, grd, zunit=u"μmol/m^3", depth=0, color=cmap, clim=(-10,10))
+
+# Note that such numerical noise generally dampens out with distance and depth
+
+cmap = :viridis
+plot(
+    plothorizontalslice(s_τ50, grd, zunit=u"μmol/m^3", depth=200, color=cmap, clim=(0,2)),
+    plothorizontalslice(s_smooth, grd, zunit=u"μmol/m^3", depth=200, color=cmap, clim=(0,2)),
+    layout=(2,1)
+)
+
+# Nevertheless, it is always good to reduce noise as much as possible!

src/AIBECS.jl

@@ -25,6 +25,7 @@ using NearestNeighbors
 @reexport using MetadataArrays
 using Shapefile
 using Bijectors
+using ImageFiltering
 
 
 

src/Rivers.jl

@@ -241,7 +241,7 @@ citation() = """
 - For the formal publication that describe the dataset: Dai, A., and K. E. Trenberth, 2002: Estimates of freshwater discharge from continents: Latitudinal and seasonal variations. J. Hydrometeorol., 3, 660–687.
 - For descriptions of updates to the dataset:
     - Dai, A., T. Qian, K. E. Trenberth, and J. D Milliman, 2009: Changes in continental freshwater discharge from 1948-2004. J. Climate, 22, 2773–2791.
-    - Dai, A., 2016: Historical and future changes in streamflow and continental runoff: A review, in Terrestrial Water Cycle and Climate Change: Natural and Human-Induced Impacts, Geophysical Monograph 221. Ed. Qiuhong Tang and Taikan Oki, AGU, John Wiley & Sons, 17–37 (DOI: 10.1002/9781118971772).
+    - Dai, A., 2016: Historical and future changes in streamflow and continental runoff: A review, in Terrestrial Water Cycle and Climate Change: Natural and Human-Induced Impacts, Geophysical Monograph 221. Ed. Qiuhong Tang and Taikan Oki, AGU, John Wiley & Sons, 17–37 (DOI: 10.1002/9781118971772).
 """
 
 uconvert(u, r::River) = River(r.name, r.lon, r.lat, uconvert(u, r.VFR))

src/diagnostics.jl

@@ -85,7 +85,25 @@ export directional_transport, directional_transports
 
 
 
+"""
+    smooth_operator(grd, T; σs=(1.0, 1.0, 0.25))
 
+return a matrix of the same size and sparsity as `T` that smoothes data using
+a Gaussian Kernel for values, but conserving mass.
+"""
+function smooth_operator(grd, T; σs=(1.0, 1.0, 0.25))
+    st = stencil(grd, T)
+    weight(dir) = Kernel.gaussian(σs)[dir]
+    Isym, Jsym, _, _ = symmetric_IJVVᵀ(T)
+    dirs = directions(Isym, Jsym, grd)
+    A = sparse(Isym, Jsym, weight.(dirs))
+    # We add mass conservation to the Gaussian smoothing matrix
+    # by modifying A column by column such that A conserves mass,
+    # i.e., such that vᵀ A x = vᵀ x
+    v = depthvec(grd)
+    A * sparse(Diagonal(v ./ (A' * v)))
+end
+export smooth_operator