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 
    [;coords_record=:tmodel] [coords_record_units="year"]
)

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

The default for coords_record 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] [; coords::AbstractVector]) -> FieldArray
get_array(output::PALEOmodel.AbstractOutputWriter, varname::AbstractString; allselectargs...) -> FieldArray

Return a PALEOmodel.FieldArray containing data values and any attached coordinates, for the PALEOmodel.FieldRecord for varname, and records and spatial region defined by selectargs

If coords is not supplied, equivalent to PALEOmodel.get_array(PB.get_field(output, varname), allselectargs).

Optional argument coords can be used to supply plot coordinates from Variable in output, to replace any default coordinates. 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")]

NB: the coordinates will be generated by applying selectargs, so the supplied coordinate Variables must have the same dimensionality as vars.

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

Return the PALEOmodel.FieldRecord for varname.

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

Get Variable varname raw data array, optionally restricting to records

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

Return grid for output Domain domainname.

FieldRecord{D <: AbstractData, S <: AbstractSpace, V, N, M, R}
FieldRecord(
    f::PB.Field{D, S, V, N, M}, attributes; 
    coords_record, 
    sizehint::Union{Nothing, Int}=nothing
) -> fr

A series of records::R each containing the values from a PALEOboxes.Field{D, S, N, V, M}.

A coords_record may be attached to provide a coordinate (eg model time) corresponding to records.

Implementation

Fields with array values are stored in records as a Vector of arrays. Fields with single values (field_single_element true) are stored as a Vector of eltype(Field.values).

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 coords_record (usually :tmodel, when storing output vs time).

Implementation

data::DataFrame contains columns of same type as FieldRecord.records for each Variable.

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"
• add_coordinates::Bool = false: true to attempt to add CF convention coords to variables (experimental, doesn't look that useful)

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

save_jld2(output::OutputMemory, filename)

Deprecated - use save_netcdf

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

load_jld2!(output::OutputMemory, filename)

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

Example

julia> output = PALEOmodel.OutputWriters.load_jld2!(PALEOmodel.OutputWriters.OutputMemory(), "savedoutput.jld2")

FieldArray

A generic xarray-like or IRIS-like 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.

get_array(f::Field [, selectargs::NamedTuple]; [attributes=nothing]) -> FieldArray

Return a FieldArray containing f::Field data values and any attached coordinates, for the spatial region defined by selectargs.

Available selectargs depend on the grid f.mesh, and are passed to PB.Grids.get_region.

attributes (if present) are added to FieldArray

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(modeldata::AbstractModelData, vars::Union{AbstractString, Vector{<:AbstractString}} [, selectargs::NamedTuple] [; coords::AbstractVector])
heatmap(modeldata::AbstractModelData, var::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")]

NB: the coordinates will be generated by applying selectargs, so the supplied coordinate Variables must have the same dimensionality as vars.

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_records can be used to supply plot coordinates from FieldRecords. Format is a Vector of Pairs of "coordname"=>(cr1::FieldRecord, cr2::FieldRecord, ...) Example: coordsrecords=["z"=>(zmid::FieldRecord, zlower::FieldRecord, zupper::FieldRecord)] to replace a 1D column default pressure coordinate with a z coordinate. NB: the coordinates will be generated by applying selectargs, so the supplied coordinate FieldRecords must have the same dimensionality as fr.

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)