From 464d27314fbfe81bec6d1b460266be769f55b8ec Mon Sep 17 00:00:00 2001
From: Kevin Phan <98072684+ph-kev@users.noreply.github.com>
Date: Wed, 4 Sep 2024 21:19:29 -0700
Subject: [PATCH] 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               |  40 +++++++++++++
 docs/src/api.md       |   3 +
 docs/src/howdoi.md    |  40 +++++++++++++
 src/ClimaAnalysis.jl  |   3 +
 src/Numerics.jl       | 132 +++++++++++++++++++++++++++++++++++++++++
 src/Var.jl            | 105 ++++++++++++++++++++++++++++++++-
 test/runtests.jl      |  17 +++---
 test/test_Numerics.jl | 113 +++++++++++++++++++++++++++++++++++
 test/test_Var.jl      | 134 ++++++++++++++++++++++++++++++++++++++++++
 9 files changed, 578 insertions(+), 9 deletions(-)
 create mode 100644 src/Numerics.jl
 create mode 100644 test/test_Numerics.jl

diff --git a/NEWS.md b/NEWS.md
index 6c8ac53..055aa41 100644
--- a/NEWS.md
+++ b/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.
diff --git a/docs/src/api.md b/docs/src/api.md
index d227772..9c7a1c6 100644
--- a/docs/src/api.md
+++ b/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
diff --git a/docs/src/howdoi.md b/docs/src/howdoi.md
index 3a42088..42cae81 100644
--- a/docs/src/howdoi.md
+++ b/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)"
+```
\ No newline at end of file
diff --git a/src/ClimaAnalysis.jl b/src/ClimaAnalysis.jl
index aff9deb..229179c 100644
--- a/src/ClimaAnalysis.jl
+++ b/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")
diff --git a/src/Numerics.jl b/src/Numerics.jl
new file mode 100644
index 0000000..2caabed
--- /dev/null
+++ b/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
diff --git a/src/Var.jl b/src/Var.jl
index d3395bb..f2f2965 100644
--- a/src/Var.jl
+++ b/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)
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 94f54bf..0e65c89 100644
--- a/test/runtests.jl
+++ b/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
 
diff --git a/test/test_Numerics.jl b/test/test_Numerics.jl
new file mode 100644
index 0000000..ad11d87
--- /dev/null
+++ b/test/test_Numerics.jl
@@ -0,0 +1,113 @@
+using Test
+import ClimaAnalysis
+
+@testset "integration weights for lon and lat" begin
+    # Integration weights for lon (not equispaced)
+    lon = [-180.0, -45.0, 100.0, 180.0]
+    lon_weights = [135.0, 145.0, 80.0, 0.0] .* (π / 180.0)
+    @test all(
+        isapprox.(
+            lon_weights,
+            ClimaAnalysis.Numerics._integration_weights_lon_left(lon),
+        ),
+    )
+
+    # Integration weights for lat (not equispaced)
+    lat = [-90.0, 20.0, 45.0, 90.0]
+    lat_weights = [110.0, 25.0, 45.0, 0.0] .* (π / 180.0) .* cosd.(lat)
+    @test all(
+        isapprox.(
+            lat_weights,
+            ClimaAnalysis.Numerics._integration_weights_lat_left(lat),
+        ),
+    )
+
+    # Integration weights for lon (not equispaced)
+    lon = collect(range(-180.0, 180.0, 5))
+    lon_weights = [90.0 for _ in lon] .* (π / 180.0)
+    @test all(
+        isapprox.(
+            lon_weights,
+            ClimaAnalysis.Numerics._integration_weights_lon_equispaced(lon),
+        ),
+    )
+
+    lat = collect(range(-90.0, 90.0, 5))
+    lat_weights = [45.0 for _ in lat] .* (π / 180.0) .* cosd.(lat)
+    @test all(
+        isapprox.(
+            lat_weights,
+            ClimaAnalysis.Numerics._integration_weights_lat_equispaced(lat),
+        ),
+    )
+end
+
+@testset "Integrating on lon and lat" begin
+    # Integrating only lon (non equispaced)
+    lon = collect(range(-180.0, 179.0, 100))
+    # Force array to be non equispaced for testing _integration_weights_lon
+    push!(lon, 180.0)
+    lon_data = ones(length(lon))
+    @test isapprox(
+        ClimaAnalysis.Numerics._integrate_lon(lon_data, lon, dims = 1)[1],
+        2.0π,
+        atol = 0.01,
+    )
+
+    # Integrating only lat (non equispaced)
+    lat = collect(range(-90.0, 89.0, 100))
+    # Force array to be non equispaced for testing _integration_weights_lat
+    push!(lat, 90.0)
+    lat_data = ones(length(lat))
+    @test isapprox(
+        ClimaAnalysis.Numerics._integrate_lat(lat_data, lat, dims = 1)[1],
+        2.0,
+        atol = 0.01,
+    )
+
+    # Integrating both lon and lat
+    data = ones(length(lat), length(lon))
+    integrated_lat = ClimaAnalysis.Numerics._integrate_lat(data, lat, dims = 1)
+    integrated_latlon =
+        ClimaAnalysis.Numerics._integrate_lon(integrated_lat, lon, dims = 1)
+
+    integrated_lon = ClimaAnalysis.Numerics._integrate_lon(data, lon, dims = 2)
+    integrated_lonlat =
+        ClimaAnalysis.Numerics._integrate_lat(integrated_lon, lat, dims = 1)
+
+    # Order of integration should not matter
+    @test isapprox(integrated_latlon[1], integrated_lonlat[1])
+    @test isapprox(integrated_latlon[1], 4π, atol = 0.01)
+
+    # Error checking for length of lon and lat and values in lon and lat
+    @test_throws "Cannot integrate when latitude is a single point" ClimaAnalysis.Numerics._integrate_lat(
+        lat_data,
+        [0.0],
+        dims = 1,
+    )
+
+    @test_throws "Cannot integrate when latitude is a single point" ClimaAnalysis.Numerics._integrate_lat(
+        lon_data,
+        [0.0],
+        dims = 1,
+    )
+
+    # Integrating only lon (non equispaced)
+    lon = collect(range(-179.5, 179.5, 360))
+    lon_data = ones(length(lon))
+    @test isapprox(
+        ClimaAnalysis.Numerics._integrate_lon(lon_data, lon, dims = 1)[1],
+        2.0π,
+        atol = 0.01,
+    )
+
+    # Integrating only lat (non equispaced)
+    lat = collect(range(-89.5, 89.5, 180))
+    lat_data = ones(length(lat))
+    @test isapprox(
+        ClimaAnalysis.Numerics._integrate_lat(lat_data, lat, dims = 1)[1],
+        2.0,
+        atol = 0.01,
+    )
+
+end
diff --git a/test/test_Var.jl b/test/test_Var.jl
index 7243b96..d256baf 100644
--- a/test/test_Var.jl
+++ b/test/test_Var.jl
@@ -687,3 +687,137 @@ end
     )
 
 end
+
+@testset "Integrating on lat and lon" begin
+    # Tests for integrate_lon
+    lon = collect(range(-179.5, 179.5, 360))
+    lon_data = ones(length(lon))
+    lon_dims = OrderedDict(["lon" => lon])
+    lon_attribs = Dict("long_name" => "hi")
+    lon_dim_attribs = OrderedDict(["lon" => Dict("units" => "deg")])
+    var = ClimaAnalysis.OutputVar(
+        lon_attribs,
+        lon_dims,
+        lon_dim_attribs,
+        lon_data,
+    )
+    var_integrated_lon = ClimaAnalysis.Var.integrate_lon(var)
+
+    @test isapprox(var_integrated_lon.data[1], 2.0 * π, atol = 0.01)
+    @test var_integrated_lon.dims == OrderedDict()
+    @test var_integrated_lon.dim_attributes == OrderedDict()
+    @test "hi integrated over lon (-179.5 to 179.5deg)" ==
+          var_integrated_lon.attributes["long_name"]
+    @test_throws "var does not has latitude as a dimension" ClimaAnalysis.Var.integrate_lat(
+        var,
+    )
+
+    # Tests for integrate_lat
+    lat = collect(range(-89.5, 89.5, 180))
+    lat_data = ones(length(lat))
+    lat_dims = OrderedDict(["lat" => lat])
+    lat_attribs = Dict("long_name" => "hi")
+    lat_dim_attribs = OrderedDict(["lat" => Dict("units" => "deg")])
+    var = ClimaAnalysis.OutputVar(
+        lat_attribs,
+        lat_dims,
+        lat_dim_attribs,
+        lat_data,
+    )
+    var_integrated_lat = ClimaAnalysis.Var.integrate_lat(var)
+
+    @test isapprox(var_integrated_lat.data[1], 2.0, atol = 0.01)
+    @test var_integrated_lat.dims == OrderedDict()
+    @test var_integrated_lat.dim_attributes == OrderedDict()
+    @test "hi integrated over lat (-89.5 to 89.5deg)" ==
+          var_integrated_lat.attributes["long_name"]
+    @test_throws "var does not has longitude as a dimension" ClimaAnalysis.Var.integrate_lon(
+        var,
+    )
+
+    # Unit checking
+    dim_attribs_no_units = OrderedDict([
+        "lon" => Dict("units" => ""),
+        "lat" => Dict("units" => ""),
+    ])
+    var_lon_no_units = ClimaAnalysis.OutputVar(
+        lon_attribs,
+        lon_dims,
+        dim_attribs_no_units,
+        lon_data,
+    )
+    @test_throws ErrorException ClimaAnalysis.Var.integrate_lon(
+        var_lon_no_units,
+    )
+
+    var_lat_no_units = ClimaAnalysis.OutputVar(
+        lat_attribs,
+        lat_dims,
+        dim_attribs_no_units,
+        lat_data,
+    )
+    @test_throws ErrorException ClimaAnalysis.Var.integrate_lat(
+        var_lat_no_units,
+    )
+end
+
+@testset "Integrating on sphere" begin
+    # Integrate out all dimensions (lat and lon) from OutputVar
+    lon = collect(range(-179.5, 179.5, 360))
+    lat = collect(range(-89.5, 89.5, 180))
+    data = ones(length(lon), length(lat))
+    dims = OrderedDict(["lon" => lon, "lat" => lat])
+    attribs = Dict("long_name" => "hi")
+    dim_attribs = OrderedDict([
+        "lon" => Dict("units" => "deg"),
+        "lat" => Dict("units" => "deg"),
+    ])
+    var = ClimaAnalysis.OutputVar(attribs, dims, dim_attribs, data)
+    integrated_var = ClimaAnalysis.Var.integrate_lonlat(var)
+    @test isapprox(integrated_var.data[1], 4 * π, atol = 0.1)
+    @test integrated_var.dims == OrderedDict()
+    @test integrated_var.dim_attributes == OrderedDict()
+
+    # Integrating out lon and lat to get time series data
+    lon = collect(range(-179.5, 179.5, 360))
+    lat = collect(range(-89.5, 89.5, 180))
+    time = collect(range(0.0, 10.0, 10))
+    data = ones(length(lat), length(time), length(lon))
+    dims = OrderedDict(["lat" => lat, "time" => time, "lon" => lon])
+    attribs = Dict("long_name" => "hi")
+    dim_attribs = OrderedDict([
+        "lat" => Dict("units" => "deg"),
+        "time" => Dict("units" => "days"),
+        "lon" => Dict("units" => "deg"),
+    ])
+    var = ClimaAnalysis.OutputVar(attribs, dims, dim_attribs, data)
+    integrated_var = ClimaAnalysis.Var.integrate_lonlat(var)
+
+    @test all(
+        isapprox.(integrated_var.data, [4 * π for _ in 1:10], atol = 0.01),
+    )
+    @test integrated_var.dims == OrderedDict(["time" => time])
+    @test integrated_var.dim_attributes ==
+          OrderedDict(["time" => Dict("units" => "days")])
+    @test "hi integrated over lon (-179.5 to 179.5deg) and integrated over lat (-89.5 to 89.5deg)" ==
+          integrated_var.attributes["long_name"]
+
+    # Unit checking
+    dim_attribs_no_lon = OrderedDict([
+        "time" => Dict("units" => "days"),
+        "lat" => Dict("units" => "deg"),
+    ])
+    var = ClimaAnalysis.OutputVar(attribs, dims, dim_attribs_no_lon, data)
+    @test_throws "The unit for lon is missing or is not degrees" ClimaAnalysis.Var.integrate_lonlat(
+        var,
+    )
+
+    dim_attribs_no_lat = OrderedDict([
+        "time" => Dict("units" => "days"),
+        "lon" => Dict("units" => "deg"),
+    ])
+    var = ClimaAnalysis.OutputVar(attribs, dims, dim_attribs_no_lat, data)
+    @test_throws "The unit for lat is missing or is not degrees" ClimaAnalysis.Var.integrate_lonlat(
+        var,
+    )
+end