Add Var.integrate* functions

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.

Project.toml

@@ -1,6 +1,6 @@
 name = "ClimaAnalysis"
 uuid = "29b5916a-a76c-4e73-9657-3c8fd22e65e6"
-authors = ["Gabriele Bozzola <gbozzola@caltech.edu>"]
+authors = ["Gabriele Bozzola <gbozzola@caltech.edu>", "Kevin Phan <kphan2@caltech.edu>"]
 version = "0.5.7"
 
 [deps]

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