diff --git a/Project.toml b/Project.toml index 53be35b5..0fc0c5f5 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "AIBECS" uuid = "ace601d6-714c-11e9-04e5-89b7fad23838" authors = ["Benoit Pasquier "] -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" diff --git a/docs/lit/tutorials/5_river_discharge.jl b/docs/lit/tutorials/5_river_discharge.jl index a8605e8e..d83e245f 100644 --- a/docs/lit/tutorials/5_river_discharge.jl +++ b/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! diff --git a/src/AIBECS.jl b/src/AIBECS.jl index 47d182c4..5669230e 100644 --- a/src/AIBECS.jl +++ b/src/AIBECS.jl @@ -25,6 +25,7 @@ using NearestNeighbors @reexport using MetadataArrays using Shapefile using Bijectors +using ImageFiltering diff --git a/src/Rivers.jl b/src/Rivers.jl index 56de7aa7..fbac7f71 100644 --- a/src/Rivers.jl +++ b/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)) diff --git a/src/diagnostics.jl b/src/diagnostics.jl index ba929597..1b28d3eb 100644 --- a/src/diagnostics.jl +++ b/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