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. See ncread to automatically create a ClimArray from a .nc file.

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)

ncread(file, var; name, kwargs...) → A

Load the variable var from the file and convert it into a ClimArray with proper dimension mapping and also containing the variable attributes as a dictionary. Dimension attributes are also given to the dimensions of A, if any exist. See keywords below for specifications for unstructured grids.

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

using Glob # for getting all files
alldata = glob("toa_fluxes_2020_*.nc")
file = NCDataset(alldata; aggdim = "time")
A = ClimArray(file, "tow_sw_all")

var is a String denoting which variable to load. For .nc data containing groups var can also be a tuple ("group_name", "var_name") that loads a specific variable from a specific group. In this case, the attributes of both the group and the CF-variable are attributed to the created ClimArray.

See also ncdetails, nckeys and ncwrite.

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).

Keywords

• name optionally rename loaded array.
• grid = nothing optionally specify whether the underlying grid is grid = LonLatGrid() or grid = UnstructuredGrid(), see Types of spatial coordinates. If nothing, we try to deduce automatically based on the names of dimensions and other keys of the NCDataset.
• lon, lat. These two keywords are useful in unstructured grid data where the grid information is provided in a separate .nc file. What we need is the user to provide vectors of the central longitude and central latitude of each grid point. This is done e.g. by

  ds = NCDataset("path/to/grid.nc");
  lon = Array(ds["clon"]);
  lat = Array(ds["clat"]);

  If lon, lat are given, grid is automatically assumed UnstructuredGrid().

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.

ncwrite(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.

See also ncread.

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). This is also done in the case where the time sampling is unknown.

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 successive entries is exactly 1 hour.
• :daily, where the temporal difference between successive entries is exactly 1 day.
• :monthly, where all dates have the same day, but different month.
• :yearly, where all dates have the same month and 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)

Space coordinates are represented by two orthogonal dimensions Lon, Lat, one being longitude and the other being latitude.

Space coordinates are represented by a single dimension Coord, whose elements are coordinate locations, i.e. 2-element SVector(longitude, latitude). Each coordinate represents an equal area polygon corresponding to the point in space. The actual limits of each polygon are not included in the dimension for performance reasons.

This dimension allows indexing according to the underlying Lon, Lat representation, e.g. you can do

A # some `ClimArray` with unstructured grid type.
A[Coord(Lon(Between(0, 30)), Lat(Between(-30, 30)))]

To use functions such as zonalmean or hemispheric_means with this grid, you must first sort the ClimArray so that the latitudes of its coordinates are sorted in ascending order. I.e.

A # some `ClimArray` with unstructured grid type.
coords = dims(A, Coord).val
si = sortperm(coords, by = reverse)
A = A[Coord(si)]

This is done automatically by ncread.

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, and optionally a lookup index. 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.

Broadcasting and most Base methods maintain and sync dimension context.

DimensionalData.jl also implements:

• comprehensive plot recipes for Plots.jl.
• a Tables.jl interface with DimTable
• multi-layered DimStacks that can be indexed together, and have base methods applied to all layers.
• the Adapt.jl interface for use on GPUs, even as GPU kernel arguments.
• traits for handling a wide range of spatial data types accurately.

Dimensions

Dimensions are wrapper types. They hold the lookup index, details about the grid, and other metadata. They are also used to index into the array. X, Y, Z and Ti are the exported defaults. A generalised Dim type is available to use arbitrary symbols to name dimensions. Custom dimension types can also be defined using the @dim macro.

Dimensions can be used to construct arrays in rand, ones, zeros and fill with either a range for a lookup index or a number for the dimension length:

julia> using DimensionalData

julia> A = rand(X(1:40), Y(50))
40×50 DimArray{Float64,2} with dimensions:
  X: 1:40 (Sampled - Ordered Regular Points)
  Y
 0.929006   0.116946  0.750017  …  0.172604  0.678835   0.495294
 0.0550038  0.100739  0.427026     0.778067  0.309657   0.831754
 ⋮                              ⋱
 0.647768   0.965682  0.049315     0.220338  0.0326206  0.36705
 0.851769   0.164914  0.555637     0.771508  0.964596   0.30265

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

julia> A[Y(1), X(1:10)]
10-element DimArray{Float64,1} with dimensions:
  X: 1:10 (Sampled - Ordered Regular Points)
and reference dimensions: Y(1) 
 0.929006
 0.0550038
 0.641773
 ⋮
 0.846251
 0.506362
 0.0492866

And this has no runtime cost:

julia> A = ones(X(3), Y(3))
3×3 DimArray{Float64,2} with dimensions: X, Y
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0

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

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

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

julia> A = rand(X(40), Y(50))
40×50 DimArray{Float64,2} with dimensions: X, Y
 0.377696  0.105445  0.543156  …  0.844973  0.163758  0.849367
 ⋮                             ⋱
 0.431454  0.108927  0.137541     0.531587  0.592512  0.598927

julia> A[Y=3]
40-element DimArray{Float64,1} with dimensions: X
and reference dimensions: Y(3)
 0.543156
 ⋮
 0.137541

julia> view(A, Y(), X(1:5))
5×50 DimArray{Float64,2} with dimensions: X, Y
 0.377696  0.105445  0.543156  …  0.844973  0.163758  0.849367
 ⋮                             ⋱
 0.875279  0.133032  0.925045     0.156768  0.736917  0.444683

And for specifying dimension number in all Base and Statistics functions that have a dims argument:

julia> using Statistics

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

julia> mean(A; dims=Ti)
3×4×1 DimArray{Float64,3} with dimensions: X, Y, Ti (Time)
[:, :, 1]
 0.168058  0.52353   0.563065  0.347025
 0.472786  0.395884  0.307846  0.518926
 0.365028  0.381367  0.423553  0.369339

You can also use symbols to create Dim{X} dimensions, although we can't use the rand method directly with Symbols, and insteadd use the regular DimArray constructor:

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

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

4×20 DimArray{Float64,2} with dimensions: Dim{:a}, Dim{:b}
and reference dimensions: Dim{:c}(9)
 0.134354  0.581673  0.422615  …  0.410222   0.687915  0.753441
 0.573664  0.547341  0.835962     0.0353398  0.794341  0.490831
 0.166643  0.133217  0.879084     0.695685   0.956644  0.698638
 0.325034  0.147461  0.149673     0.560843   0.889962  0.75733

Selectors

Selectors find indices in the lookup index for each dimension:

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

(Between and Contains exlude the upper boundary so that adjacent selections never contain the same index)

Selectors can be used in getindex, setindex! and view to select indices matching the passed in value(s)

We can use selectors inside dim wrappers:

julia> using Dates

julia> timespan = DateTime(2001,1):Month(1):DateTime(2001,12)
DateTime("2001-01-01T00:00:00"):Month(1):DateTime("2001-12-01T00:00:00")

julia> A = DimArray(rand(12,10), (Ti(timespan), X(10:10:100)))
12×10 DimArray{Float64,2} with dimensions:
  Ti (Time): DateTime("2001-01-01T00:00:00"):Month(1):DateTime("2001-12-01T00:00:00") (Sampled - Ordered Regular Points)
  X: 10:10:100 (Sampled - Ordered Regular Points)
 0.14106   0.476176  0.311356  0.454908  …  0.464364  0.973193  0.535004
 ⋮                                       ⋱
 0.522759  0.390414  0.797637  0.686718     0.901123  0.704603  0.0740788

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

Without dim wrappers selectors must be in the right order:

using Unitful

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

julia> A[Between(10.5u"m", 50.5u"m"), Near(23u"s")]
4-element DimArray{Float64,1} with dimensions:
  X: (11:10:41) m (Sampled - Ordered Regular Points)
and reference dimensions:
  Ti(21 s) (Time): 21 s (Sampled - Ordered Regular Points)
 0.584028
 ⋮
 0.716715

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 = rand(X([: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 = rand(X(Val((:a, :b, :c))), Y(Val((5.0, 6.0, 7.0))))
3×3 DimArray{Float64,2} with dimensions:
  X: Val{(:a, :b, :c)}() (Categorical - Unordered)
  Y: Val{(5.0, 6.0, 7.0)}() (Categorical - Unordered)
 0.5808   0.835037  0.528461
 0.8924   0.431394  0.506915
 0.66386  0.955305  0.774132

julia> @btime $A[:c, 6.0]
  2.777 ns (0 allocations: 0 bytes)
0.9553052910459472

julia> @btime $A[Val(:c), Val(6.0)]
  1.288 ns (0 allocations: 0 bytes)
0.9553052910459472

Methods where dims can be used containing indices or Selectors

getindex, setindex! view

Methods where dims, dim types, or Symbols can be used to indicate the array dimension:

• 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

Methods where dims can be used to construct DimArrays:

• fill, ones, zeros, rand

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.