Update functions computing bias/RMSE to use mask

All functions that compute bias or RMSE now have a mask parameter which
take in `apply_landmask` or `apply_oceanmask`.

Author: ph-kev

NEWS.md

@@ -1,5 +1,25 @@
 ClimaAnalysis.jl Release Notes
 ===============================
+v0.5.11
+-------
+
+## Features
+
+### Applying a land/sea mask to `OutputVar`s
+
+There is now support to applying a land or sea mask to `OutputVar`s though
+`apply_landmask(var)` and `apply_oceanmask(var)` respectively. A land mask sets all the
+coordinates corresponding to land to zeros in the data of the `OutputVar` and a sea mask
+sets all the coordinates corresponding to ocean to zeros in the data of the `OutputVar`.
+Furthermore, the parameter `mask` is added to the functions `bias`, `global_bias`,
+`squared_error`,`global_mse`, and `global_rmse` which takes either `apply_landmask` and
+`apply_oceanmask`. This is useful for computing these quantities only over the land or sea.
+
+```julia
+var_no_land = ClimaAnalysis.apply_landmask(var)
+var_no_ocean = ClimaAnalysis.apply_oceanmask(var)
+```
+
 v0.5.10
 -------
 

docs/src/var.md

@@ -226,7 +226,7 @@ longitude and latitude. Furthermore, units must be supplied for data and dimensi
 and `obs` and the units for longitude and latitude should be degrees.
 
 Consider the following example, where we compute the bias and RMSE between our simulation
-and some observations stored in "ta\_1d\_average.nc".
+and some observations stored in `ta_1d_average.nc`.
 
 ```@julia bias_and_mse
 julia> obs_var = OutputVar("ta_1d_average.nc"); # load in observational data
@@ -255,3 +255,19 @@ julia> global_rmse(sim, obs)
 julia> units(se_var)
 "K^2"
 ```
+
+### Masking
+Bias and squared error can be computed only over the land or ocean through the `mask` parameter.
+As of now, the mask parameter takes in `apply_oceanmask` or `apply_oceanmask`. See the
+example below of this usage.
+
+```julia
+# Do not consider the ocean when computing the bias
+ClimaAnalysis.bias(sim_var, obs_var; mask = apply_oceanmask)
+ClimaAnalysis.global_bias(sim_var, obs_var; mask = apply_oceanmask)
+
+# Do not consider the land when computing the squared error
+ClimaAnalysis.squared_error(sim_var, obs_var; mask = apply_landmask)
+ClimaAnalysis.global_mse(sim_var, obs_var; mask = apply_landmask)
+ClimaAnalysis.global_rmse(sim_var, obs_var; mask = apply_landmask)
+```

src/Var.jl

@@ -1163,7 +1163,7 @@ function _check_sim_obs_units_consistent(sim::OutputVar, obs::OutputVar)
 end
 
 """
-    bias(sim::OutputVar, obs::OutputVar)
+    bias(sim::OutputVar, obs::OutputVar; mask = nothing)
 
 Return a `OutputVar` whose data is the bias (`sim.data - obs.data`) and compute the global
 bias of `data` in `sim` and `obs` over longitude and latitude. The result is stored in
@@ -1175,10 +1175,13 @@ for longitude and latitude should be degrees. Resampling is done automatically b
 `obs` on `sim`. Attributes in `sim` and `obs` will be thrown away. The long name and short
 name of the returned `OutputVar` will be updated to reflect that a bias is computed.
 
+The parameter `mask` is a function that masks a `OutputVar`. See [`apply_landmask`](@ref)
+and [`apply_oceanmask`](@ref).
+
 See also [`global_bias`](@ref), [`squared_error`](@ref), [`global_mse`](@ref),
 [`global_rmse`](@ref).
 """
-function bias(sim::OutputVar, obs::OutputVar)
+function bias(sim::OutputVar, obs::OutputVar; mask = nothing)
     _check_sim_obs_units_consistent(sim, obs)
 
     # Resample obs on sim to ensure the size of data in sim and obs are the same and the
@@ -1203,24 +1206,24 @@ function bias(sim::OutputVar, obs::OutputVar)
     end
 
     # Compute global bias and store it as an attribute
+    !isnothing(mask) && (bias = mask(bias))
     integrated_bias = integrate_lonlat(bias).data
-    normalization =
-        integrate_lonlat(
-            OutputVar(
-                bias.attributes,
-                bias.dims,
-                bias.dim_attributes,
-                ones(size(bias.data)),
-            ),
-        ).data
+    ones_var = OutputVar(
+        bias.attributes,
+        bias.dims,
+        bias.dim_attributes,
+        ones(size(bias.data)),
+    )
+    !isnothing(mask) && (ones_var = mask(ones_var))
+    normalization = integrate_lonlat(ones_var).data
     # Do ./ instead of / because we are dividing between zero dimensional arrays
     global_bias = integrated_bias ./ normalization
     ret_attributes["global_bias"] = global_bias
     return OutputVar(ret_attributes, bias.dims, bias.dim_attributes, bias.data)
 end
 
 """
-    global_bias(sim::OutputVar, obs::OutputVar)
+    global_bias(sim::OutputVar, obs::OutputVar; mask = nothing)
 
 Return the global bias of `data` in `sim` and `obs` over longitude and latitude.
 
@@ -1229,16 +1232,19 @@ longitude and latitude. Units must be supplied for data and dimensions in `sim`
 The units for longitude and latitude should be degrees. Resampling is done automatically by
 resampling `obs` on `sim`.
 
+The parameter `mask` is a function that masks a `OutputVar`. See [`apply_landmask`](@ref)
+and [`apply_oceanmask`](@ref).
+
 See also [`bias`](@ref), [`squared_error`](@ref), [`global_mse`](@ref),
 [`global_rmse`](@ref).
 """
-function global_bias(sim::OutputVar, obs::OutputVar)
-    bias_var = bias(sim, obs)
+function global_bias(sim::OutputVar, obs::OutputVar; mask = nothing)
+    bias_var = bias(sim, obs, mask = mask)
     return bias_var.attributes["global_bias"]
 end
 
 """
-    squared_error(sim::OutputVar, obs::OutputVar)
+    squared_error(sim::OutputVar, obs::OutputVar; mask = nothing)
 
 Return a `OutputVar` whose data is the squared error (`(sim.data - obs.data)^2`) and compute
 the global mean squared error (MSE) and global root mean squared error (RMSE) of `data` in
@@ -1251,9 +1257,12 @@ for longitude and latitude should be degrees. Resampling is done automatically b
 `obs` on `sim`. Attributes in `sim` and `obs` will be thrown away. The long name and short
 name of the returned `OutputVar` will be updated to reflect that a squared error is computed.
 
+The parameter `mask` is a function that masks a `OutputVar`. See [`apply_landmask`](@ref)
+and [`apply_oceanmask`](@ref).
+
 See also [`global_mse`](@ref), [`global_rmse`](@ref), [`bias`](@ref), [`global_bias`](@ref).
 """
-function squared_error(sim::OutputVar, obs::OutputVar)
+function squared_error(sim::OutputVar, obs::OutputVar; mask = nothing)
     _check_sim_obs_units_consistent(sim, obs)
 
     # Resample obs on sim to ensure the size of data in sim and obs are the same and the
@@ -1282,16 +1291,16 @@ function squared_error(sim::OutputVar, obs::OutputVar)
     end
 
     # Compute global mse and global rmse and store it as an attribute
+    !isnothing(mask) && (squared_error = mask(squared_error))
     integrated_squared_error = integrate_lonlat(squared_error).data
-    normalization =
-        integrate_lonlat(
-            OutputVar(
-                squared_error.attributes,
-                squared_error.dims,
-                squared_error.dim_attributes,
-                ones(size(squared_error.data)),
-            ),
-        ).data
+    ones_var = OutputVar(
+        squared_error.attributes,
+        squared_error.dims,
+        squared_error.dim_attributes,
+        ones(size(squared_error.data)),
+    )
+    !isnothing(mask) && (ones_var = mask(ones_var))
+    normalization = integrate_lonlat(ones_var).data
     # Do ./ instead of / because we are dividing between zero dimensional arrays
     mse = integrated_squared_error ./ normalization
     ret_attributes["global_mse"] = mse
@@ -1305,7 +1314,7 @@ function squared_error(sim::OutputVar, obs::OutputVar)
 end
 
 """
-    global_mse(sim::OutputVar, obs::OutputVar)
+    global_mse(sim::OutputVar, obs::OutputVar; mask = nothing)
 
 Return the global mean squared error (MSE) of `data` in `sim` and `obs` over longitude and
 latitude.
@@ -1315,15 +1324,18 @@ and latitude. Units must be supplied for data and dimensions in `sim` and `obs`.
 for longitude and latitude should be degrees. Resampling is done automatically by resampling
 `obs` on `sim`.
 
+The parameter `mask` is a function that masks a `OutputVar`. See [`apply_landmask`](@ref)
+and [`apply_oceanmask`](@ref).
+
 See also [`squared_error`](@ref), [`global_rmse`](@ref), [`bias`](@ref), [`global_bias`](@ref).
 """
-function global_mse(sim::OutputVar, obs::OutputVar)
-    squared_error_var = squared_error(sim, obs)
+function global_mse(sim::OutputVar, obs::OutputVar; mask = nothing)
+    squared_error_var = squared_error(sim, obs, mask = mask)
     return squared_error_var.attributes["global_mse"]
 end
 
 """
-    global_rmse(sim::OutputVar, obs::OutputVar)
+    global_rmse(sim::OutputVar, obs::OutputVar; mask = nothing)
 
 Return the global root mean squared error (RMSE) of `data` in `sim` and `obs` over longitude
 and latitude.
@@ -1333,10 +1345,13 @@ and latitude. Units must be supplied for data and dimensions in `sim` and `obs`.
 for longitude and latitude should be degrees. Resampling is done automatically by resampling
 `obs` on `sim`.
 
+The parameter `mask` is a function that masks a `OutputVar`. See [`apply_landmask`](@ref)
+and [`apply_oceanmask`](@ref).
+
 See also [`squared_error`](@ref), [`global_mse`](@ref), [`bias`](@ref), [`global_bias`](@ref).
 """
-function global_rmse(sim::OutputVar, obs::OutputVar)
-    squared_error_var = squared_error(sim, obs)
+function global_rmse(sim::OutputVar, obs::OutputVar; mask = nothing)
+    squared_error_var = squared_error(sim, obs, mask = mask)
     return squared_error_var.attributes["global_rmse"]
 end
 

test/test_Var.jl

@@ -1566,3 +1566,114 @@ end
     @test_throws ErrorException ClimaAnalysis.apply_landmask(var)
     @test_throws ErrorException ClimaAnalysis.apply_oceanmask(var)
 end
+
+@testset "Bias and RMSE with masks" begin
+    # Test bias and global_bias
+    land_var = ClimaAnalysis.OutputVar(ClimaAnalysis.Var.LAND_MASK)
+    ocean_var = ClimaAnalysis.OutputVar(ClimaAnalysis.Var.OCEAN_MASK)
+    data_zero = zeros(land_var.data |> size)
+    zero_var = ClimaAnalysis.OutputVar(
+        land_var.attributes,
+        land_var.dims,
+        land_var.dim_attributes,
+        data_zero,
+    )
+
+    # Trim data because periodic boundary condition on the edges
+    @test ClimaAnalysis.bias(
+        land_var,
+        zero_var,
+        mask = ClimaAnalysis.apply_oceanmask,
+    ).data[
+        begin:(end - 1),
+        :,
+    ] == data_zero[begin:(end - 1), :]
+    @test ClimaAnalysis.bias(
+        ocean_var,
+        zero_var,
+        mask = ClimaAnalysis.apply_landmask,
+    ).data[
+        begin:(end - 1),
+        :,
+    ] == data_zero[begin:(end - 1), :]
+
+    # Not exactly zero because of the periodic boundary condition on the edges
+    # which results in some ones in the data
+    @test isapprox(
+        ClimaAnalysis.global_bias(
+            land_var,
+            zero_var,
+            mask = ClimaAnalysis.apply_oceanmask,
+        ),
+        0.0,
+        atol = 1e-5,
+    )
+    @test isapprox(
+        ClimaAnalysis.global_bias(
+            ocean_var,
+            zero_var,
+            mask = ClimaAnalysis.apply_landmask,
+        ),
+        0.0,
+        atol = 1e-5,
+    )
+
+    # Test squared error, global_mse, and global_rmse
+    # Trim data because periodic boundary condition on the edges
+    @test ClimaAnalysis.squared_error(
+        land_var,
+        zero_var,
+        mask = ClimaAnalysis.apply_oceanmask,
+    ).data[
+        begin:(end - 1),
+        :,
+    ] == data_zero[begin:(end - 1), :]
+    @test ClimaAnalysis.squared_error(
+        ocean_var,
+        zero_var,
+        mask = ClimaAnalysis.apply_landmask,
+    ).data[
+        begin:(end - 1),
+        :,
+    ] == data_zero[begin:(end - 1), :]
+
+    # Not exactly zero because of the periodic boundary condition on the edges
+    # which results in some ones in the data
+    @test isapprox(
+        ClimaAnalysis.global_mse(
+            land_var,
+            zero_var,
+            mask = ClimaAnalysis.apply_oceanmask,
+        ),
+        0.0,
+        atol = 1e-5,
+    )
+    @test isapprox(
+        ClimaAnalysis.global_mse(
+            ocean_var,
+            zero_var,
+            mask = ClimaAnalysis.apply_landmask,
+        ),
+        0.0,
+        atol = 1e-5,
+    )
+
+    @test isapprox(
+        ClimaAnalysis.global_rmse(
+            land_var,
+            zero_var,
+            mask = ClimaAnalysis.apply_oceanmask,
+        ),
+        0.0,
+        atol = 10^(-2.5),
+    )
+    @test isapprox(
+        ClimaAnalysis.global_rmse(
+            ocean_var,
+            zero_var,
+            mask = ClimaAnalysis.apply_landmask,
+        ),
+        0.0,
+        atol = 10^(-2.5),
+    )
+end