Add Visualize.plot_bias!

This functionality is already in ClimaCoupler and should belong in
ClimaAnalysis. This commit ports over plot_bias! to ClimaAnalysis.

NEWS.md

@@ -253,6 +253,29 @@ julia> units(se_var)
 "K^2"
 ```
 
+### Plotting bias
+
+Building upon the other features introduced in this release, you can now directly plot bias
+and root mean squared error between two variables with the `plot_bias_on_globe!` function.
+Typically, this is done to compare simulated data against observations.
+
+In the example below, we plot the bias between our simulation and some observations stored
+in `ta_1d_average.nc`.
+
+```julia
+import ClimaAnalysis
+import ClimaAnalysis.Visualize: plot_bias_on_globe!
+import GeoMakie
+import CairoMakie
+
+obs_var = ClimaAnalysis.OutputVar("ta_1d_average.nc")
+sim_var = ClimaAnalysis.get(ClimaAnalysis.simdir("simulation_output"), "ta")
+
+fig = CairoMakie.Figure()
+plot_bias_on_globe!(fig, sim_var, obs_var)
+CairoMakie.save("myfigure.pdf", fig)
+```
+
 ## Bug fixes
 
 - Increased the default value for `warp_string` to 72.

docs/src/api.md

@@ -106,4 +106,5 @@ Visualize.oceanmask
 Visualize.landmask
 Visualize.contour2D_on_globe!
 Visualize.heatmap2D_on_globe!
+Visualize.plot_bias_on_globe!
 ```

docs/src/visualize.md

@@ -45,3 +45,29 @@ more easily.
 The output might look something like:
 
 ![oceanmask](./assets/oceanmask.png)
+
+### Plotting bias
+
+After [computing the bias](@ref bias) between observational and simulation data, you may
+want to plot the bias and display information such as the root mean squared error (RMSE) and
+the global bias in the plot. To do this, you use the function [`plot_bias_on_globe!(fig, sim,
+obs)`](@ref Visualize.plot_bias_on_globe!). In the example below, we plot the bias between our
+simulation and some observations stored in `ta_1d_average.nc`.
+
+```julia
+import ClimaAnalysis
+import ClimaAnalysis.Visualize: plot_bias_on_globe!
+import GeoMakie
+import CairoMakie
+
+obs_var = ClimaAnalysis.OutputVar("ta_1d_average.nc")
+sim_var = ClimaAnalysis.get(ClimaAnalysis.simdir("simulation_output"), "ta")
+
+fig = CairoMakie.Figure()
+plot_bias_on_globe!(fig, sim_var, obs_var)
+CairoMakie.save("myfigure.pdf", fig)
+```
+
+The output produces something like:
+
+![biasplot](./assets/bias_plot.png)

ext/ClimaAnalysisGeoMakieExt.jl

@@ -254,4 +254,125 @@ function Visualize.contour2D_on_globe!(
     )
 end
 
+"""
+    plot_bias_on_globe!(fig::Makie.Figure,
+                        sim::ClimaAnalysis.OutputVar,
+                        obs::ClimaAnalysis.OutputVar;
+                        cmap_extrema = extrema(ClimaAnalysis.bias(sim, obs).data),
+                        p_loc = (1, 1),
+                        plot_coastline = true,
+                        plot_colorbar = true,
+                        mask = nothing,
+                        more_kwargs)
+    plot_bias_on_globe!(grid_layout::Makie.GridLayout,
+                        sim::ClimaAnalysis.OutputVar,
+                        obs::ClimaAnalysis.OutputVar;
+                        cmap_extrema = extrema(ClimaAnalysis.bias(sim, obs).data),
+                        p_loc = (1, 1),
+                        plot_coastline = true,
+                        plot_colorbar = true,
+                        mask = nothing,
+                        more_kwargs)
+
+
+Plot the bias (`sim.data - var.data`) on a projected geoid. The gloal bias and root mean
+squared error (RMSE) are computed and can be found in the title of the plot. This function
+plots the returned `OutputVar` of `ClimaAnalysis.bias(sim, obs)`. See also
+[`ClimaAnalysis.bias`](@ref).
+
+The plot comes with labels, units, and a colorbar. This function uses a constrained colormap
+based on the values of `cmap_extrema`.
+
+The dimensions have to be longitude and latitude.
+
+`mask` has to be an object that can be plotted by `Makie.poly`. `ClimaAnalysis` comes with
+predefined masks, check out [`Visualize.oceanmask`](@ref) and [`Visualize.landmask`](@ref).
+
+!!! note
+    Masking does not affect the colorbar. If you have values defined beneath the map, they
+    can still affect the colorbar.
+
+Additional arguments to the plotting and axis functions
+=======================================================
+
+`more_kwargs` can be a dictionary that maps symbols to additional options for:
+- the axis (`:axis`)
+- the plotting function (`:plot`)
+- the colorbar (`:cb`)
+- the coastline (`:coast`)
+- the mask (`:mask`)
+
+The coastline is plotted from `GeoMakie.coastline` using the `lines!` plotting function.
+
+The values are splatted in the relevant functions. Populate them with a
+Dictionary of `Symbol`s => values to pass additional options.
+"""
+function Visualize.plot_bias_on_globe!(
+    place::MakiePlace,
+    sim::ClimaAnalysis.OutputVar,
+    obs::ClimaAnalysis.OutputVar;
+    cmap_extrema = extrema(ClimaAnalysis.bias(sim, obs).data),
+    p_loc = (1, 1),
+    plot_coastline = true,
+    plot_colorbar = true,
+    mask = nothing,
+    more_kwargs = Dict(
+        :plot => Dict(),
+        :cb => Dict(),
+        :axis => Dict(),
+        :coast => Dict(:color => :black),
+        :mask => Dict(),
+    ),
+)
+    bias_var = ClimaAnalysis.bias(sim, obs)
+    global_bias = round(bias_var.attributes["global_bias"], sigdigits = 3)
+    rmse = round(ClimaAnalysis.global_rmse(sim, obs), sigdigits = 3)
+    units = ClimaAnalysis.units(bias_var)
+
+    bias_var.attributes["long_name"] *= " (RMSE: $rmse $units, Global bias: $global_bias $units)"
+    min_level, max_level = cmap_extrema
+
+    # Make sure that 0 is at the center
+    cmap = Visualize._constrained_cmap(
+        Makie.cgrad(:vik).colors,
+        min_level,
+        max_level;
+        categorical = true,
+    )
+    nlevels = 11
+    # Offset so that it covers 0
+    levels = collect(range(min_level, max_level, length = nlevels))
+    offset = levels[argmin(abs.(levels))]
+    levels = levels .- offset
+    ticklabels = map(x -> string(round(x; digits = 0)), levels)
+    ticks = (levels, ticklabels)
+
+    default_kwargs = Dict(
+        :plot => Dict(
+            :colormap => cmap,
+            :levels => levels,
+            :extendhigh => :auto,
+            :extendlow => :auto,
+        ),
+        :cb => Dict(:ticks => ticks),
+    )
+
+    # Function for recursively merging two dictionaries if the values of the dictionaries
+    # are dictionaries and the values of those are also dictionaries and so on
+    # See: https://discourse.julialang.org/t/multi-layer-dict-merge/27261/6
+    recursive_merge(x::AbstractDict...) = merge(recursive_merge, x...)
+    recursive_merge(x...) = x[end]
+    default_and_more_kwargs = recursive_merge(default_kwargs, more_kwargs)
+
+    return Visualize.contour2D_on_globe!(
+        place,
+        bias_var;
+        p_loc = p_loc,
+        plot_coastline = plot_coastline,
+        plot_colorbar = plot_colorbar,
+        mask = mask,
+        more_kwargs = default_and_more_kwargs,
+    )
+end
+
 end

src/Visualize.jl

@@ -28,4 +28,6 @@ function contour2D_on_globe! end
 
 function heatmap2D_on_globe! end
 
+function plot_bias_on_globe! end
+
 end

test/test_GeoMakieExt.jl

@@ -100,4 +100,40 @@ using OrderedCollections
 
     output_name = joinpath(tmp_dir, "test_contours2D_globe_with_oceanmask.png")
     Makie.save(output_name, fig5)
+
+    # Test plot_bias
+    fig6 = Makie.Figure()
+
+    lon = collect(range(-179.5, 179.5, 360))
+    lat = collect(range(-89.5, 89.5, 180))
+    data = collect(reshape(-32400:32399, (360, 180))) ./ (32399.0 / 5.0)
+    dims = OrderedDict(["lon" => lon, "lat" => lat])
+    attribs =
+        Dict("long_name" => "idk", "short_name" => "short", "units" => "kg")
+    dim_attribs = OrderedDict([
+        "lon" => Dict("units" => "deg"),
+        "lat" => Dict("units" => "deg"),
+    ])
+    var = ClimaAnalysis.OutputVar(attribs, dims, dim_attribs, data)
+
+    data_zero = zeros(length(lon), length(lat))
+    var_zero = ClimaAnalysis.OutputVar(attribs, dims, dim_attribs, data_zero)
+
+    ClimaAnalysis.Visualize.plot_bias_on_globe!(fig6, var, var_zero)
+    output_name = joinpath(tmp_dir, "plot_bias.png")
+    Makie.save(output_name, fig6)
+
+    # Test plot bias with keyword arguments
+    fig7 = Makie.Figure()
+    ClimaAnalysis.Visualize.plot_bias_on_globe!(
+        fig7,
+        var,
+        var_zero,
+        more_kwargs = Dict(
+            :axis => Dict(:title => "no title"),
+            :plot => Dict(:extendhigh => nothing),
+        ),
+    )
+    output_name = joinpath(tmp_dir, "plot_bias_kwargs.png")
+    Makie.save(output_name, fig7)
 end