PALEOmodel output

Run

Container for model and output.

Fields

• model::Union{Nothing, PB.Model}: The model instance.
• output::Union{Nothing, AbstractOutputWriter}: model output

AbstractOutputWriter

Interface implemented by containers for PALEO model output.

Implementations should define methods for:

Writing output

• initialize!
• add_record!

Modifying output

• PB.add_field!

Querying output

• PB.get_table
• PB.show_variables
• PB.has_variable

Accessing output data

• PALEOmodel.get_array
• PB.get_field
• PB.get_mesh
• PB.get_data

initialize!(
    output::PALEOmodel.AbstractOutputWriter, model, modeldata, [nrecords] 
    [;record_dim_name=:tmodel] [record_coord_units="yr"]
)

Initialize from a PALEOboxes::Model, optionally reserving memory for an assumed output dataset of nrecords.

The default for record_dim_name is :tmodel, for a sequence of records following the time evolution of the model.

add_record!(output::PALEOmodel.AbstractOutputWriter, model, modeldata, rec_coord)

Add an output record for current state of model at record coordinate rec_coord. The usual case (set by initialize!) is that the record coordinate is model time tmodel.

add_field!(output::PALEOmodel.AbstractOutputWriter, fr::PALEOmodel.FieldRecord)

Add PALEOmodel.FieldRecord fr to output, with Domain and name defined by fr.attributes[:var_domain] and fr.attributes[:var_name].

get_table(output::PALEOmodel.AbstractOutputWriter, domainname::AbstractString) -> Table
get_table(output::PALEOmodel.AbstractOutputWriter, varnames::Vector{<:AbstractString}) -> Table

Return a DataFrame with raw model output data for Domain domainname, or for Variables varnames

show_variables(output::PALEOmodel.AbstractOutputWriter; kwargs...) -> Table
show_variables(output::PALEOmodel.AbstractOutputWriter, domainname::AbstractString; kwargs...) -> Table

Keywords:

• attributes=[:units, :vfunction, :space, :field_data, :description]: Variable attributes to include
• filter = attrb->true: function to filter by Variable attributes. Example: filter=attrb->attrb[:vfunction]!=PB.VF_Undefined to show state Variables and derivatives.

Examples:

julia> vscodedisplay(PB.show_variables(run.output))  # show all output Variables in VS code table viewer
julia> vscodedisplay(PB.show_variables(run.output, ["atm.pCO2PAL", "fluxOceanBurial.flux_P_total"]))  # show subset of output Variables in VS code table viewer

has_variable(output::PALEOmodel.AbstractOutputWriter, varname::AbstractString)  -> Bool

True if model output contains Variable varname.

get_array(output::PALEOmodel.AbstractOutputWriter, varname::AbstractString [, allselectargs::NamedTuple]; kwargs...) -> FieldArray
[deprecated] get_array(output::PALEOmodel.AbstractOutputWriter, varname::AbstractString; allselectargs...) -> FieldArray

Return a PALEOmodel.FieldArray containing data values and any attached coordinates.

Equivalent to PALEOmodel.get_array(PB.get_field(output, varname), allselectargs; kwargs), see PALEOmodel.get_array(fr::PALEOmodel.FieldRecord).

get_field(output::PALEOmodel.AbstractOutputWriter, varname::AbstractString) -> FieldRecord

Return the PALEOmodel.FieldRecord for varname.

get_data(output::PALEOmodel.AbstractOutputWriter, varname; records=nothing, kwargs...) -> values

Get Variable varname raw data array, optionally restricting to records.

Equivalent to PB.get_data(PB.get_field(output, varname); records, kwargs...), see PB.get_data(fr::PALEOmodel.FieldRecord).

PB.get_data(fr::FieldRecord; records=nothing, squeeze_all_single_dims=true)

Get data records in raw format. Only recommended for variables with scalar data ie one value per record.

records may be nothing to get all records, an Int to select a single record, or a range to select multiple records.

If squeeze_all_single_dims=true (the default), if each record represents a scalar (eg a PALEO Variable with Space PB.ScalarSpace, or a PB.CellSpace variable in a Domain with a single cell), then data is returned as a Vector. NB: if records is an Int, the single record requested is returned as a length-1 Vector.

Non-scalar data (eg a non-ScalarSpace variable from a Domain with > 1 cell) is returned in internal format as a Vector-of-Vectors.

get_mesh(output::PALEOmodel.AbstractOutputWriter, domainname::AbstractString) -> grid::Union{AbstractMesh, Nothing}

Return grid for output Domain domainname.

FieldRecord{FieldData <: AbstractData, Space <: AbstractSpace, V, N, Mesh <: AbstractMeshOrNothing, R}
FieldRecord(dataset, f::PB.Field, attributes; [sizehint=nothing])

A series of records::R each containing values from a PALEOboxes.Field{FieldData, Space, N, V, Mesh}.

Access stored values:

• As a FieldArray, using FieldArray(fr::FieldRecord) or get_array.
• For scalar data only, as a Vector using PALEOboxes.get_data.

Implementation

Storage in records::R is an internal format that may have:

• An optimisation to store Fields with single values as a Vector of eltype(Field.values), cf Fields with array values are stored in records as a Vector of arrays.
• A spatial cartesian grid stored as a linear Vector (eg mesh isa PALEOboxes.Grids.CartesianLinearGrid)

OutputMemory(; user_data=Dict{String, UserDataTypes}())

In-memory container for model output, organized by model Domains.

Implements the PALEOmodel.AbstractOutputWriter interface, with additional methods save_netcdf and load_netcdf! to save and load data.

Implementation

• Field domains::Dict{String, OutputMemoryDomain} contains per-Domain model output.
• Field user_data::Dict{String, UserDataTypes} contains optional user data NB:
  • available types are restricted to those that are compatible with NetCDF attribute types, ie Float64, Int64, String, Vector{Float64}, Vector{Int64}, Vector{String}
  • Vectors with a single element are read back from netcdf as scalars, see https://alexander-barth.github.io/NCDatasets.jl/dev/issues/#Corner-cases

OutputMemoryDomain

In-memory model output, for one model Domain.

Includes an additional record coordinate variable not present in Domain (usually :tmodel, when storing output vs time).

save_netcdf(output::OutputMemory, filename; kwargs...)

Save to filename in netcdf4 format (NB: filename must either have no extension or have extension .nc)

Notes on structure of netcdf output

• Each PALEO Domain is written to a netcdf4 group. These can be read into a Python xarray using the group=<domainname> argument to open_dataset.
• Isotope-valued variables (field_data = PB.IsotopeLinear) are written with an extra isotopelinear netCDF dimension, containing the variable total and delta.
• Any '/' characters in PALEO variables are substited for '%' in the netcdf name.

Keyword arguments

• check_ext::Bool = true: check that filename ends in ".nc"

load_netcdf!(output::OutputMemory, filename)

Load from filename in netCDF format, replacing any existing content in output. (NB: filename must either have no extension or have extension .nc).

Example

julia> output = PALEOmodel.OutputWriters.load_netcdf!(PALEOmodel.OutputWriters.OutputMemory(), "savedoutput.nc")

FieldArray

A generic xarray-like or IRIS-like n-dimensional Array with named dimensions and optional coordinates.

NB: this aims to be simple and generic, not efficient !!! Intended for representing model output, not for numerically-intensive calculations.

Fields

• name::String: variable name

• values::AbstractArray: n-dimensional Array of values

• dims_coords::Vector{Pair{PALEOboxes.NamedDimension, Vector{PALEOmodel.FieldArray}}}: Names of dimensions with optional attached coordinates

• attributes::Union{Nothing, Dict{Symbol, Any}}: variable attributes

FieldArray(
    fr::FieldRecord; 
    expand_cartesian=true, 
    omit_recorddim_if_constant=true,
    lookup_coords=true,
    coords=nothing,
)

Construct FieldArray from all records in fr::FieldRecord

Provides a view of internal storage format of FieldRecord as an n-dimensional Array.

Keyword arguments

• expand_cartesian: (spatially resolved output using a CartesianLinearGrid only), true to expand compressed internal data (with spatial dimension cells) to a cartesian array (with spatial dimensions eg lon, lat, zt)

• omit_recorddim_if_constant: Specify whether to include multiple (identical) records and record dimension for constant variables (with attribute is_constant = true). PALEO currently always stores these records in the input fr::FieldRecord; set false include them in FieldArray output, true to omit duplicate records and record dimension from output.

• lookup_coords: true to include coordinates, false to omit coordinates.

• coords: replace the attached coordinates from the input fr::FieldRecord for one or more dimensions. Format is a Vector of Pairs of "<dim_name>"=>("<var_name1>", "<var_name2>", ...), eg to replace an nD column atmosphere model default pressure coordinate with a z coordinate:

  coords=["cells"=>("atm.zmid", "atm.zlower", "atm.zupper")]

select(fa::FieldArray [, allselectargs::NamedTuple] ; kwargs...) -> fa_out::FieldArray

Select a region of a FieldArray defined by allselectargs.

Selecting records and regions

allselectargs is a NamedTuple of form:

(<dimcoordname> = <filter>, <dimcoordname> = <filter>,  ... [, squeeze_all_single_dims=true])

where <dimcoordname> is of form:

• <dimname>_isel to select by array indices: <filter> may then be a single Int to select a single index, or a range first:last to select a range of indices.
• <coordname> to select by coordinate values using the coordinates attached to each dimension: <filter> may then be a single number to select a single index corresponding to the nearest value of the corresponding coordinate, or (first::Float64, last::Float64) (a Tuple) to select a range starting at the index with the nearest value of fr.coords_record before first and ending at the nearest index after last.

Available dimensions and coordinates <dimcoordname> depend on the FieldArray dimensions (as returned by get_dimensions, which will be a subset of grid spatial dimensions and Domain data dimensions) and corresponding attached coordinates (as returned by get_coordinates).

Some synonyms are defined for commonly used <dimnamecoordname>, these may require optional record_dim_name and mesh to be supplied in kwargs:

NB: Dimensions corresponding to a selection for a single index or coordinate value are always squeezed out from the returned FieldArray. Optional argument squeeze_all_single_dims (default true) controls whether all output dimensions that contain a single index are squeezed out (including eg a selection for a range that results in a dimension with one index, or where the input FieldArray contains a dimension with a single index).

Selection arguments used are optionally returned as strings in fa_out.attributes[:filter_records] and fa_out.attributes[:filter_region] for use in plot labelling etc.

Keyword arguments

• record_dim_name=nothing: optionally supply name of record dimension as a String, to substitute for <dimcoordname> records and to define dimension to use to populate filter_records attribute.
• mesh=nothing: optionally supply a grid from PALEOboxes.Grids, to use for substituting cell names, looking up column indices.
• add_attributes=true: true to transfer attributes from input fr::FieldRecord to output FieldArray adding :filter_region and :filter_records, false to omit.
• update_name=true: true to update output FieldArray name to add Domain name prefix, and a suffix generated from allselectargs (NB: requires add_attributes=true), false to use name from input FieldRecord.

Limitations

• it is only possible to select either a single slice or a contiguous range for each dimension, not a set of slices for a Vector of index or coordinate values.

get_array(
    fr::FieldRecord [, allselectargs::NamedTuple];
    coords=nothing,
    lookup_coords=true,
    add_attributes=true,
    update_name=true,
    omit_recorddim_if_constant=false, 
) -> fa_out::FieldArray

[deprecated] get_array(fr::FieldRecord [; kwargs...] [; allselectargs...]) -> fa_out::FieldArray

Return a FieldArray containing an Array of values and any attached coordinates, for records and spatial region defined by allselectargs.

Combines FieldArray(fr::FieldRecord) and select(fa::FieldArray)

Optional allselectargs field expand_cartesian (default false) is only needed for spatially resolved output using a CartesianLinearGrid, set to true to expand compressed internal data (with spatial dimension cells) to a cartesian array (with spatial dimensions eg lon, lat, zt)

Keyword arguments

See FieldArray(fr::FieldRecord) and select(fa::FieldArray)

Examples

• select a timeseries for single cell index 3 :

  get_array(fr, (cell=3, ))

• select first column at nearest available time to model time 10.0 :

  get_array(fr, (column=1, tmodel=10.0))

• set first column at nearest available time to model time 10.0, replacing atmosphere model pressure coordinate with z coordinate:

  get_array(
      fr, (column=1, tmodel=10.0);
      coords=["cells"=>("atm.zmid", "atm.zlower", "atm.zupper")]
  )

• select surface 2D array (dimension zt, index 1) from 3D output at nearest available time to model time 10.0, expanding CartesianLinearGrid:

  get_array(fr, (zt_isel=1, tmodel=10.0, expand_cartesian=true))

• select section 2D array at nearest index to longitude 200 degrees from 3D output at nearest available time to model time 10.0, expanding CartesianLinearGrid:

  get_array(fr, (lon=200.0, tmodel=10.0, expand_cartesian=true))

• get full data cube as used for netcdf output, omitting coordinates and attributes, retaining singleton dimensions, and omitting record dimension if variable is a constant:

  get_array(
      fr, (expand_cartesian=true, squeeze_all_single_dims=false);
      lookup_coordinates=false, add_attributes=false, omit_recorddim_if_constant=true
  )

Limitations

• it is only possible to select either a single slice or a contiguous range for each dimension, not a set of slices for a Vector of index or coordinate values.

plot(output::AbstractOutputWriter, vars::Union{AbstractString, Vector{<:AbstractString}} [, selectargs::NamedTuple] [; coords::AbstractVector])
heatmap(output::AbstractOutputWriter, var::AbstractString [, selectargs::NamedTuple] [; coords::AbstractVector])
plot(outputs::Vector{<:AbstractOutputWriter}, vars::Union{AbstractString, Vector{<:AbstractString}} [, selectargs::NamedTuple] [; coords::AbstractVector])

Plot recipe that calls PB.get_field(output, var), and passes on to plot(fr::FieldRecord, selectargs) (see RecipesBase.apply_recipe(::Dict{Symbol, Any}, fr::FieldRecord, selectargs::NamedTuple))

Vector-valued outputs or vars are "broadcast" to create a plot series for each element. A labelprefix (index in outputs Vector) is added to identify each plot series produced.

If var is of form <domain>.<name>.<structfield>, then sets the structfield keyword argument to take a single field from a struct Variable.

Optional argument coords can be used to supply plot coordinates from Variable in output. Format is a Vector of Pairs of "coordname"=>("varname1", "var_name2", ...)

Example: to replace a 1D column default pressure coordinate with a z coordinate:

coords=["z"=>("atm.zmid", "atm.zlower", "atm.zupper")]

plot(fr::FieldRecord, selectargs::NamedTuple)
heatmap(fr::FieldRecord, selectargs::NamedTuple)

Plot recipe to plot a single FieldRecord

Calls get_array(fr, selectargs) and passes on to plot(fa::FieldArray) (see RecipesBase.apply_recipe(::Dict{Symbol, Any}, fa::FieldArray)).

Vector-valued fields in selectargs are "broadcasted" (generating a separate plot series for each combination)

Optional argument coords can be used to supply plot coordinates from Variable in output. Format is a Vector of Pairs of "coordname"=>("varname1", "var_name2", ...)

Example: to replace a 1D column default pressure coordinate with a z coordinate:

coords=["z"=>("atm.zmid", "atm.zlower", "atm.zupper")]

plot(fa::FieldArray; kwargs...)
heatmap(fa::FieldArray; kwargs...)
plot(fas::Vector{<:FieldArray}; kwargs...)

Plot recipe that plots a single [FieldArray] or Vector of [FieldArray].

If fa has a single dimension, this is suitable for a line-like plot, if two dimensions, a heatmap.

If fas::Vector is supplied, this is "broadcast" generating one plot series for each element, with the Vector index prepended to labelprefix to identify the plot series (unless overridden by labellist or labelattribute)

Keywords

• swap_xy::Bool=false: true to swap x and y axes
• mult_y_coord=1.0: workaround for bugs in Plots.jl heatmap yflip - multiply y coordinate by constant factor.
• structfield::Union{Symbol, Nothing}=nothing: use field structfield from a struct-valued array.
• map_values=PB.get_total: function to apply to y (for a 1D series) or z (for a 2D heatmap etc) before plotting
• labelprefix="": prefix for plot label.
• labellist=[]: list of labels to override defaults
• labelattribute=nothing: FieldArray attribute to use as label

PlotPager(layout [, kwargs=NamedTuple()][; displayfunc=(plot, nplot)->display(plot)])

Accumulate plots into subplots.

layout is supplied to Plots.jl layout keyword, may be an Int or a Tuple (ny, nx), see https://docs.juliaplots.org/latest/.

Optional argument kwargs::NamedTuple provides additional keyword arguments passed through to plot (eg (legend_background_color=nothing, ) to set all subplot legends to transparent backgrounds)

Optional keyword argument displayfunc::(plot, nplot)->display(plot) provides the function used to display a screen of plots, where plot is a Plot object and nplot::Integer is the sequential number of this screen.

Usage

julia> pp = PlotPager((2,2))  # 4 panels per screen (2 down, 2 across)
julia> pp(plot(1:3))  # Accumulate
julia> pp(:skip, plot(1:4), plot(1:5), plot(1:6))  # add multiple panels in one command
julia> pp(:newpage) # flush any partial screen and start new page (NB: always add this at end of a sequence!)

julia> pp = PlotPager((2, 2); displayfunc=(plot, nplot)->savefig(plot, "plot_$nplot.png")) # save to file instead of default display

Commands

• pp(p): accumulate plot p
• pp(:skip): leave a blank panel
• pp(:newpage): fill with blank panels and start new page
• pp(p1, p2, ...): multiple plots/commands in one call

DefaultPlotPager()

Dummy version of PlotPager that just calls display for each plot.

ReactionNetwork

Functions to analyze a PALEOboxes.Model or PALEOmodel output that contains a reaction network

Compiles reaction stoichiometry and rate information from attributes attached to reaction rate variables:

• rate_processname::String: a process name (eg "photolysis", "reaction", ...)
• rate_species::Vector{String} names of reactant and product species
• rate_stoichiometry::Vector{Float64} stoichiometry of reactant and product species

get_ratetable(model, domainname) -> DataFrame
get_ratetable(output, domainname) -> DataFrame

Get table of rate Variables and reactions

Returns a DataFrame with columns :name, :rate_processname, :rate_species, :rate_stoichiometry

get_all_species_ratevars(model, domainname) -> OrderedDict(speciesname => [(stoich, ratevarname, processname), ...])
get_all_species_ratevars(output, domainname) -> OrderedDict(speciesname => [(stoich, ratevarname, processname), ...])
get_all_species_ratevars(ratetable::DataFrame) -> OrderedDict(speciesname => [(stoich, ratevarname, processname), ...])

Get all species and contributing reaction rate Variable names as Dict of Tuples (stoich, ratevarname, processname) where ratevarname is the name of an output Variable with a reaction rate, stoich is the stoichiometry of that rate applied to species, and processname identifies the nature of the reaction.

get_rates(output, domainname [, outputrec] [, indices] [, scalefac] [, add_equations] [, ratetable_source]) -> DataFrame

Get all reaction rates as column rate_total for domainname from output record outputrec (defaults to last time record), for subset of cells in indices (defaults to whole domain).

Set optional ratetable_source = model to use with older output that doesn't include rate variable attributes.

get_all_species_ratesummaries(output, domainname [, outputrec] [, indices] [, scalefac] [, ratetable_source]) 
    -> OrderedDict(speciesname => (source, sink, net, source_rxs, sink_rxs))

Get source, sink, net rates and rates of source_rxs and sink_rxs for all species in domainname from output record outputrec (defaults to last record), cells in indices (defaults to whole domain),

Optional scalefac to convert units, eg scalefac=1.90834e12 to convert mol m-2 yr-1 to molecule cm-2 s-1

Set optional ratetable_source = model to use with older output that doesn't include rate variable attributes.

show_ratesummaries([io::IO = stdout], ratesummaries [; select_species=[]])

Print per-species reaction rates from ratesummaries to output stream io (defaults to stdout), optionally selecting species to print.

Example

ratesummaries = PALEOmodel.ReactionNetwork.get_all_species_ratesummaries(output, "atm")
PALEOmodel.ReactionNetwork.show_ratesummaries(ratesummaries)