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,50 @@ 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 unit of both longitude and latitude should be degree.
+
+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. 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/var.md

@@ -54,3 +54,49 @@ new_var = ClimaAnalysis.convert_units(var, "kg m/s", conversion_function = (x) -
 ```
 
 !!! note If you find some unparseable units, please open an issue. We can fix them!
+
+## Integration
+
+`OutputVar`s can be integrated with respect to longitude, latitude, or both using
+`integrate_lon(var)`, `integrate_lat(var)`, or `integrate_lonlat(var)` respectively. The
+bounds of integration are determined by the range of the dimensions longitude and latitude
+in `var`. The unit of both longitude and latitude should be degree.
+
+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.
+
+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,8 @@ import Reexport: @reexport
 include("Utils.jl")
 import .Utils
 
+include("Numerics.jl")
+
 include("Var.jl")
 @reexport using .Var
 include("Sim.jl")

src/Numerics.jl

@@ -0,0 +1,128 @@
+module Numerics
+
+import ..Utils: _isequispaced
+
+"""
+    _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.
+"""
+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.
+"""
+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,
+)
+    # Reshape to add extra dimensions for int_weights for broadcasting if needed
+    size_to_reshape =
+        (i == angle_idx ? length(int_weights) : 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)
+
+Return integration weights for rectangular integration using left endpoints for integrating
+along longitude.
+"""
+function _integration_weights_lon_left(lon)
+    # This is where we use the assumption that units are degrees
+    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)
+
+Return integration weights for rectangular integration using left endpoints for integrating
+along latitude.
+"""
+function _integration_weights_lat_left(lat)
+    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)
+
+Return integration weights for rectangular integration when the points are equispaced for
+integrating along longitude.
+"""
+function _integration_weights_lon_equispaced(lon)
+    # This is where we use the assumption that units are degrees
+    # Use fill to make a zero dimensional array so reshaping is possible
+    return fill(deg2rad(lon[begin + 1] - lon[begin]))
+end
+
+"""
+    _integration_weights_lat_equispaced(lat)
+
+Return integration weights for rectangular integration when the points are equispaced for
+integrating along latitude.
+"""
+function _integration_weights_lat_equispaced(lat)
+    d_lat = deg2rad.(lat[begin + 1] - lat[begin])
+    cos_lat = cosd.(lat)
+    return d_lat .* cos_lat
+end
+
+end

src/Utils.jl

@@ -296,4 +296,24 @@ function split_by_season(dates::AbstractArray{<:Dates.DateTime})
     return (MAM, JJA, SON, DJF)
 end
 
+"""
+    _isequispaced(arr::Vector)
+
+Return whether the array is equispaced or not.
+
+Examples
+=========
+
+```jldoctest
+julia> Utils._isequispaced([1.0, 2.0, 3.0])
+true
+
+julia> Utils._isequispaced([0.0, 2.0, 3.0])
+false
+```
+"""
+function _isequispaced(arr::Vector)
+    return all(diff(arr) .≈ arr[begin + 1] - arr[begin])
+end
+
 end