Add Var.integrate_* functions

This commit adds Var.integrate_lonlat, Var.integrate_lon, and
Var.integrate_lat for integrating `OutputVar` with respect to longitude,
latitude, or both. Additionally, this commit adds a Numerics module to
decouple the numerical implementation of integration and the
construction of an OutputVar. The numerical implementation is
rectangular integration using the midpoint rule when the points are
equispaced which assume each point correspond to a midpoint of a cell
and using the left endpoint rule when the points are not equispaced.

NEWS.md

@@ -120,6 +120,46 @@ contour2D_on_globe!(fig,
 CairoMakie.save("myfigure.pdf", fig)
 ```
 
+## Integrating `OutputVar` with respect to longitude or latitude
+
+You can use the `integrate_lon(var)`, `integrate_lat(var)`, or `integrate_lonlat(var)`
+functions for integrating along longitude, latitude, or both respectively. The bounds of
+integration are determined by the range of the dimensions longitude and latitude in `var`.
+The convention is that latitude ranges from -90° to 90° and longitude ranges from -180° to
+180° and the unit of both longitude and latitude should be degree. See the example of
+integrating over a sphere where the data is all ones to find the surface area of a sphere.
+
+```@julia integrate_lonlat
+julia> lon = collect(range(-179.5, 179.5, 360));
+
+julia> lat = collect(range(-89.5, 89.5, 180));
+
+julia> data = ones(length(lon), length(lat));
+
+julia> dims = OrderedDict(["lon" => lon, "lat" => lat]);
+
+julia> dim_attribs = OrderedDict([
+           "lon" => Dict("units" => "degrees_east"),
+           "lat" => Dict("units" => "degrees_north"),
+       ]);
+
+julia> attribs = Dict("long_name" => "f");
+
+julia> var = ClimaAnalysis.OutputVar(attribs, dims, dim_attribs, data);
+
+julia> integrated_var = integrate_lonlat(var);
+
+julia> integrated_var.dims # no dimensions since longitude and latitude are integrated out
+OrderedDict{String, Vector{Float64}}()
+
+julia> integrated_var.data # approximately 4π (the surface area of a sphere)
+0-dimensional Array{Float64, 0}:
+12.566530113084296
+
+julia> long_name(integrated_var) # updated long name to reflect the data being integrated
+"f integrated over lon (-179.5 to 179.5degrees_east) and integrated over lat (-89.5 to 89.5degrees_north)"
+```
+
 ## Bug fixes
 
 - Increased the default value for `warp_string` to 72.

docs/src/api.md

@@ -55,6 +55,9 @@ Var.dim_units
 Var.range_dim
 Var.resampled_as
 Var.convert_units
+Var.integrate_lonlat
+Var.integrate_lat
+Var.integrate_lon
 ```
 
 ## Utilities

docs/src/howdoi.md

@@ -107,3 +107,43 @@ OrderedDict{String, Vector{Float64}} with 2 entries:
   "lon"      => [0.0, 1.0]
   "latitude" => [0.0, 1.0, 2.0]
 ```
+
+## How do I integrate the data in a `OutputVar` with respect to longitude or latitude?
+
+You can use the `integrate_lon(var)`, `integrate_lat(var)`, or `integrate_lonlat(var)`
+functions for integrating along longitude, latitude, or both respectively. The bounds of
+integration are determined by the range of the dimensions longitude and latitude in `var`.
+The convention is that latitude ranges from -90° to 90° and longitude ranges from -180° to
+180° and the unit of both longitude and latitude should be degree. See the example of
+integrating over a sphere where the data is all ones to find the surface area of a sphere.
+
+```@julia integrate_lonlat
+julia> lon = collect(range(-179.5, 179.5, 360));
+
+julia> lat = collect(range(-89.5, 89.5, 180));
+
+julia> data = ones(length(lon), length(lat));
+
+julia> dims = OrderedDict(["lon" => lon, "lat" => lat]);
+
+julia> dim_attribs = OrderedDict([
+           "lon" => Dict("units" => "degrees_east"),
+           "lat" => Dict("units" => "degrees_north"),
+       ]);
+
+julia> attribs = Dict("long_name" => "f");
+
+julia> var = ClimaAnalysis.OutputVar(attribs, dims, dim_attribs, data);
+
+julia> integrated_var = integrate_lonlat(var);
+
+julia> integrated_var.dims # no dimensions since longitude and latitude are integrated out
+OrderedDict{String, Vector{Float64}}()
+
+julia> integrated_var.data # approximately 4π (the surface area of a sphere)
+0-dimensional Array{Float64, 0}:
+12.566530113084296
+
+julia> long_name(integrated_var) # updated long name to reflect the data being integrated
+"f integrated over lon (-179.5 to 179.5degrees_east) and integrated over lat (-89.5 to 89.5degrees_north)"
+```

src/ClimaAnalysis.jl

@@ -4,6 +4,9 @@ import Reexport: @reexport
 include("Utils.jl")
 import .Utils
 
+include("Numerics.jl")
+@reexport using .Numerics
+
 include("Var.jl")
 @reexport using .Var
 include("Sim.jl")

src/Numerics.jl

@@ -0,0 +1,132 @@
+module Numerics
+
+"""
+    isequispaced(arr::Vector)
+
+Return whether the array is equispaced or not.
+"""
+function _isequispaced(arr::Vector)
+    return all(diff(arr) .≈ arr[begin + 1] - arr[begin])
+end
+
+"""
+    _integrate_lon(data::AbstractArray, lon::AbstractVector; dims)
+
+Integrate out longitude from `data`. `data` has to be discretized on `lon`.
+
+`dims` indicates which axis of `data` is longitude.
+
+If the points are equispaced, it is assumed that each point correspond to the midpoint of a
+cell which results in rectangular integration using the midpoint rule. Otherwise, the
+integration being done is rectangular integration using the left endpoints. The unit for
+longitude should be degrees. The convention is that longitude ranges from -180° and 180°.
+"""
+function _integrate_lon(data::AbstractArray, lon::AbstractVector; dims)
+    length(lon) == 1 &&
+        error("Cannot integrate when longitude is a single point")
+    _isequispaced(lon) ?
+    int_weights = _integration_weights_lon_equispaced(lon) :
+    int_weights = _integration_weights_lon_left(lon)
+    return _integrate_over_angle(data, lon, dims, int_weights)
+end
+
+"""
+    _integrate_lat(data::AbstractArray, lat::AbstractVector; dims)
+
+Integrate out latitude from `data`. `data` has to be discretized on `lat`.
+
+`dims` indicates which axis of `data` is latitude.
+
+If the points are equispaced, it is assumed that each point correspond to the midpoint of a
+cell which results in rectangular integration using the midpoint rule. Otherwise, the
+integration being done is rectangular integration using the left endpoints. The unit for
+latitude should be degrees. The convention is that latitude ranges from -90° to 90°.
+"""
+function _integrate_lat(data::AbstractArray, lat::AbstractVector; dims)
+    length(lat) == 1 &&
+        error("Cannot integrate when latitude is a single point")
+    _isequispaced(lat) ?
+    int_weights = _integration_weights_lat_equispaced(lat) :
+    int_weights = _integration_weights_lat_left(lat)
+    return _integrate_over_angle(data, lat, dims, int_weights)
+end
+
+"""
+    _integrate_over_angle(
+    data::AbstractArray,
+    angle_arr::AbstractVector,
+    angle_idx,
+    int_weights::AbstractVector,
+)
+
+Integrate out angle (latitude or longitude) from `data` using the weights `int_weights`.
+`data` has to be discretized on `angle_arr`.
+
+`angle_idx` indicates which axis of `data` is angle.
+"""
+function _integrate_over_angle(
+    data::AbstractArray,
+    angle_arr::AbstractVector,
+    angle_idx,
+    int_weights::AbstractVector,
+)
+    # Reshape to add extra dimensions for int_weights for broadcasting if needed
+    size_to_reshape =
+        (i == angle_idx ? length(angle_arr) : 1 for i in 1:ndims(data))
+    int_weights = reshape(int_weights, size_to_reshape...)
+    int_on_angle = sum(data .* int_weights, dims = angle_idx)
+    return int_on_angle
+end
+
+"""
+    _integration_weights_lon_left(lon::AbstractVector)
+
+Return integration weights for rectangular integration using left endpoints for integrating
+along longitude.
+"""
+function _integration_weights_lon_left(lon::AbstractVector)
+    d_lon = deg2rad.(diff(lon))
+    # We are doing integration using the left endpoints, so we weight the rightmost endpoint
+    # zero so that it make no contribution to the integral
+    push!(d_lon, zero(eltype(d_lon)))
+    return d_lon
+end
+
+"""
+    _integration_weights_lat_left(lat::AbstractVector)
+
+Return integration weights for rectangular integration using left endpoints for integrating
+along latitude.
+"""
+function _integration_weights_lat_left(lat::AbstractVector)
+    d_lat = deg2rad.(diff(lat))
+    # We are doing integration using the left endpoints, so we weight the rightmost endpoint
+    # zero so that it make no contribution to the integral
+    push!(d_lat, zero(eltype(d_lat)))
+    cos_lat = cosd.(lat)
+    return d_lat .* cos_lat
+end
+
+"""
+    _integration_weights_lon_equispaced(lon::AbstractVector)
+
+Return integration weights for rectangular integration when the points are equispaced for
+integrating along longitude.
+"""
+function _integration_weights_lon_equispaced(lon::AbstractVector)
+    return deg2rad.([lon[begin + 1] - lon[begin] for _ in lon])
+end
+
+"""
+    _integration_weights_lat_equispaced(lat::AbstractVector)
+
+Return integration weights for rectangular integration when the points are equispaced for
+integrating along latitude.
+"""
+function _integration_weights_lat_equispaced(lat::AbstractVector)
+    d_lat = deg2rad.([lat[begin + 1] - lat[begin] for _ in lat])
+    cos_lat = cosd.(lat)
+    return d_lat .* cos_lat
+end
+
+end

src/Var.jl

@@ -7,6 +7,7 @@ import Interpolations as Intp
 import Statistics: mean
 import NaNStatistics: nanmean
 
+import ..Numerics
 import ..Utils: nearest_index, seconds_to_prettystr, squeeze
 
 export OutputVar,
@@ -30,7 +31,10 @@ export OutputVar,
     range_dim,
     resampled_as,
     has_units,
-    convert_units
+    convert_units,
+    integrate_lonlat,
+    integrate_lon,
+    integrate_lat
 
 """
     Representing an output variable
@@ -785,6 +789,105 @@ function resampled_as(src_var::OutputVar, dest_var::OutputVar)
     )
 end
 
+"""
+    integrate_lonlat(var::OutputVar)
+
+Integrate `data` in `var` on longitude and latitude with a first-order scheme. `data` has to
+be discretized on longitude and latitude.
+
+If the points are equispaced, it is assumed that each point correspond to the midpoint of a
+cell which results in rectangular integration using the midpoint rule. Otherwise, the
+integration being done is rectangular integration using the left endpoints for integrating
+longitude and latitude. The units for longitude and latitude should be degrees. The
+convention is that longitude vary from -180° and 180° and latitude vary from -90° to 90°.
+"""
+function integrate_lonlat(var::OutputVar)
+    var_integrate_lon = var |> integrate_lon
+    # Update long name so that we get "...integrated over lat... and integrated over lon..."
+    # instead of "...integrated over lat... integrated over lon..."
+    if haskey(var_integrate_lon.attributes, "long_name")
+        var_integrate_lon.attributes["long_name"] *= " and"
+    end
+    return var_integrate_lon |> integrate_lat
+end
+
+"""
+    integrate_lon(var::OutputVar)
+
+Integrate `data` in `var` on longitude with a first-order scheme. `data` has to be
+discretized on longitude.
+
+If the points are equispaced, it is assumed that each point correspond to the midpoint of a
+cell which results in rectangular integration using the midpoint rule. Otherwise, the
+integration being done is rectangular integration using the left endpoints. The unit for
+longitude should be degrees. The convention is that longitude ranges from -180° and 180°.
+"""
+function integrate_lon(var::OutputVar)
+    has_longitude(var) || error("var does not has longitude as a dimension")
+    lon_name = longitude_name(var)
+    return _integrate_over_angle(var, Numerics._integrate_lon, lon_name)
+end
+
+"""
+    integrate_lat(var::OutputVar)
+
+Integrate `data` in `var` on latitude with a first-order scheme. `data` has to be
+discretized on latitude.
+
+If the points are equispaced, it is assumed that each point correspond to the midpoint of a
+cell which results in rectangular integration using the midpoint rule. Otherwise, the
+integration being done is rectangular integration using the left endpoints. The unit for
+latitude should be degrees. The convention is that latitude ranges from -90° to 90°.
+"""
+function integrate_lat(var::OutputVar)
+    has_latitude(var) || error("var does not has latitude as a dimension")
+    lat_name = latitude_name(var)
+    return _integrate_over_angle(var, Numerics._integrate_lat, lat_name)
+end
+
+"""
+    _integrate_over_angle(var::OutputVar, integrate_on, angle_dim_name)
+
+Integrate `data` in `var` on `angle_dim_name` with a first-order scheme. `data` has to be
+discretized on `angle_dim_name`.
+
+`angle_dim_name` is the name of the angle that is being integrated over. `integrate_on` is a
+function that compute the integration for data over `angle_dim_name`.
+"""
+function _integrate_over_angle(var::OutputVar, integrate_on, angle_dim_name)
+    # Enforce constraint that unit is degree because we compute integration weights assuming
+    # degrees (see functions _integration_weights_lon_left and
+    # _integration_weights_lon_equispaced for examples in Numerics.jl)
+    deg_unit_names = [
+        "degrees",
+        "degree",
+        "deg",
+        "degs",
+        "°",
+        "degrees_north",
+        "degrees_east",
+    ]
+    angle_dim_unit = dim_units(var, angle_dim_name)
+    angle_dim_unit in deg_unit_names ||
+        lowercase(angle_dim_unit) |> (x -> occursin("deg", x)) ||
+        error("The unit for $angle_dim_name is missing or is not degree")
+
+    integrated_var = _reduce_over(
+        integrate_on,
+        angle_dim_name,
+        var,
+        var.dims[angle_dim_name],
+    )
+
+    _update_long_name_generic!(
+        integrated_var,
+        var,
+        angle_dim_name,
+        "integrated",
+    )
+    return integrated_var
+end
+
 """
     overload_binary_op(op)
 

test/runtests.jl

@@ -2,15 +2,16 @@ using SafeTestsets
 using Test
 
 #! format: off
-# @safetestset "Aqua" begin @time include("aqua.jl") end
-# @safetestset "Docstrings" begin @time include("doctest.jl") end
-# @safetestset "Format" begin @time include("format.jl") end
+@safetestset "Aqua" begin @time include("aqua.jl") end
+@safetestset "Docstrings" begin @time include("doctest.jl") end
+@safetestset "Format" begin @time include("format.jl") end
 
-# @safetestset "Utils" begin @time include("test_Utils.jl") end
-# @safetestset "SimDir" begin @time include("test_Sim.jl") end
-# @safetestset "Atmos" begin @time include("test_Atmos.jl") end
-# @safetestset "OutputVar" begin @time include("test_Var.jl") end
-# @safetestset "MakieExt" begin @time include("test_MakieExt.jl") end
+@safetestset "Utils" begin @time include("test_Utils.jl") end
+@safetestset "Numerics" begin @time include("test_Numerics.jl") end
+@safetestset "SimDir" begin @time include("test_Sim.jl") end
+@safetestset "Atmos" begin @time include("test_Atmos.jl") end
+@safetestset "OutputVar" begin @time include("test_Var.jl") end
+@safetestset "MakieExt" begin @time include("test_MakieExt.jl") end
 @safetestset "GeoMakieExt" begin @time include("test_GeoMakieExt.jl") end
 #! format: on