From 1e4d9e51393f77c704d52366a87066fa4c7d0bed Mon Sep 17 00:00:00 2001 From: Stuart Daines Date: Fri, 29 Jul 2022 20:27:31 +0100 Subject: [PATCH 1/4] get_array can now modify / replace coordinates Main use is to support time-varying coordinates (default coordinates can only be time-independent). 'get_array(output, ...)' and 'plot(output, ...)' can now take an optional argument `coords` to supply plot coordinates from Variables in output, to replace any default coordinates. Format is a Vector of Pairs of "coord_name"=>("var_name1", "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`. Also rationalize 'get_array' syntax to match that of 'plot' so that we have: get_array(output, "varname", (tmodel=1.0, column=1)) and plot(output, "varname", (tmodel=1.0, column=1)) --- src/FieldRecord.jl | 63 +++++++++++++++++++++++++++++++++++++++----- src/OutputWriters.jl | 45 ++++++++++++++++++++++++++----- src/PlotRecipes.jl | 51 ++++++++++++++++++++++++++++------- 3 files changed, 137 insertions(+), 22 deletions(-) diff --git a/src/FieldRecord.jl b/src/FieldRecord.jl index f6e9d25..b61c029 100644 --- a/src/FieldRecord.jl +++ b/src/FieldRecord.jl @@ -128,12 +128,14 @@ function Base.copy(fr::FieldRecord{D, S, V, N, M, R}) where {D, S, V, N, M, R} ) end -""" - get_array(fr::FieldRecord; [recordarg] [,selectargs...]) -> FieldArray +""" + get_array(fr::FieldRecord [, allselectargs::NamedTuple] [,coords_records::AbstractVector]) -> FieldArray + get_array(fr::FieldRecord; [allselectargs...]) -> FieldArray Return a [`FieldArray`](@ref) containing `fr::FieldRecord` data values and -any attached coordinates, for records defined by `recordarg` and the -spatial region defined by `selectargs`. +any attached coordinates, for records and spatial region defined by `allselectargs`. + +`allselectarg` may include `recordarg` to define records, `selectargs` to define a spatial region. `recordarg` can be one of: - `records=r::Int` to select a single record, or `records = first:last` to select a range. @@ -144,16 +146,37 @@ spatial region defined by `selectargs`. Available `selectargs` depend on the grid `fr.mesh`, and are passed to `PB.Grids.get_region`. + + +Optional argument `coords_records` can be used to supply coordinates from additional `FieldRecords`, replacing any coordinates +attached to `fr`. Format is a Vector of Pairs of "coord_name"=>(cr1::FieldRecord, cr2::FieldRecord, ...) + +Example: to replace a 1D column default pressure coordinate with a z coordinate: + + coords_records=["z"=>(zmid::FieldRecord, zlower::FieldRecord, zupper::FieldRecord)] + +NB: the coordinates will be generated by applying `allselectargs`, +so the supplied coordinate FieldRecords must have the same dimensionality as `fr`. """ function get_array( - fr::FieldRecord{D, S, V, N, M, R}; selectargs... + fr::FieldRecord; allselectargs... +) + Base.depwarn( + "allselectargs... will be deprecated in a future release. Please use allselectargs::NamedTuple instead", + :get_array, + ) + return get_array(fr, NamedTuple(allselectargs)) +end + +function get_array( + fr::FieldRecord{D, S, V, N, M, R}, allselectargs::NamedTuple=NamedTuple(), ) where {D, S, V, N, M, R} # select records to use and create PB.NamedDimension ridx = 1:length(fr) selectargs_region = Dict() selectargs_records = NamedTuple() - for (k, v) in selectargs + for (k, v) in zip(keys(allselectargs), allselectargs) if k==:records if v isa Integer ridx = [v] @@ -256,3 +279,31 @@ function get_array( end +# substitute or add coords_records to FieldArray(fr; selectargs...) +function get_array(fr::FieldRecord, allselectargs::NamedTuple, coords_records::AbstractVector) + + # get data NB: with original coords or no coords + varray = get_array(fr; allselectargs...) + + # generate Vector of NamedDimensions to use as new coordinates + named_dimensions = PB.NamedDimension[] + for (coord_name, coord_fieldrecords) in coords_records + fixed_coords = [] + for coord_fr in coord_fieldrecords + coord_array = get_array(coord_fr; allselectargs...) + push!(fixed_coords, PB.FixedCoord(get(coord_array.attributes, :var_name, ""), coord_array.values, coord_array.attributes)) + end + push!( + named_dimensions, PB.NamedDimension( + coord_name, + length(first(fixed_coords).values), + fixed_coords, + ) + ) + end + + # replace coordinates + varray_coords = FieldArray(varray.name, varray.values, Tuple(named_dimensions), varray.attributes) + + return varray_coords +end \ No newline at end of file diff --git a/src/OutputWriters.jl b/src/OutputWriters.jl index 134c43a..7c57d2f 100644 --- a/src/OutputWriters.jl +++ b/src/OutputWriters.jl @@ -13,7 +13,7 @@ import DataFrames import FileIO import JLD2 -# using Infiltrator # Julia debugger +import Infiltrator # Julia debugger ################################## @@ -110,18 +110,51 @@ function PB.has_variable(output::PALEOmodel.AbstractOutputWriter, varname::Abstr function PB.show_variables(output::PALEOmodel.AbstractOutputWriter) end """ - get_array(output::PALEOmodel.AbstractOutputWriter, varname::AbstractString; kwargs...) -> FieldArray + 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`](@ref) containing data values and any attached coordinates, for the -[`PALEOmodel.FieldRecord`](@ref) for `varname`, and records and spatial region defined by `kwargs` +[`PALEOmodel.FieldRecord`](@ref) for `varname`, and records and spatial region defined by `selectargs` -Equivalent to `PALEOmodel.get_array(PB.get_field(output, varname), kwargs...)` +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 "coord_name"=>("var_name1", "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`. """ -function PALEOmodel.get_array(output::PALEOmodel.AbstractOutputWriter, varname::AbstractString; kwargs...) +function PALEOmodel.get_array( + output::PALEOmodel.AbstractOutputWriter, varname::AbstractString; + allselectargs... +) + Base.depwarn( + "allselectargs... will be deprecated in a future release. Please use allselectargs::NamedTuple instead", + PALEOmodel.get_array + ) + return PALEOmodel.get_array(output, varname, NamedTuple(allselectargs)) +end +function PALEOmodel.get_array( + output::PALEOmodel.AbstractOutputWriter, varname::AbstractString, allselectargs::NamedTuple=NamedTuple(); + coords::AbstractVector=[], +) fr = PB.get_field(output, varname) - fa = PALEOmodel.get_array(fr; kwargs...) + coords_records = [ + coord_name => (PB.get_field(output, cvn) for cvn in coord_varnames) + for (coord_name, coord_varnames) in coords + ] + + if isempty(coords_records) + fa = PALEOmodel.get_array(fr, allselectargs) + else + fa = PALEOmodel.get_array(fr, allselectargs, coords_records) + end return fa end diff --git a/src/PlotRecipes.jl b/src/PlotRecipes.jl index a2ee3d6..4621e20 100644 --- a/src/PlotRecipes.jl +++ b/src/PlotRecipes.jl @@ -115,9 +115,9 @@ end ############################### """ - plot(output::AbstractOutputWriter, vars::Union{AbstractString, Vector{<:AbstractString}}, selectargs::NamedTuple=NamedTuple()) - heatmap(output::AbstractOutputWriter, vars::Union{AbstractString, Vector{<:AbstractString}}, selectargs::NamedTuple=NamedTuple()) - plot(outputs::Vector{<:AbstractOutputWriter}, vars::Union{AbstractString, Vector{<:AbstractString}}, selectargs::NamedTuple=NamedTuple()) + plot(output::AbstractOutputWriter, vars::Union{AbstractString, Vector{<:AbstractString}} [, selectargs::NamedTuple] [; coords::AbstractVector]) + heatmap(output::AbstractOutputWriter, vars::Union{AbstractString, Vector{<: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)` @@ -128,17 +128,34 @@ A `labelprefix` (index in `outputs` Vector) is added to identify each plot serie If `var` is of form `..`, 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 "coord_name"=>("var_name1", "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`. """ RecipesBase.@recipe function f( output::AbstractOutputWriter, vars::Union{AbstractString, Vector{<:AbstractString}}, - selectargs::NamedTuple=NamedTuple() + selectargs::NamedTuple=NamedTuple(); + coords::AbstractVector = [], ) if isa(vars, AbstractString) vars = [vars] end + coords_records = [ + coord_name => (PB.get_field(output, cvn) for cvn in coord_varnames) + for (coord_name, coord_varnames) in coords + ] + delete!(plotattributes, :coords) + for var in vars RecipesBase.@series begin varsplit = split(var, ".") @@ -148,7 +165,7 @@ RecipesBase.@recipe function f( end var = join(varsplit, ".") - PB.get_field(output, var), selectargs + PB.get_field(output, var), selectargs, coords_records end end @@ -165,14 +182,16 @@ adding a `labelprefix` (index in `outputs` Vector) to identify each plot series RecipesBase.@recipe function f( outputs::Vector{<:AbstractOutputWriter}, vars::Union{AbstractString, Vector{<:AbstractString}}, - selectargs::NamedTuple=NamedTuple() + selectargs::NamedTuple=NamedTuple(); + coords::AbstractVector=[], ) for (i, output) in enumerate(outputs) RecipesBase.@series begin labelprefix --> "$i: " - output, vars, selectargs + output, vars, selectargs, coords end end + delete!(plotattributes, :coords) return nothing end @@ -184,19 +203,31 @@ end Plot recipe to plot a single [`FieldRecord`](@ref) -Calls `get_array(fr; selectargs...)` and passes on to `plot(fa::FieldArray)` +Calls [`get_array(fr, selectargs)`](@ref) and passes on to `plot(fa::FieldArray)` (see [`RecipesBase.apply_recipe(::Dict{Symbol, Any}, fa::FieldArray)`](@ref)). 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 "coord_name"=>(cr1::FieldRecord, cr2::FieldRecord, ...) +Example: + coords=["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`. """ -RecipesBase.@recipe function f(fr::FieldRecord, selectargs::NamedTuple) +RecipesBase.@recipe function f(fr::FieldRecord, selectargs::NamedTuple, coords_records::AbstractVector=[]) # broadcast any Vector-valued argument in selectargs bcastargs = broadcast_dict([Dict{Symbol, Any}(pairs(selectargs))]) for sa in bcastargs + sa_nt = NamedTuple(sa) RecipesBase.@series begin - get_array(fr; sa...) + if isempty(coords_records) + get_array(fr, sa_nt) + else + get_array(fr, sa_nt, coords_records) + end end end end From 7d4da818fd2977d8cf06037cad602f8a04d8b328 Mon Sep 17 00:00:00 2001 From: Stuart Daines Date: Fri, 29 Jul 2022 20:36:39 +0100 Subject: [PATCH 2/4] update Project.toml for v0.15.4 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 3696cad..655fed2 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "PALEOmodel" uuid = "bf7b4fbe-ccb1-42c5-83c2-e6e9378b660c" authors = ["Stuart Daines "] -version = "0.15.3" +version = "0.15.4" [deps] DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" From de48b44503269fe99a67cb6e57da91213f2a41a5 Mon Sep 17 00:00:00 2001 From: Stuart Daines Date: Sun, 31 Jul 2022 16:15:47 +0100 Subject: [PATCH 3/4] bugfixes, update documentation --- docs/src/HOWTOshowmodelandoutput.md | 24 +++-- src/FieldArray.jl | 138 ++++++++++++++++++++++++---- src/FieldRecord.jl | 99 ++++++++++---------- src/OutputWriters.jl | 39 ++++---- src/PlotRecipes.jl | 95 ++++++++++++++----- 5 files changed, 277 insertions(+), 118 deletions(-) diff --git a/docs/src/HOWTOshowmodelandoutput.md b/docs/src/HOWTOshowmodelandoutput.md index 97c2757..be724ea 100644 --- a/docs/src/HOWTOshowmodelandoutput.md +++ b/docs/src/HOWTOshowmodelandoutput.md @@ -120,7 +120,7 @@ The output can be plotted using the Julia Plots.jl package, see [Plotting output ## Spatial or wavelength-dependent output -To analyze spatial or eg wavelength-dependent output (eg time series from a 1D column or 3D general circulation model, or quantities that are a function of wavelength or frequency), [`PALEOmodel.get_array`](@ref) takes additional arguments to take 1D or 2D slices from the spatial, spectral and timeseries data. The [`PALEOmodel.FieldArray`](@ref) returned includes coordinates to plot column (1D) and heatmap (2D) data. +To analyze spatial or eg wavelength-dependent output (eg time series from a 1D column or 3D general circulation model, or quantities that are a function of wavelength or frequency), [`PALEOmodel.get_array`](@ref) takes an additional `selectargs::NamedTuple` argument to take 1D or 2D slices from the spatial, spectral and timeseries data. The [`PALEOmodel.FieldArray`](@ref) returned includes default coordinates to plot column (1D) and heatmap (2D) data, these can be overridden by supplying the optional `coords` keyword argument. ### Examples for a column-based model @@ -128,18 +128,28 @@ Visualisation of spatial and wavelength-dependent output from the PALEOdev.jl oz #### 1D column data julia> plot(title="O3 mixing ratio", output, "atm.O3_mr", (tmodel=[0.0, 0.1, 1.0, 10.0, 100.0, 1000.0], column=1), - swap_xy=true, xaxis=:log, labelattribute=:filter_records) # plots O3 vs height + swap_xy=true, xscale=:log10, labelattribute=:filter_records) # plots O3 vs default height coordinate -Here the `labelattribute=:filter_records` keyword argument is used to generate plot labels from the `:filter_records` FieldArray attribute, which contains the `tmodel` values used to select the timeseries records. The plot recipe expands +Here the optional `labelattribute=:filter_records` keyword argument is used to generate plot labels from the `:filter_records` FieldArray attribute, which contains the `tmodel` values used to select the timeseries records. The plot recipe expands the Vector-valued `tmodel` argument to overlay a sequence of plots. This is equivalent to first creating and then plotting a sequence of `FieldArray` objects: - - julia> O3_mr = PALEOmodel.get_array(run.output, "atm.O3_mr", tmodel=0.0, column=1) - julia> plot(title="O3 mixing ratio", O3_mr, swap_xy=true, xaxis=:log, labelattribute=:filter_records) - julia> O3_mr = PALEOmodel.get_array(run.output, "atm.O3_mr", tmodel=0.1, column=1) + julia> O3_mr = PALEOmodel.get_array(run.output, "atm.O3_mr", (tmodel=0.0, column=1)) + julia> plot(title="O3 mixing ratio", O3_mr, swap_xy=true, xscale=:log10, labelattribute=:filter_records) + julia> O3_mr = PALEOmodel.get_array(run.output, "atm.O3_mr", (tmodel=0.1, column=1)) julia> plot!(O3_mr, swap_xy=true, labelattribute=:filter_records) +The default height coordinate from the model grid can be replaced using the optional `coords` keyword argument, eg + julia> plot(title="O3 mixing ratio", output, "atm.O3_mr", (tmodel=[0.0, 0.1, 1.0, 10.0, 100.0, 1000.0], column=1), + coords=["p"=>("atm.pmid", "atm.plower", "atm.pupper")], + swap_xy=true, xscale=:log10, yflip=true, yscale=:log10, labelattribute=:filter_records) # plots O3 vs pressure + +The current values in the `modeldata` struct can also be plotted using the same syntax, eg + + julia> plot(title="O2, O3 mixing ratios", modeldata, ["atm.O2_mr", "atm.O3_mr"], (column=1,), + swap_xy=true, xscale=:log10) # plot current value of O2, O3 vs height + + #### Wavelength-dependent data julia> plot(title="direct transmittance", output, ["atm.direct_trans"], (tmodel=1e12, column=1, cell=[1, 80]), ylabel="fraction", labelattribute=:filter_region) # plots vs wavelength diff --git a/src/FieldArray.jl b/src/FieldArray.jl index 0f9d290..7d42968 100644 --- a/src/FieldArray.jl +++ b/src/FieldArray.jl @@ -35,6 +35,23 @@ function Base.:*(fa_in::FieldArray, a::Real) end Base.:*(a::Real, fa_in::FieldArray) = fa_in*a +# default name from attributes +default_fieldarray_name(attributes::Nothing) = "" + +function default_fieldarray_name(attributes::Dict) + name = get(attributes, :domain_name, "") + name *= isempty(name) ? "" : "." + name *= get(attributes, :var_name, "") + + selectargs_records = get(attributes, :filter_records, NamedTuple()) + selectargs_region = get(attributes, :filter_region, NamedTuple()) + if !isempty(selectargs_region) || !isempty(selectargs_records) + name *= "(" * join(["$k=$v" for (k, v) in Dict(pairs(merge(selectargs_records, selectargs_region)))], ", ") * ")" + end + + return name +end + ############################################################# # Create from PALEO objects ############################################################# @@ -47,7 +64,7 @@ Get FieldArray from PALEO object `obj` function get_array end """ - get_array(f::Field; [attributes=nothing], [selectargs...]) -> FieldArray + get_array(f::Field [, selectargs::NamedTuple]; [attributes=nothing]) -> FieldArray Return a [`FieldArray`](@ref) containing `f::Field` data values and any attached coordinates, for the spatial region defined by `selectargs`. @@ -62,26 +79,29 @@ function get_array( attributes=nothing, ) where {D, V, N, M} - return FieldArray("", f.values, f.data_dims, attributes) + return FieldArray(default_fieldarray_name(attributes), f.values, f.data_dims, attributes) end function get_array( - f::PB.Field{D, PB.CellSpace, V, 0, M}; - attributes=nothing, - selectargs... + f::PB.Field{D, PB.CellSpace, V, 0, M}, selectargs::NamedTuple=NamedTuple(); + attributes=nothing, ) where {D, V, M} values, dims = PB.Grids.get_region(f.mesh, f.values; selectargs...) - return FieldArray("", values, dims, attributes) + if !isempty(selectargs) + attributes = isnothing(attributes) ? Dict{Symbol, Any}() : copy(attributes) + attributes[:filter_region] = selectargs + end + + return FieldArray(default_fieldarray_name(attributes), values, dims, attributes) end # single data dimension # TODO generalize this to arbitrary data dimensions function get_array( - f::PB.Field{D, PB.CellSpace, V, 1, M}; + f::PB.Field{D, PB.CellSpace, V, 1, M}, selectargs::NamedTuple=NamedTuple(); attributes=nothing, - selectargs... ) where {D, V, M} f_space_dims_colons = ntuple(i->Colon(), ndims(f.values) - 1) @@ -105,23 +125,107 @@ function get_array( end end - return FieldArray("", values, (dims..., f.data_dims...), attributes) + if !isempty(selectargs) + attributes = isnothing(attributes) ? Dict{Symbol, Any}() : copy(attributes) + attributes[:filter_region] = selectargs + end + + return FieldArray(default_fieldarray_name(attributes), values, (dims..., f.data_dims...), attributes) end -""" - get_array(model::PB.Model, modeldata, varnamefull; selectargs...) -> FieldArray +""" + get_array(modeldata, varnamefull [, selectargs::NamedTuple] [; coords::AbstractVector]) -> FieldArray + Get [`FieldArray`](@ref) by Variable name, for spatial region defined by `selectargs` (which are passed to `PB.Grids.get_region`). + +Optional argument `coords` can be used to supply plot coordinates from Variable in output. +Format is a Vector of Pairs of "coord_name"=>("var_name1", "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`. """ -function get_array(model::PB.Model, modeldata::PB.AbstractModelData, varnamefull::AbstractString; selectargs...) - var = PB.get_variable(model, varnamefull) +function get_array( + modeldata::PB.AbstractModelData, varnamefull::AbstractString, selectargs::NamedTuple=NamedTuple(); + coords=nothing, +) + var = PB.get_variable(modeldata.model, varnamefull) !isnothing(var) || throw(ArgumentError("Variable $varnamefull not found")) f = PB.get_field(var, modeldata) - attrbs = deepcopy(var.attributes) - attrbs[:var_name] = var.name - attrbs[:domain_name] = var.domain.name - return get_array(f; attributes=attrbs, selectargs...) + attributes = copy(var.attributes) + attributes[:var_name] = var.name + attributes[:domain_name] = var.domain.name + + varray = get_array(f, selectargs; attributes=attributes) + + if isnothing(coords) + # keep original coords (if any) + else + check_coords_argument(coords) || + error("argument coords should be a Vector of Pairs of \"coord_name\"=>(\"var_name1\", \"var_name2\", ...), eg: [\"z\"=>(\"atm.zmid\", \"atm.zlower\", \"atm.zupper\"), ...]") + + vec_coords_arrays = [ + coord_name => Tuple(get_array(modeldata, cvn, selectargs) for cvn in coord_varnames) + for (coord_name, coord_varnames) in coords + ] + + varray = update_coordinates(varray, vec_coords_arrays) + end + + return varray +end + +# check 'coords' of form [] or ["z"=>[ ... ], ] or ["z"=>(...),] +check_coords_argument(coords) = + isa(coords, AbstractVector) && ( + isempty(coords) || ( + isa(coords, AbstractVector{<:Pair}) && + isa(first(first(coords)), AbstractString) && + isa(last(first(coords)), Union{AbstractVector, Tuple}) + ) + ) + +""" + update_coordinates(varray::FieldArray, vec_coord_arrays::AbstractVector) -> FieldArray + +Replace or add coordinates `vec_coord_arrays` to `varray`. + +`new_coord_arrays` is a Vector of Pairs of "coord_name"=>(var1::FieldArray, var2::FieldArray, ...) + +Example: to replace a 1D column default pressure coordinate with a z coordinate: + + coords=["z"=>(zmid::FieldArray, zlower::FieldArray, atm.zupper::FieldArray)] +""" +function update_coordinates(varray::FieldArray, vec_coord_arrays::AbstractVector) + + check_coords_argument(vec_coord_arrays) || + error("argument vec_coords_arrays should be a Vector of Pairs of \"coord_name\"=>(var1::FieldArray, var2::FieldArray, ...), eg: [\"z\"=>(zmid::FieldArray, zlower::FieldArray, atm.zupper::FieldArray)]") + + # generate Vector of NamedDimensions to use as new coordinates + named_dimensions = PB.NamedDimension[] + for (coord_name, coord_arrays) in vec_coord_arrays + fixed_coords = [] + for coord_array in coord_arrays + push!(fixed_coords, PB.FixedCoord(get(coord_array.attributes, :var_name, ""), coord_array.values, coord_array.attributes)) + end + push!( + named_dimensions, PB.NamedDimension( + coord_name, + length(first(fixed_coords).values), + fixed_coords, + ) + ) + end + + # replace coordinates + varray_newcoords = FieldArray(varray.name, varray.values, Tuple(named_dimensions), varray.attributes) + + return varray_newcoords end \ No newline at end of file diff --git a/src/FieldRecord.jl b/src/FieldRecord.jl index b61c029..4096d41 100644 --- a/src/FieldRecord.jl +++ b/src/FieldRecord.jl @@ -129,8 +129,8 @@ function Base.copy(fr::FieldRecord{D, S, V, N, M, R}) where {D, S, V, N, M, R} end """ - get_array(fr::FieldRecord [, allselectargs::NamedTuple] [,coords_records::AbstractVector]) -> FieldArray - get_array(fr::FieldRecord; [allselectargs...]) -> FieldArray + get_array(fr::FieldRecord [, allselectargs::NamedTuple] [; coords::AbstractVector]) -> FieldArray + [deprecated] get_array(fr::FieldRecord [; coords::AbstractVector] [; allselectargs...]) -> FieldArray Return a [`FieldArray`](@ref) containing `fr::FieldRecord` data values and any attached coordinates, for records and spatial region defined by `allselectargs`. @@ -147,34 +147,62 @@ any attached coordinates, for records and spatial region defined by `allselectar Available `selectargs` depend on the grid `fr.mesh`, and are passed to `PB.Grids.get_region`. - -Optional argument `coords_records` can be used to supply coordinates from additional `FieldRecords`, replacing any coordinates +Optional argument `coords` can be used to supply coordinates from additional `FieldRecords`, replacing any coordinates attached to `fr`. Format is a Vector of Pairs of "coord_name"=>(cr1::FieldRecord, cr2::FieldRecord, ...) Example: to replace a 1D column default pressure coordinate with a z coordinate: - coords_records=["z"=>(zmid::FieldRecord, zlower::FieldRecord, zupper::FieldRecord)] + coords=["z"=>(zmid::FieldRecord, zlower::FieldRecord, zupper::FieldRecord)] NB: the coordinates will be generated by applying `allselectargs`, so the supplied coordinate FieldRecords must have the same dimensionality as `fr`. """ function get_array( - fr::FieldRecord; allselectargs... + fr::FieldRecord; + coords=nothing, + allselectargs... ) - Base.depwarn( - "allselectargs... will be deprecated in a future release. Please use allselectargs::NamedTuple instead", - :get_array, - ) - return get_array(fr, NamedTuple(allselectargs)) + isempty(allselectargs) || + Base.depwarn( + "allselectargs... will be deprecated in a future release. Please use allselectargs::NamedTuple instead", + :get_array, + ) + + return get_array(fr, NamedTuple(allselectargs); coords=coords) end function get_array( - fr::FieldRecord{D, S, V, N, M, R}, allselectargs::NamedTuple=NamedTuple(), + fr::FieldRecord, allselectargs::NamedTuple; # allselectargs::NamedTuple=NamedTuple() creates a method ambiguity with deprecated form above + coords::Union{Nothing, AbstractVector}=nothing +) + + # get data NB: with original coords or no coords + varray = _get_array(fr, allselectargs) + + if !isnothing(coords) + check_coords_argument(coords) || + error("coords argument should be of form 'coords=[\"z\"=>zmid::FieldRecord, zlower::FieldRecord, zupper::FieldRecord), ...]") + + # get arrays for replacement coordinates + vec_coord_arrays = [ + coord_name => Tuple(_get_array(coord_fr, allselectargs) for coord_fr in coord_fieldrecords) + for (coord_name, coord_fieldrecords) in coords + ] + + # replace coordinates + varray = update_coordinates(varray, vec_coord_arrays) + end + + return varray +end + +function _get_array( + fr::FieldRecord{D, S, V, N, M, R}, allselectargs::NamedTuple, ) where {D, S, V, N, M, R} # select records to use and create PB.NamedDimension ridx = 1:length(fr) - selectargs_region = Dict() + selectargs_region = NamedTuple() selectargs_records = NamedTuple() for (k, v) in zip(keys(allselectargs), allselectargs) if k==:records @@ -193,23 +221,18 @@ function get_array( end end else - selectargs_region[k] = v + selectargs_region = merge(selectargs_region, NamedTuple((k=>v,))) end end records_dim = PB.NamedDimension("records", length(ridx), PB.get_region(fr.coords_record, ridx)) # add attributes for selection used attributes = copy(fr.attributes) - attributes[:filter_records] = selectargs_records - attributes[:filter_region] = NamedTuple(selectargs_region) + isempty(selectargs_records) || (attributes[:filter_records] = selectargs_records;) + isempty(selectargs_region) || (attributes[:filter_region] = selectargs_region;) # Generate name from attributes - name = get(attributes, :domain_name, "") - name *= isempty(name) ? "" : "." - name *= get(attributes, :var_name, "") - if !isempty(selectargs_region) || !isempty(selectargs_records) - name *= "(" * join(["$k=$v" for (k, v) in merge(Dict(pairs(selectargs_records)), selectargs_region)], ", ") * ")" - end + name = default_fieldarray_name(attributes) # Select selectargs_region if field_single_element(fr) @@ -236,7 +259,7 @@ function get_array( # pass through to Field # get FieldArray from first Field record and use this to figure out array shapes etc - far = get_array(fr[first(ridx)]; selectargs_region...) + far = get_array(fr[first(ridx)], selectargs_region) if length(ridx) > 1 # add additional (last) dimension for records favalues = Array{eltype(far.values), length(far.dims)+1}(undef, size(far.values)..., length(ridx)) @@ -265,7 +288,7 @@ function get_array( # fill with values from FieldArrays for Fields for remaining records if length(ridx) > 1 for (i, r) in enumerate(ridx[2:end]) - far = get_array(fr[r]; selectargs_region...) + far = get_array(fr[r], selectargs_region) if isempty(far.dims) fa.values[i+1] = far.values else @@ -279,31 +302,3 @@ function get_array( end -# substitute or add coords_records to FieldArray(fr; selectargs...) -function get_array(fr::FieldRecord, allselectargs::NamedTuple, coords_records::AbstractVector) - - # get data NB: with original coords or no coords - varray = get_array(fr; allselectargs...) - - # generate Vector of NamedDimensions to use as new coordinates - named_dimensions = PB.NamedDimension[] - for (coord_name, coord_fieldrecords) in coords_records - fixed_coords = [] - for coord_fr in coord_fieldrecords - coord_array = get_array(coord_fr; allselectargs...) - push!(fixed_coords, PB.FixedCoord(get(coord_array.attributes, :var_name, ""), coord_array.values, coord_array.attributes)) - end - push!( - named_dimensions, PB.NamedDimension( - coord_name, - length(first(fixed_coords).values), - fixed_coords, - ) - ) - end - - # replace coordinates - varray_coords = FieldArray(varray.name, varray.values, Tuple(named_dimensions), varray.attributes) - - return varray_coords -end \ No newline at end of file diff --git a/src/OutputWriters.jl b/src/OutputWriters.jl index 7c57d2f..dbc4755 100644 --- a/src/OutputWriters.jl +++ b/src/OutputWriters.jl @@ -130,33 +130,36 @@ so the supplied coordinate Variables must have the same dimensionality as `vars` """ function PALEOmodel.get_array( output::PALEOmodel.AbstractOutputWriter, varname::AbstractString; + coords=nothing, allselectargs... ) - Base.depwarn( - "allselectargs... will be deprecated in a future release. Please use allselectargs::NamedTuple instead", - PALEOmodel.get_array - ) - return PALEOmodel.get_array(output, varname, NamedTuple(allselectargs)) + isempty(allselectargs) || + Base.depwarn( + "allselectargs... will be deprecated in a future release. Please use allselectargs::NamedTuple instead", + :get_array, + ) + return PALEOmodel.get_array(output, varname, NamedTuple(allselectargs); coords=coords) end function PALEOmodel.get_array( - output::PALEOmodel.AbstractOutputWriter, varname::AbstractString, allselectargs::NamedTuple=NamedTuple(); - coords::AbstractVector=[], + output::PALEOmodel.AbstractOutputWriter, varname::AbstractString, allselectargs::NamedTuple; # allselectargs::NamedTuple=NamedTuple() creates a method ambiguity with deprecated form above + coords=nothing, ) - fr = PB.get_field(output, varname) - - coords_records = [ - coord_name => (PB.get_field(output, cvn) for cvn in coord_varnames) - for (coord_name, coord_varnames) in coords - ] + fr = PB.get_field(output, varname) - if isempty(coords_records) - fa = PALEOmodel.get_array(fr, allselectargs) + if isnothing(coords) + coords_records=nothing else - fa = PALEOmodel.get_array(fr, allselectargs, coords_records) - end + PALEOmodel.check_coords_argument(coords) || + error("argument coords should be a Vector of Pairs of \"coord_name\"=>(\"var_name1\", \"var_name2\", ...), eg: [\"z\"=>(\"atm.zmid\", \"atm.zlower\", \"atm.zupper\"), ...]") + + coords_records = [ + coord_name => Tuple(PB.get_field(output, cvn) for cvn in coord_varnames) + for (coord_name, coord_varnames) in coords + ] + end - return fa + return PALEOmodel.get_array(fr, allselectargs; coords=coords_records) end """ diff --git a/src/PlotRecipes.jl b/src/PlotRecipes.jl index 4621e20..c1e63ca 100644 --- a/src/PlotRecipes.jl +++ b/src/PlotRecipes.jl @@ -116,9 +116,11 @@ end """ plot(output::AbstractOutputWriter, vars::Union{AbstractString, Vector{<:AbstractString}} [, selectargs::NamedTuple] [; coords::AbstractVector]) - heatmap(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)`](@ref)) @@ -141,30 +143,36 @@ so the supplied coordinate Variables must have the same dimensionality as `vars` """ RecipesBase.@recipe function f( output::AbstractOutputWriter, - vars::Union{AbstractString, Vector{<:AbstractString}}, + vars::Union{AbstractString, AbstractVector{<:AbstractString}}, selectargs::NamedTuple=NamedTuple(); - coords::AbstractVector = [], + coords=nothing, # NB: PlotRecipes doesn't support type annotation ) if isa(vars, AbstractString) vars = [vars] end - coords_records = [ - coord_name => (PB.get_field(output, cvn) for cvn in coord_varnames) - for (coord_name, coord_varnames) in coords - ] + if isnothing(coords) + coords_records=nothing + else + check_coords_argument(coords) || + error("coords argument should be of form 'coords=[\"z\"=>(\"atm.zmid\", \"atm.zlower\", \"atm.zupper\"), ...]") + coords_records = [ + coord_name => Tuple(PB.get_field(output, cvn) for cvn in coord_varnames) + for (coord_name, coord_varnames) in coords + ] + end delete!(plotattributes, :coords) for var in vars RecipesBase.@series begin - varsplit = split(var, ".") - if length(varsplit) == 3 - structfield := Symbol(varsplit[3]) - varsplit = varsplit[1:2] + # if var is of form .., set :structfield to take a single field from struct-valued Var + var, var_structfield = _split_var_structfield(var) + if !isnothing(var_structfield) + structfield := var_structfield end - var = join(varsplit, ".") + # pass through to plot recipe for FieldRecord PB.get_field(output, var), selectargs, coords_records end end @@ -172,6 +180,51 @@ RecipesBase.@recipe function f( return nothing end +RecipesBase.@recipe function f( + modeldata::PB.AbstractModelData, + vars::Union{AbstractString, AbstractVector{<:AbstractString}}, + selectargs::NamedTuple=NamedTuple(); + coords=nothing, # NB: PlotRecipes doesn't support type annotation +) + delete!(plotattributes, :coords) + + if isa(vars, AbstractString) + vars = [vars] + end + + # broadcast any Vector-valued argument in selectargs + bcastargs = broadcast_dict([Dict{Symbol, Any}(pairs(selectargs))]) + + for var in vars + # if var is of form .., set :structfield to take a single field from struct-valued Var + var, var_structfield = _split_var_structfield(var) + if !isnothing(var_structfield) + structfield := var_structfield + end + + # broadcast selectargs + for sa in bcastargs + sa_nt = NamedTuple(sa) + RecipesBase.@series begin + # pass through to plot recipe for FieldArray + PALEOmodel.get_array(modeldata, var, sa_nt; coords=coords) + end + end + end + + return nothing +end + +function _split_var_structfield(var) + # if var is of form .., split off + varsplit = split(var, ".") + if length(varsplit) == 3 + return join(varsplit[1:2], "."), Symbol(varsplit[3]) + else + return var, nothing + end +end + """ plot(outputs::Vector{<:AbstractOutputWriter}, vars::Union{AbstractString, Vector{<:AbstractString}}, selectargs::NamedTuple=NamedTuple()) @@ -182,17 +235,15 @@ adding a `labelprefix` (index in `outputs` Vector) to identify each plot series RecipesBase.@recipe function f( outputs::Vector{<:AbstractOutputWriter}, vars::Union{AbstractString, Vector{<:AbstractString}}, - selectargs::NamedTuple=NamedTuple(); - coords::AbstractVector=[], + selectargs::NamedTuple=NamedTuple() ) for (i, output) in enumerate(outputs) RecipesBase.@series begin labelprefix --> "$i: " - output, vars, selectargs, coords + output, vars, selectargs end end - delete!(plotattributes, :coords) - + return nothing end @@ -211,11 +262,11 @@ Vector-valued fields in `selectargs` are "broadcasted" (generating a separate pl Optional argument `coords_records` can be used to supply plot coordinates from `FieldRecords`. Format is a Vector of Pairs of "coord_name"=>(cr1::FieldRecord, cr2::FieldRecord, ...) Example: - coords=["z"=>(zmid::FieldRecord, zlower::FieldRecord, zupper::FieldRecord)] + coords_records=["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`. """ -RecipesBase.@recipe function f(fr::FieldRecord, selectargs::NamedTuple, coords_records::AbstractVector=[]) +RecipesBase.@recipe function f(fr::FieldRecord, selectargs::NamedTuple, coords_records=nothing) # broadcast any Vector-valued argument in selectargs bcastargs = broadcast_dict([Dict{Symbol, Any}(pairs(selectargs))]) @@ -223,11 +274,7 @@ RecipesBase.@recipe function f(fr::FieldRecord, selectargs::NamedTuple, coords_r for sa in bcastargs sa_nt = NamedTuple(sa) RecipesBase.@series begin - if isempty(coords_records) - get_array(fr, sa_nt) - else - get_array(fr, sa_nt, coords_records) - end + get_array(fr, sa_nt; coords=coords_records) end end end From 2bc7397eaea0cda965d909578ede18b223dfc2b8 Mon Sep 17 00:00:00 2001 From: Stuart Daines Date: Sun, 31 Jul 2022 20:25:00 +0100 Subject: [PATCH 4/4] Workaround for apparent Julia compiler bug (?) triggered by get_array changes See L262 in FieldRecord.jl : --- src/FieldArray.jl | 8 ++++---- src/FieldRecord.jl | 3 +++ 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/FieldArray.jl b/src/FieldArray.jl index 7d42968..a15c676 100644 --- a/src/FieldArray.jl +++ b/src/FieldArray.jl @@ -89,8 +89,8 @@ function get_array( values, dims = PB.Grids.get_region(f.mesh, f.values; selectargs...) - if !isempty(selectargs) - attributes = isnothing(attributes) ? Dict{Symbol, Any}() : copy(attributes) + if !isnothing(attributes) && !isempty(selectargs) + attributes = copy(attributes) attributes[:filter_region] = selectargs end @@ -125,8 +125,8 @@ function get_array( end end - if !isempty(selectargs) - attributes = isnothing(attributes) ? Dict{Symbol, Any}() : copy(attributes) + if !isnothing(attributes) && !isempty(selectargs) + attributes = copy(attributes) attributes[:filter_region] = selectargs end diff --git a/src/FieldRecord.jl b/src/FieldRecord.jl index 4096d41..d7c3ea1 100644 --- a/src/FieldRecord.jl +++ b/src/FieldRecord.jl @@ -260,6 +260,9 @@ function _get_array( # get FieldArray from first Field record and use this to figure out array shapes etc far = get_array(fr[first(ridx)], selectargs_region) + # TODO - Julia bug ? length(far.dims) returning wrong value, apparently triggered by this line + # attributes = isnothing(attributes) ? Dict{Symbol, Any}() : copy(attributes) + # in get_array if length(ridx) > 1 # add additional (last) dimension for records favalues = Array{eltype(far.values), length(far.dims)+1}(undef, size(far.values)..., length(ridx))