Introduction

ClimArray(A::Array, dims::Tuple; name = "", attrib = nothing)

ClimArray is a structure that contains numerical array data bundled with dimensional information, a name and an attrib field (typically a dictionary) that holds general attributes. You can think of ClimArray as a in-memory representation of a CFVariable.

At the moment, a ClimArray is using DimensionalArray from DimensionalData.jl, and all basic handling of ClimArray is offered by DimensionalData (see below).

ClimArray is created by passing in standard array data A and a tuple of dimensions dims.

Example

using ClimateBase, Dates
Time = ClimateBase.Ti # more intuitive name for time dimension
lats = -90:5:90
lons = 0:10:359
t = Date(2000, 3, 15):Month(1):Date(2020, 3, 15)
# dimensional information:
dimensions = (Lon(lons), Lat(lats), Time(t))
data = rand(36, 37, 241) # numeric data
A = ClimArray(data, dimensions)

ClimArray(file::Union{String,NCDataset}, var::String, name = var) -> A

Load the variable var from the file and convert it into a ClimArray which also contains the variable attributes as a dictionary. Dimension attributes are also given to the dimensions of A, if any exist.

Notice that file can be an NCDataset, which allows you to lazily combine different .nc data (typically split by time), e.g.

alldata = ["toa_fluxes_2020_$(i).nc" for i in 1:12]
file = NCDataset(alldata; aggdim = "time")
A = ClimArray(file, "tow_sw_all")

(but you can also directly give the string to a single file "file.nc" in ClimArray if data are contained in a single file for single files).

We do two performance improvements while loading the data:

1. If there are no missing values in the data (according to CF standards), the returned array is automatically converted to a concrete type (i.e. Union{Float32, Missing} becomes Float32).
2. Dimensions that are ranges (i.e. sampled with constant step size) are automatically transformed to a standard Julia Range type (which makes sub-selecting faster).

nckeys(file::String)

Return all keys of the .nc file in file.

ncdetails(file::String, io = stdout)

Print details about the .nc file in file on io.

globalattr(file::String) → Dict

Return the global attributes of the .nc file.

climarrays_to_nc(file::String, Xs; globalattr = Dict())

Write the given ClimArray instances (any iterable of ClimArrays or a single ClimArray) to a .nc file following CF standard conventions using NCDatasets.jl. Optionally specify global attributes for the .nc file.

The metadata of the arrays in Xs, as well as their dimensions, are properly written in the .nc file and any necessary type convertions happen automatically.

WARNING: We assume that any dimensions shared between the Xs are identical.

timemean(A::ClimArray [, w]) = timeagg(mean, A, w)

Temporal average of A, see timeagg.

timeagg(f, A::ClimArray, W = nothing)

Perform a proper temporal aggregation of the function f (e.g. mean, std) on A where:

• Only full year spans of A are included, see maxyearspan (because most processes are affected by yearly cycle, and averaging over an uneven number of cycles typically results in artifacts)
• Each month in A is weighted with its length in days (for monthly sampled data)

If you don't want these features, just do dropagg(f, A, Time, W).

W are possible statistical weights that are used in conjuction to the temporal weighting, to weight each time point differently. If they are not a vector (a weight for each time point), then they have to be a dimensional array of same dimensional layout as A (a weight for each data point).

See also monthlyagg, yearlyagg, seasonalyagg.

timeagg(f, t::Vector, x::Vector, w = nothing)

Same as above, but for arbitrary vector x accompanied by time vector t.

monthlyagg(A::ClimArray, f = mean; mday = 15) -> B

Create a new array where the temporal information has been aggregated into months using the function f. The dates of the new array always have day number of mday.

yearlyagg(A::ClimArray, f = mean) -> B

Create a new array where the temporal information has been aggregated into years using the function f. By convention, the dates of the new array always have month and day number of 1.

seasonalyagg(A::ClimArray, f = mean) -> B

Create a new array where the temporal information has been aggregated into seasons using the function f. By convention, seasons are represented as Dates spaced 3-months apart, where only the months December, March, June and September are used to specify the date, with day 1.

temporalrange(A::ClimArray, f = Dates.month) → r
temporalrange(t::AbstractVector{<:TimeType}}, f = Dates.month) → r

Return a vector of ranges so that each range of indices are values of t that belong in either the same month, year, day, or season, depending on f. f can take the values: Dates.year, Dates.month, Dates.day or season (all are functions).

Used in e.g. monthlyagg, yearlyagg or seasonalyagg.

maxyearspan(A::ClimArray) = maxyearspan(dims(A, Time))
maxyearspan(t::Vector{<:DateTime}) → i

Find the maximum index i of t so that t[1:i] covers exact(*) multiples of years.

(*) For monthly spaced data i is a multiple of 12 while for daily data we find the largest possible multiple of DAYS_IN_ORBIT.

temporal_sampling(x) → symbol

Return the temporal sampling type of x, which is either an array of Dates or a dimensional array (with Time dimension).

Possible return values are:

• :hourly, where the temporal difference between entries is exactly 1 hour.
• :daily, where the temporal difference between dates is exactly 1 day.
• :monthly, where all dates have the same day, but different month.
• :yearly, where all dates have the same month+day, but different year.
• :other, which means that x doesn't fall to any of the above categories.

time_in_days(t::AbstractArray{<:TimeType}, T = Float32)

Convert a given date time array into measurement units of days: a real-valued array which counts time in days, always increasing (cumulative).

zonalmean(A::ClimArray)

Return the zonal mean of A.

latmean(A::ClimArray)

Return the latitude-mean A (mean across dimension Lat). This function properly weights by the cosine of the latitude.

spacemean(A::ClimArray [, W]) = spaceagg(mean, A, W)

Average given A over its spatial coordinates. Optionally provide statistical weights in W.

spaceagg(f, A::ClimArray, W = nothing)

Aggregate A using function f (e.g. mean, std) over all available space (i.e. longitude and latitude) of A, weighting every part of A by its spatial area.

W can be extra weights, to weight each spatial point with. W can either be just a ClimArray with same space as A, or of exactly same shape as A.

hemispheric_means(A) → nh, sh

Return the (proper) averages of A over the northern and southern hemispheres. Notice that this function explicitly does both zonal as well as meridional averaging. Use hemispheric_functions to just split A into two hemispheres.

hemispheric_functions(A::ClimArray) → north, south

Return two arrays north, south, by splitting A to its northern and southern hemispheres, appropriately translating the latitudes of south so that both arrays have the same latitudinal dimension (and thus can be compared and do opperations between them).

lonlatfirst(A::ClimArray, args...) → B

Permute the dimensions of A to make a new array B that has first dimension longitude, second dimension latitude, with the remaining dimensions of A following (useful for most plotting functions). Optional extra dimensions can be given as args..., specifying a specific order for the remaining dimensions.

Example:

B = lonlatfirst(A)
C = lonlatfirst(A, Time)

spatialidxs(A::ClimArray) → idxs

Return an iterable that can be used to access all spatial points of A with the syntax

idxs = spatialidxs(A)
for i in idxs
    slice_at_give_space_point = A[i...]
end

Works for all types of space (... is necessary because i is a Tuple).

dropagg(f, A, dims)

Apply statistics/aggregating function f (e.g. sum or mean) on array A across dimension(s) dims and drop the corresponding dimension(s) from the result (Julia inherently keeps singleton dimensions).

If A is one dimensional, dropagg will return the single number of applying f(A).

collapse(f, A, dim)

Reduce A towards dimension dim using the collapsing function f (e.g. mean). This means that f is applied across all other dimensions of A, each of which are subsequently dropped, leaving only the collapsed result of A vs. the remaining dimension.

sinusoidal_continuation(T::ClimArray, fs = [1, 2]; Tmin = -Inf, Tmax = Inf)

Fill-in the missing values of spatiotemporal field T, by fitting sinusoidals to the non-missing values, and then using the fitted sinusoidals for the missing values.

Frequencies are given per year (frequency 2 means 1/2 of a year).

Tmin, Tmax limits are used to clamp the result into this range (no clamping in the default case).

seasonal_decomposition(A::ClimArray, fs = [1, 2]) → seasonal, residual

Decompose A into a seasonal and residual components, where the seasonal contains the periodic parts of A, with frequencies given in fs, and residual contains what's left.

fs is measured in 1/year. This function works even for non-equispaced time axis (e.g. monthly averages) and uses LPVSpectral.jl and SignalDecomposition.jl.

insolation(t, ϕ; kwargs...)

Calculate daily averaged insolation in W/m² at given time and latitude t, φ. φ is given in degrees, and t in days (real number or date).

Keywords:

Ya = DAYS_IN_ORBIT # = 365.26 # days
t_VE = 76.0 # days of vernal equinox
S_0 = 1362.0 # W/m^2
γ=23.44
ϖ=282.95
e=0.017 # eccentricity

surface_atmosphere_contributions(I, F_toa_⬆, F_s_⬆, F_s_⬇) → α_ATM, α_SFC

Calculate the atmospheric and surface contributions of the planetary albedo, so that the TOA albedo is α = α_ATM + α_SFC, using the simple 1-layer radiative transfer model by Donohoe & Battisti (2011) or G. Stephens (2015). Stephens' formulas are incorrect and I have corrected them!

total_toa_albedo(a, s, t) = a + s*t^2/(1-a*s)

Combine given atmosphere albedo a, surface albedo s and atmosphere transmittance t into a total top-of-the-atmosphere albedo α according to the model of Donohoe & Battisti (2011).

DimensionalData

DimensionalData.jl provides tools and abstractions for working with datasets that have named dimensions. It's a pluggable, generalised version of AxisArrays.jl with a cleaner syntax, and additional functionality found in NamedDims.jl. It has similar goals to pythons xarray, and is primarily written for use with spatial data in GeoData.jl.

Dimensions

Dimensions are just wrapper types. They store the dimension index and define details about the grid and other metadata, and are also used to index into the array, wrapping a value or a Selector. X, Y, Z and Ti are the exported defaults.

A generalised Dim type is available to use arbitrary symbols to name dimensions. Custom dimensions can be defined using the @dim macro.

We can use dim wrappers for indexing, so that the dimension order in the underlying array does not need to be known:

julia> using DimensionalData

julia> A = DimArray(rand(40, 50), (X, Y));

julia> A[Y(1), X(1:10)]
DimArray with dimensions:
 X: 1:10 (NoIndex)
and referenced dimensions:
 Y: 1 (NoIndex)
and data: 10-element Array{Float64,1}
[0.515774, 0.575247, 0.429075, 0.234041, 0.4484, 0.302562, 0.911098, 0.541537, 0.267234, 0.370663]

And this has no runtime cost:

julia> A = DimArray(rand(40, 50), (X, Y));

julia> @btime $A[X(1), Y(2)]
  2.092 ns (0 allocations: 0 bytes)
0.27317596504655417

julia> @btime parent($A)[1, 2]
  2.092 ns (0 allocations: 0 bytes)
0.27317596504655417

Dims can be used for indexing and views without knowing dimension order:

julia> A[X(10)]
DimArray with dimensions:
 Y (type Y): Base.OneTo(50) (NoIndex)
and referenced dimensions:
 X (type X): 10 (NoIndex)
and data: 50-element Array{Float64,1}
[0.0850249, 0.313408, 0.0762157, 0.549103, 0.297763, 0.309075, 0.854535, 0.659537, 0.392969, 0.89998  …  0.63791, 0.875881, 0.437688, 0.925918, 0.291636, 0.358024, 0.692283, 0.606932, 0.629122, 0.284592]

julia> view(A, Y(30:40), X(1:5))
DimArray with dimensions:
 X (type X): 1:5 (NoIndex)
 Y (type Y): 30:40 (NoIndex)
and data: 5×11 view(::Array{Float64,2}, 1:5, 30:40) with eltype Float64
 0.508793   0.721117  0.558849  …  0.505518   0.532322
 0.869126   0.754219  0.328315     0.0148934  0.778308
 0.0596468  0.458492  0.250458     0.980508   0.524938
 0.446838   0.659638  0.632399     0.33478    0.549402
 0.292962   0.995038  0.26026      0.526124   0.589176

And for indicating dimensions to reduce or permute in julia Base and Statistics functions that have dims arguments:

julia> using Statistics

julia> A = DimArray(rand(3, 4, 5), (X, Y, Ti));

julia> mean(A, dims=Ti)
DimArray with dimensions:
 X (type X): Base.OneTo(3) (NoIndex)
 Y (type Y): Base.OneTo(4) (NoIndex)
 Time (type Ti): 1 (NoIndex)
and data: 3×4×1 Array{Float64,3}
[:, :, 1]
 0.495295  0.650432  0.787521  0.502066
 0.576573  0.568132  0.770812  0.504983
 0.39432   0.5919    0.498638  0.337065
[and 0 more slices...]

julia> permutedims(A, [Ti, Y, X])
DimArray with dimensions:
 Time (type Ti): Base.OneTo(5) (NoIndex)
 Y (type Y): Base.OneTo(4) (NoIndex)
 X (type X): Base.OneTo(3) (NoIndex)
and data: 5×4×3 Array{Float64,3}
[:, :, 1]
 0.401374  0.469474  0.999326  0.265688
 0.439387  0.57274   0.493883  0.88678
 0.425845  0.617372  0.998552  0.650999
 0.852777  0.954702  0.928367  0.0045136
 0.357095  0.637873  0.517476  0.702351
[and 2 more slices...]

You can also use symbols to create Dim{X} dimensions:

julia> A = DimArray(rand(10, 20, 30), (:a, :b, :c));

julia> A[a=2:5, c=9]

DimArray with dimensions:
 Dim{:a}: 2:5 (NoIndex)
 Dim{:b}: Base.OneTo(20) (NoIndex)
and referenced dimensions:
 Dim{:c}: 9 (NoIndex)
and data: 4×20 Array{Float64,2}
 0.868237   0.528297   0.32389   …  0.89322   0.6776    0.604891
 0.635544   0.0526766  0.965727     0.50829   0.661853  0.410173
 0.732377   0.990363   0.728461     0.610426  0.283663  0.00224321
 0.0849853  0.554705   0.594263     0.217618  0.198165  0.661853

Selectors

Selectors find indices in the dimension based on values At, Near, or Between the index value(s). They can be used in getindex, setindex! and view to select indices matching the passed in value(s)

• At(x): get indices exactly matching the passed in value(s)
• Near(x): get the closest indices to the passed in value(s)
• Where(f::Function): filter the array axis by a function of dimension index values.
• Between(a, b): get all indices between two values (inclusive)
• Contains(x): get indices where the value x falls in the interval. Only used for Sampled Intervals, for Points us At.

We can use selectors with dim wrappers:

using Dates, DimensionalData
timespan = DateTime(2001,1):Month(1):DateTime(2001,12)
A = DimArray(rand(12,10), (Ti(timespan), X(10:10:100)))

julia> A[X(Near(35)), Ti(At(DateTime(2001,5)))]
0.658404535807791

julia> A[Near(DateTime(2001, 5, 4)), Between(20, 50)]
DimArray with dimensions:
 X: 20:10:50
and referenced dimensions:
 Time (type Ti): 2001-05-01T00:00:00
and data: 4-element Array{Float64,1}
[0.456175, 0.737336, 0.658405, 0.520152]

Without dim wrappers selectors must be in the right order:

using Unitful

julia> A = DimArray(rand(10, 20), (X((1:10:100)u"m"), Ti((1:5:100)u"s")))

julia> A[Between(10.5u"m", 50.5u"m"), Near(23u"s")]
DimArray with dimensions:
 X: (11:10:41) m (Sampled: Ordered Regular Points)
and referenced dimensions:
 Time (type Ti): 21 s (Sampled: Ordered Regular Points)
and data: 4-element Array{Float64,1}
[0.819172, 0.418113, 0.461722, 0.379877]

For values other than Int/AbstractArray/Colon (which are set aside for regular indexing) the At selector is assumed, and can be dropped completely:

julia> A = DimArray(rand(3, 3), (X(Val((:a, :b, :c))), Y([25.6, 25.7, 25.8])));

julia> A[:b, 25.8]
0.61839141062599

Compile-time selectors

Using all Val indexes (only recommended for small arrays) you can index with named dimensions At arbitrary values with no runtime cost:

julia> A = DimArray(rand(3, 3), (cat=Val((:a, :b, :c)),
                                 val=Val((5.0, 6.0, 7.0))));

julia> @btime $A[:a, 7.0]
  2.094 ns (0 allocations: 0 bytes)
0.25620608873275397

julia> @btime $A[cat=:a, val=7.0]
  2.091 ns (0 allocations: 0 bytes)
0.25620608873275397

Methods where dims can be used containing indices or Selectors

getindex, setindex! view

Methods where dims can be used

• size, axes, firstindex, lastindex
• cat
• reverse
• dropdims
• reduce, mapreduce
• sum, prod, maximum, minimum,
• mean, median, extrema, std, var, cor, cov
• permutedims, adjoint, transpose, Transpose
• mapslices, eachslice
• fill

Warnings

Indexing with unordered or reverse order arrays has undefined behaviour. It will trash the dimension index, break searchsorted and nothing will make sense any more. So do it at you own risk. However, indexing with sorted vectors of Int can be useful. So it's allowed. But it will still do strange things to your interval sizes if the dimension span is Irregular.

Alternate Packages

There are a lot of similar Julia packages in this space. AxisArrays.jl, NamedDims.jl, NamedArrays.jl are registered alternative that each cover some of the functionality provided by DimensionalData.jl. DimensionalData.jl should be able to replicate most of their syntax and functionality.

AxisKeys.jl and AbstractIndices.jl are some other interesting developments. For more detail on why there are so many similar options and where things are headed, read this thread.