From 1b9daa5b9021fb11fe6b59488ed7ef3ba5a0ba34 Mon Sep 17 00:00:00 2001 From: Stuart Daines Date: Sun, 29 Dec 2024 22:42:36 +0000 Subject: [PATCH] Tidy up and remove old code - remove old unused code - fix doc build --- docs/src/PALEOmodelOutput.md | 7 +- src/CoordsDims.jl | 175 +------------------------------ src/FieldArray.jl | 187 --------------------------------- src/FieldRecord.jl | 181 ++++---------------------------- src/PALEOmodel.jl | 4 +- src/PlotRecipes.jl | 144 ++++++++++++++++---------- src/Regions.jl | 194 ----------------------------------- 7 files changed, 116 insertions(+), 776 deletions(-) delete mode 100644 src/Regions.jl diff --git a/docs/src/PALEOmodelOutput.md b/docs/src/PALEOmodelOutput.md index 3ca4c53..275a9f1 100644 --- a/docs/src/PALEOmodelOutput.md +++ b/docs/src/PALEOmodelOutput.md @@ -69,8 +69,6 @@ OutputMemoryDomain ```@docs save_netcdf load_netcdf! -save_jld2 -load_jld2! ``` ## Field Array @@ -78,11 +76,12 @@ load_jld2! ```@meta CurrentModule = PALEOmodel ``` -[`FieldArray`](@ref) provides a generic array type with named dimensions `PALEOboxes.NamedDimension` and optional coordinates `PALEOboxes.FixedCoord` for processing of model output. +[`FieldArray`](@ref) provides a generic array type with named dimensions `PALEOboxes.NamedDimension` and optional coordinates [`PALEOmodel.FixedCoord`](@ref) for processing of model output. ```@docs FieldArray -get_array(f::PB.Field{D, PB.ScalarSpace, V, N, M}; attributes=nothing) where {D, V, N, M} +get_array(fr::FieldRecord; coords=nothing) +FixedCoord ``` ## Plotting output diff --git a/src/CoordsDims.jl b/src/CoordsDims.jl index 79a9e67..f7d7def 100644 --- a/src/CoordsDims.jl +++ b/src/CoordsDims.jl @@ -8,6 +8,8 @@ A fixed (state independent) coordinate +These are generated from coordinate variables for use in output visualisation. + N = 1: a cell-centre coordinate, size(values) = (ncells, ) N = 2: a boundary coordinate, size(values) = (2, ncells) """ @@ -19,90 +21,11 @@ end is_boundary_coordinate(fc::FixedCoord) = ndims(fc.values) == 2 && size(fc.values, 1) == 2 -""" - append_units(name::AbstractString, attributes) -> "name (units)" - -Utility function to append variable units string to a variable name for display. -""" -function append_units(name::AbstractString, attributes::Dict{Symbol, Any}) - units = get(attributes, :units, "") - if isempty(units) - return name - else - return name*" ($units)" - end -end - -append_units(name::AbstractString, attributes::Nothing) = name - -""" - build_coords_edges(coords_vec::Vector{FixedCoord}) -> Vector{Float64} - -Build a vector of coordinate edges (length `n+1``) from `coords_vec`, assuming the PALEO -convention that `coords_vec` contains three elements with -cell midpoints, lower edges, upper edges each of length `n`, in that order. - -Falls back to just returning the first entry in `coords_vec` for other cases. -""" -function build_coords_edges(dim::PB.NamedDimension, coords_vec::Vector{FixedCoord}) - - if isempty(coords_vec) - # no coordinate - just return indices - co_values = collect(1:(dim.size+1)) - co_label = dim.name - elseif length(coords_vec) == 1 && is_boundary_coordinate(coords_vec[1]) - co = coords_vec[1] - co_values = vcat(co.values[1, :], co.values[2, end]) - co_label = append_units(co.name, co.attributes) - elseif length(coords_vec) == 1 || length(coords_vec) > 3 - # 1 coordinate or something we don't understand - take first - co = first(coords_vec) - co_values = co.values - co_label = append_units(co.name, co.attributes) - elseif length(coords_vec) == 2 && is_boundary_coordinate(coords_vec[2]) - # assume mid, boundary - co = coords_vec[2] - co_values = vcat(co.values[1, :], co.values[2, end]) - co_label = append_units(co.name, co.attributes) - elseif length(coords_vec) in (2, 3) - # 2 coordinates assume lower, upper edges - # 3 coordinates assume mid, lower, upper - co_lower = coords_vec[end-1] - co_upper = coords_vec[end] - co_label = append_units(co_lower.name*", "*co_upper.name, co_lower.attributes) - first(co_lower.values) < first(co_upper.values) || - @warn "build_coords_edges: $co_label co_lower is > co_upper - check model grid" - if co_lower.values[end] > co_lower.values[1] # ascending order - co_lower.values[2:end] == co_upper.values[1:end-1] || - @warn "build_coords_edges: $co_label lower and upper edges don't match" - co_values = [co_lower.values; co_upper.values[end]] - else # descending order - co_lower.values[1:end-1] == co_upper.values[2:end] || - @warn "build_coords_edges: $co_label lower and upper edges don't match" - co_values = [co_upper.values[1]; co_lower.values] - end - - end - return co_values, co_label -end - -"guess coordinate edges from midpoints, assuming uniform spacing" -function guess_coords_edges(x_midpoints) - first_x = x_midpoints[1] - 0.5*(x_midpoints[2] - x_midpoints[1]) - last_x = x_midpoints[end] + 0.5*(x_midpoints[end] - x_midpoints[end-1]) - return [first_x; 0.5.*(x_midpoints[1:end-1] .+ x_midpoints[2:end]); last_x] -end - - -function get_region(fc::FixedCoord, indices::AbstractVector) - return FixedCoord(fc.name, fc.values[indices], fc.attributes) -end - -function get_region(fcv::Vector{FixedCoord}, indices::AbstractVector) - return [FixedCoord(fc.name, fc.values[indices], fc.attributes) for fc in fcv] -end +################################################################## +# Coordinate filtering and selection +################################################################ "find indices of coord from first before range[1] to first after range[2]" function find_indices(coord::AbstractVector, range) @@ -180,91 +103,3 @@ function dimscoord_subset(dim::PB.NamedDimension, coords::Vector{FixedCoord}, se return cidx, dim_subset, coords_subset, coords_used end - -#= -################################################# -# Dimensions -##################################################### - -""" - NamedDimension - -A named dimension, with optional attached fixed coordinates `coords` - -PALEO convention is that where possible `coords` contains three elements, for cell -midpoints, lower edges, upper edges, in that order. -""" -mutable struct NamedDimension - name::String - size::Int64 - coords::Vector{FixedCoord} # may be empty -end - -"create from size only (no coords)" -function NamedDimension(name, size::Integer) - return NamedDimension( - name, - size, - FixedCoord[], - ) -end - -"create from coord mid-points" -function NamedDimension(name, coord::AbstractVector) - return NamedDimension( - name, - length(coord), - [ - FixedCoord(name, coord, Dict{Symbol, Any}()), - ] - ) -end - -"create from coord mid-points and edges" -function NamedDimension(name, coord::AbstractVector, coord_edges::AbstractVector) - if coord[end] > coord[1] - # ascending order - coord_lower = coord_edges[1:end-1] - coord_upper = coord_edges[2:end] - else - # descending order - coord_lower = coord_edges[2:end] - coord_upper = coord_edges[1:end-1] - end - return NamedDimension( - name, - length(coord), - [ - FixedCoord(name, coord, Dict{Symbol, Any}()), - FixedCoord(name*"_lower", coord_lower, Dict{Symbol, Any}()), - FixedCoord(name*"_upper", coord_upper, Dict{Symbol, Any}()), - ] - ) -end - -function get_region(nd::NamedDimension, indices::AbstractVector) - return NamedDimension(nd.name, length(indices), get_region(nd.coords, indices)) -end - -""" - build_coords_edges(nd::NamedDimension) -> Vector{Float64} - -Call [`build_coords_edges`](@ref)(nd.coords), or fallback to just returning indices -if no coords present. -""" -function build_coords_edges(nd::NamedDimension) - if !isempty(nd.coords) - return build_coords_edges(nd.coords) - else - @warn "no coords for NamedDimension $(nd.name), returning indices" - return collect(1:nd.size), nd.name*" (indices)" - end -end - -function Base.show(io::IO, nd::NamedDimension) - print(io, "NamedDimension(name=", nd.name, ", size=", nd.size, ", coords=(") - join(io, [c.name for c in nd.coords], ", ") - print(io, "))") - return nothing -end - =# \ No newline at end of file diff --git a/src/FieldArray.jl b/src/FieldArray.jl index 6bdbea8..7ff952c 100644 --- a/src/FieldArray.jl +++ b/src/FieldArray.jl @@ -1,4 +1,3 @@ -import Infiltrator """ FieldArray @@ -76,189 +75,3 @@ Get FieldArray from PALEO object `obj` """ function get_array end -""" - get_array(dataset, recordidx, 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`. - -Available `selectargs` depend on the grid `f.mesh`, and -are passed to `get_region`. - -`attributes` (if present) are added to `FieldArray` -""" -function get_array( - f::PB.Field{FieldData, PB.ScalarSpace, V, N, Mesh}, selectargs::NamedTuple=NamedTuple(); - attributes=nothing, coord_source, -) where {FieldData, V, N, Mesh} - isempty(selectargs) || - error("get_array on Field f defined on ScalarSpace with non-empty selectargs=$selectargs") - - data_dims_coordinates = get_data_dims_coords(coord_source, f.data_dims) - return FieldArray(default_fieldarray_name(attributes), f.values, data_dims_coordinates, attributes) -end - -function get_array( - f::PB.Field{FieldData, PB.CellSpace, V, 0, Mesh}, selectargs::NamedTuple=NamedTuple(); - attributes=nothing, coord_source=nothing, -) where {FieldData, V, Mesh} - - values, dims_coords = get_region(f.mesh, f.values; coord_source, selectargs...) - - if !isnothing(attributes) && !isempty(selectargs) - attributes = copy(attributes) - attributes[:filter_region] = selectargs - end - - return FieldArray(default_fieldarray_name(attributes), values, dims_coords, attributes) -end - -# single data dimension -# TODO generalize this to arbitrary data dimensions -function get_array( - f::PB.Field{FieldData, PB.CellSpace, V, 1, Mesh}, selectargs::NamedTuple=NamedTuple(); - attributes=nothing, coord_source, -) where {FieldData, V, Mesh} - - f_space_dims_colons = ntuple(i->Colon(), ndims(f.values) - 1) - f_size_datadim = size(f.values)[end] - - dvalues, dims_coords = get_region(f.mesh, f.values[f_space_dims_colons..., 1]; coord_source, selectargs...) - - d = (size(dvalues)..., f_size_datadim) - values = Array{eltype(dvalues), length(d)}(undef, d...) - - if length(d) == 1 - # single cell - space dimension squeezed out - for i in 1:f_size_datadim - values[i], dims_coords = get_region(f.mesh, f.values[f_space_dims_colons..., i]; coord_source, selectargs...) - end - else - dvalues_colons = ntuple(i->Colon(), ndims(dvalues)) - for i in 1:f_size_datadim - dvalues, dims_coords = get_region(mesh, f.values[f_space_dims_colons..., i]; coord_source, selectargs...) - values[dvalues_colons..., i] .= dvalues - end - end - - if !isnothing(attributes) && !isempty(selectargs) - attributes = copy(attributes) - attributes[:filter_region] = selectargs - end - - data_dims_coords = get_data_dims_coords(coord_source, f.data_dims) - - return FieldArray(default_fieldarray_name(attributes), values, vcat(dims_coords, data_dims_coords), attributes) -end - - - -""" - 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 `get_region`). - -Optional argument `coords` can be used to supply plot coordinates from Variable in output. -Format is a Vector of Pairs of "dim_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( - modeldata::PB.AbstractModelData, varnamefull::AbstractString, selectargs::NamedTuple=NamedTuple(); - coords=nothing, -) - varsplit = split(varnamefull, ".") - length(varsplit) == 2 || - throw(ArgumentError("varnamefull $varnamefull is not of form .")) - domainname, varname = varsplit - - coord_source = (modeldata, domainname) - f, attributes = _get_field_attributes(coord_source, varname) - - varray = get_array(f, selectargs; attributes, coord_source) - - if isnothing(coords) - # keep original coords (if any) - else - check_coords_argument(coords) || - error("argument coords should be a Vector of Pairs of \"dim_name\"=>(\"var_name1\", \"var_name2\", ...), eg: [\"z\"=>(\"atm.zmid\", \"atm.zlower\", \"atm.zupper\"), ...]") - - vec_coords_arrays = [ - dim_name => Tuple(get_array(modeldata, cvn, selectargs) for cvn in coord_varnames) - for (dim_name, coord_varnames) in coords - ] - - varray = update_coordinates(varray, vec_coords_arrays) - end - - return varray -end - -function get_data_dims_coords(coord_source, data_dims) - - data_dims_coordinates = Pair{PB.NamedDimension, Vector{PALEOmodel.FixedCoord}}[] - - for dd in data_dims - coordinates = PALEOmodel.FixedCoord[] - for coord_name in PB.get_coordinates(dd.name) - ddf, attributes = _get_field_attributes(coord_source, coord_name) - push!(coordinates, FixedCoord(coord_name, ddf.values, attributes)) - end - push!(data_dims_coordinates, dd => coordinates) - end - - return data_dims_coordinates -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 "dim_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 \"dim_name\"=>(var1::FieldArray, var2::FieldArray, ...), eg: [\"z\"=>(zmid::FieldArray, zlower::FieldArray, atm.zupper::FieldArray)]") - - # generate Vector of NamedDimensions to use as new coordinates - new_dims_coords = Pair{PB.NamedDimension, Vector{FixedCoord}}[nd => copy(coords) for (nd, coords) in varray.dims_coords] - for (dim_name, coord_arrays) in vec_coord_arrays - dc_idx = findfirst(dc->first(dc).name == dim_name, new_dims_coords) - !isnothing(dc_idx) || error("FieldArray $varray has no dimension $dim_name") - nd, coords = new_dims_coords[dc_idx] - empty!(coords) - for coord_array in coord_arrays - coord_array_name = get(coord_array.attributes, :var_name, "") - nd.size == length(coord_array.values) || - error("FieldArray $varray dimension $dim_name size $(nd.size) != new coordinate $coord_array_name length(values) $(coord_array.values)") - push!(coords, PB.FixedCoord(coord_array_name, coord_array.values, coord_array.attributes)) - end - end - - # replace coordinates - varray_newcoords = FieldArray(varray.name, varray.values, new_dims_coords, varray.attributes) - - return varray_newcoords -end \ No newline at end of file diff --git a/src/FieldRecord.jl b/src/FieldRecord.jl index f039238..843c1af 100644 --- a/src/FieldRecord.jl +++ b/src/FieldRecord.jl @@ -14,7 +14,6 @@ struct FieldRecord{FieldData <: PB.AbstractData, Space <: PB.AbstractSpace, V, N data_dims::NTuple{N, PB.NamedDimension} mesh::Mesh attributes::Dict{Symbol, Any} - # coords_record::Vector{FixedCoord} # coordinates attached to record dimension end # create empty FieldRecord @@ -201,21 +200,15 @@ any attached coordinates, for records and spatial region defined by `allselectar starting at the record with the nearest value of `fr.coords_record` before `first` and ending at the nearest record after `last`. -Available `selectargs` depend on the grid `fr.mesh` and `fr.data_dims`, options include: -- grid spatial dimensions -[- TODO grid spatial coordinates (if defined)] -- `fr.data_dims` +Available `selectargs` depend on the FieldRecord 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`. +Optional argument `coords` can be used to replace the attached coordinates for one or more dimensions. +Format is a Vector of Pairs of `""=>("", "", ...)`, +eg to replace a 1D column default pressure coordinate with a z coordinate: -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 "dim_name"=>(cr1::FieldRecord, cr2::FieldRecord, ...) + coords=["cells"=>("atm.zmid", "atm.zlower", "atm.zupper")] -Example: to replace a 1D column default pressure coordinate with a z coordinate: - - 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; @@ -231,163 +224,17 @@ function get_array( return get_array(fr, NamedTuple(allselectargs); coords=coords) end -function get_array_orig( - fr::FieldRecord, allselectargs::NamedTuple; # allselectargs::NamedTuple=NamedTuple() creates a method ambiguity with deprecated form above - coords::Union{Nothing, AbstractVector}=nothing -) - - varray = nothing - - try - # 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 - catch - frname = get(fr.attributes, :domain_name, "")*"."*get(fr.attributes, :var_name, "") - @warn "PALEOmodel.get_array(fr::FieldRecord) failed for variable $frname" - rethrow() - end - - return varray -end - -function _get_array( - fr::FieldRecord{FieldData, Space, V, N, Mesh, R}, allselectargs::NamedTuple, -) where {FieldData, Space, V, N, Mesh, R} - - # select records to use and create PB.NamedDimension - ridx = 1:length(fr) - selectargs_region = NamedTuple() - selectargs_records = NamedTuple() - for (k, v) in zip(keys(allselectargs), values(allselectargs)) - if k==:records - if v isa Integer - ridx = [v] - else - ridx = v - end - selectargs_records = (records=v,) - elseif String(k) in fr.dataset.record_dim_coordinates - # find ridx corresponding to a coordinate - rfr = PB.get_field(fr.dataset, String(k)) - ridx, cvalue = find_indices(rfr.records, v) - selectargs_records = NamedTuple((k=>cvalue,)) - else - selectargs_region = merge(selectargs_region, NamedTuple((k=>v,))) - end - end - records_coords = FixedCoord[] - for rcn in fr.dataset.record_dim_coordinates - rfr = PB.get_field(fr.dataset, rcn) - push!(records_coords, FixedCoord(rcn, rfr.records[ridx], rfr.attributes)) - end - records_dim_coords = PB.NamedDimension(fr.dataset.record_dim.name, length(ridx)) => records_coords - - # add attributes for selection used - attributes = copy(fr.attributes) - isempty(selectargs_records) || (attributes[:filter_records] = selectargs_records;) - isempty(selectargs_region) || (attributes[:filter_region] = selectargs_region;) - - # Generate name from attributes - name = default_fieldarray_name(attributes) - - # Select selectargs_region - if field_single_element(fr) - # create FieldArray directly from FieldRecord - data_dims_coords = get_data_dims_coords((fr.dataset, ridx), fr.data_dims) - - isempty(selectargs_region) || - throw(ArgumentError("invalid index $selectargs_region")) - if length(ridx) > 1 - return FieldArray( - name, - fr.records[ridx], - vcat(data_dims_coords, records_dim_coords), - attributes - ) - else - # squeeze out records dimension - return FieldArray( - name, - fr.records[first(ridx)], - data_dims_coords, - attributes - ) - end - else - # 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; coord_source=(fr.dataset, first(ridx))) - # 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_coords)+1}(undef, size(far.values)..., length(ridx)) - fa = FieldArray( - name, - favalues, - vcat(far.dims_coords, records_dim_coords), - attributes, - ) - # copy values for first record - if isempty(far.dims_coords) - fa.values[1] = far.values - else - fa.values[fill(Colon(), length(far.dims_coords))..., 1] .= far.values - end - else - # squeeze out record dimension so this is just fa with attributes added - fa = FieldArray( - name, - far.values, - far.dims_coords, - attributes - ) - end - - # 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; coord_source=(fr.dataset, r)) - if isempty(far.dims_coords) - fa.values[i+1] = far.values - else - fa.values[fill(Colon(), length(far.dims_coords))..., i+1] .= far.values - end - end - end - - return fa - end - -end - function get_array( fr::FieldRecord, allselectargs::NamedTuple; # allselectargs::NamedTuple=NamedTuple() creates a method ambiguity with deprecated form above coords=nothing, expand_cartesian=false, # can be overridden in allselectargs squeeze_all_single_dims=true, # can be overridden in allselectargs - verbose=true, + verbose=false, ) frname = default_fieldarray_name(fr.attributes) - @info "get_array (begin): $frname, allselectargs $allselectargs" + verbose && @info "get_array (begin): $frname, allselectargs $allselectargs" if !isnothing(coords) check_coords_argument(coords) || @@ -561,7 +408,7 @@ function get_array( # Generate name from attributes, including selection suffixes name = default_fieldarray_name(attributes) - @info "get_array (end): $frname -> $name, allselectargs $allselectargs" + verbose && @info "get_array (end): $frname -> $name, allselectargs $allselectargs" return FieldArray( name, @@ -571,6 +418,16 @@ function get_array( ) 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}) + ) + ) + # Read coordinates attached to dimension 'dim' # For record coordinate, ridx_to_use = nothing # For other coordinates, ridx_to_use = record index to use for coordinate (ie limitation only one record in use, or coordinate is constant) diff --git a/src/PALEOmodel.jl b/src/PALEOmodel.jl index 24be4a1..cfb5a91 100644 --- a/src/PALEOmodel.jl +++ b/src/PALEOmodel.jl @@ -62,9 +62,7 @@ include("SteadyState.jl") include("SteadyStateKinsol.jl") -include("CoordsDims.jl") # TODO - -include("Regions.jl") # TODO +include("CoordsDims.jl") include("FieldArray.jl") diff --git a/src/PlotRecipes.jl b/src/PlotRecipes.jl index 8cc012e..024f9e7 100644 --- a/src/PlotRecipes.jl +++ b/src/PlotRecipes.jl @@ -133,9 +133,6 @@ end 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)) @@ -152,8 +149,6 @@ 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, @@ -166,16 +161,6 @@ RecipesBase.@recipe function f( vars = [vars] end - # 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 @@ -194,41 +179,6 @@ 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, ".") @@ -273,12 +223,13 @@ Calls [`get_array(fr, selectargs)`](@ref) and passes on to `plot(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 "coord_name"=>(cr1::FieldRecord, cr2::FieldRecord, ...) -Example: - 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`. +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")] + """ RecipesBase.@recipe function f(fr::FieldRecord, selectargs::NamedTuple, coords=nothing) @@ -527,3 +478,84 @@ function broadcast_dict(dv::Vector{<:Dict}) return dvout end end + + +######################################################## +# Helper functions for output visualisation +######################################################### + +""" + build_coords_edges(coords_vec::Vector{FixedCoord}) -> Vector{Float64} + +Build a vector of coordinate edges (length `n+1``) from `coords_vec`, assuming the PALEO +convention that `coords_vec` contains three elements with +cell midpoints, lower edges, upper edges each of length `n`, in that order. + +Falls back to just returning the first entry in `coords_vec` for other cases. +""" +function build_coords_edges(dim::PB.NamedDimension, coords_vec::Vector{FixedCoord}) + + if isempty(coords_vec) + # no coordinate - just return indices + co_values = collect(1:(dim.size+1)) + co_label = dim.name + elseif length(coords_vec) == 1 && is_boundary_coordinate(coords_vec[1]) + co = coords_vec[1] + co_values = vcat(co.values[1, :], co.values[2, end]) + co_label = append_units(co.name, co.attributes) + elseif length(coords_vec) == 1 || length(coords_vec) > 3 + # 1 coordinate or something we don't understand - take first + co = first(coords_vec) + co_values = co.values + co_label = append_units(co.name, co.attributes) + elseif length(coords_vec) == 2 && is_boundary_coordinate(coords_vec[2]) + # assume mid, boundary + co = coords_vec[2] + co_values = vcat(co.values[1, :], co.values[2, end]) + co_label = append_units(co.name, co.attributes) + elseif length(coords_vec) in (2, 3) + # 2 coordinates assume lower, upper edges + # 3 coordinates assume mid, lower, upper + co_lower = coords_vec[end-1] + co_upper = coords_vec[end] + co_label = append_units(co_lower.name*", "*co_upper.name, co_lower.attributes) + first(co_lower.values) < first(co_upper.values) || + @warn "build_coords_edges: $co_label co_lower is > co_upper - check model grid" + if co_lower.values[end] > co_lower.values[1] # ascending order + co_lower.values[2:end] == co_upper.values[1:end-1] || + @warn "build_coords_edges: $co_label lower and upper edges don't match" + co_values = [co_lower.values; co_upper.values[end]] + else # descending order + co_lower.values[1:end-1] == co_upper.values[2:end] || + @warn "build_coords_edges: $co_label lower and upper edges don't match" + co_values = [co_upper.values[1]; co_lower.values] + end + + end + + return co_values, co_label +end + +"guess coordinate edges from midpoints, assuming uniform spacing" +function guess_coords_edges(x_midpoints) + first_x = x_midpoints[1] - 0.5*(x_midpoints[2] - x_midpoints[1]) + last_x = x_midpoints[end] + 0.5*(x_midpoints[end] - x_midpoints[end-1]) + return [first_x; 0.5.*(x_midpoints[1:end-1] .+ x_midpoints[2:end]); last_x] +end + + +""" + append_units(name::AbstractString, attributes) -> "name (units)" + +Utility function to append variable units string to a variable name for display. +""" +function append_units(name::AbstractString, attributes::Dict{Symbol, Any}) + units = get(attributes, :units, "") + if isempty(units) + return name + else + return name*" ($units)" + end +end + +append_units(name::AbstractString, attributes::Nothing) = name diff --git a/src/Regions.jl b/src/Regions.jl deleted file mode 100644 index 924db7a..0000000 --- a/src/Regions.jl +++ /dev/null @@ -1,194 +0,0 @@ -""" - get_region(grid::Union{PB.AbstractMesh, Nothing}, values; coord_source, selectargs...) - -> values_subset, dim_subset::Vector{Pair{PB.NamedDimension, Vector{PALEOmodel.FixedCoord}}} - -Return the subset of `values` given by `selectargs` (Grid-specific keywords eg cell=, column=, ...) -and corresponding dimensions (with attached coordinates). -""" -function get_region(grid::Union{PB.AbstractMesh, Nothing}, values) end - -""" - get_region(dataset, recordidx, grid::Nothing, values) -> values[], [] - -Fallback for Domain with no grid, assumed 1 cell -""" -function get_region(grid::Nothing, values; coord_source=nothing) - length(values) == 1 || - throw(ArgumentError("grid==Nothing and length(values) != 1")) - return values[], Pair{PB.NamedDimension, Vector{PALEOmodel.FixedCoord}}[] -end - - -""" - get_region(dataset, recordidx, grid::PB.Grids.UnstructuredVectorGrid, values; cell=Nothing) -> - values_subset, dim_subset - -# Keywords for region selection: -- `cell::Union{Nothing, Int, Symbol}`: an Int for cell number (first cell = 1), or a Symbol to look up in `cellnames` - `cell = Nothing` is also allowed for a single-cell Domain. -""" -function get_region( - grid::PB.Grids.UnstructuredVectorGrid, values; - cell::Union{Nothing, Int, Symbol}=nothing, - coord_source=nothing, -) - if isnothing(cell) - grid.ncells == 1 || throw(ArgumentError("'cell' argument (an Int or Symbol) is required for an UnstructuredVectorGrid with > 1 cell")) - idx = 1 - elseif cell isa Int - idx = cell - else - idx = get(grid.cellnames, cell, nothing) - !isnothing(idx) || - throw(ArgumentError("cell ':$cell' not present in grid.cellnames=$(grid.cellnames)")) - end - - return ( - values[idx], - Pair{PB.NamedDimension, Vector{PALEOmodel.FixedCoord}}[], # no dimensions (ie squeeze out a dimension length 1 for single cell) - ) -end - - - -""" - get_region(dataset, recordidx, grid::PB.Grids.UnstructuredColumnGrid, values; column, [cell=nothing]) -> - values_subset, (dim_subset::NamedDimension, ...) - -# Keywords for region selection: -- `column::Union{Int, Symbol}`: (may be an Int, or a Symbol to look up in `columnames`) -- `cell::Int`: optional cell index within `column`, highest cell is cell 1 -""" -function get_region( - grid::PB.Grids.UnstructuredColumnGrid, values; - column, - cell::Union{Nothing, Int}=nothing, - coord_source=nothing -) - - if column isa Int - column in 1:length(grid.Icolumns) || - throw(ArgumentError("column index $column out of range")) - colidx = column - else - colidx = findfirst(isequal(column), grid.columnnames) - !isnothing(colidx) || - throw(ArgumentError("columnname '$column' not present in grid.columnnames=$(grid.columnnames)")) - end - - if isnothing(cell) - indices = grid.Icolumns[colidx] - coordnames = PB.get_coordinates(grid, "cells") - coordinates = FixedCoord[] - for cn in coordnames - cf, attributes = _get_field_attributes(coord_source, cn) - push!(coordinates, FixedCoord(cn, cf.values[indices], attributes)) - end - return ( - values[indices], - [PB.NamedDimension("cells", length(indices)) => coordinates], - ) - else - # squeeze out z dimension - idx = grid.Icolumns[colidx][cell] - return ( - values[idx], - Pair{PB.NamedDimension, Vector{PALEOmodel.FixedCoord}}[], # no dimensions (ie squeeze out a dimension length 1 for single cell) - ) - end - -end - - -""" - get_region(dataset, recordidx, rid::Union{PB.Grids.CartesianLinearGrid{2}, PB.Grids.CartesianArrayGrid{2}} , internalvalues; [i=i_idx], [j=j_idx]) -> - arrayvalues_subset, (dim_subset::NamedDimension, ...) - -# Keywords for region selection: -- `i::Int`: optional, slice along first dimension -- `j::Int`: optional, slice along second dimension - -`internalvalues` are transformed if needed from internal Field representation as a Vector length `ncells`, to -an Array (2D if neither i, j arguments present, 1D if i or j present, 0D ie one cell if both present) -""" -function get_region( - grid::Union{PB.Grids.CartesianLinearGrid{2}, PB.Grids.CartesianArrayGrid{2}}, internalvalues; - i::Union{Integer, Colon}=Colon(), - j::Union{Integer, Colon}=Colon(), - coord_source, -) - return _get_region(grid, internalvalues, [i, j], coord_source) -end - -""" - get_region(grid::Union{PB.Grids.CartesianLinearGrid{3}, PB.Grids.CartesianArrayGrid{3}}, internalvalues; [i=i_idx], [j=j_idx]) -> - arrayvalues_subset, (dim_subset::NamedDimension, ...) - -# Keywords for region selection: -- `i::Int`: optional, slice along first dimension -- `j::Int`: optional, slice along second dimension -- `k::Int`: optional, slice along third dimension - -`internalvalues` are transformed if needed from internal Field representation as a Vector length `ncells`, to -an Array (3D if neither i, j, k arguments present, 2D if one of i, j or k present, 1D if two present, -0D ie one cell if i, j, k all specified). -""" -function get_region( - grid::Union{PB.Grids.CartesianLinearGrid{3}, PB.Grids.CartesianArrayGrid{3}}, internalvalues; - i::Union{Integer, Colon}=Colon(), - j::Union{Integer, Colon}=Colon(), - k::Union{Integer, Colon}=Colon(), - coord_source, -) - return _get_region(grid, internalvalues, [i, j, k], coord_source) -end - -function _get_region( - grid::Union{PB.Grids.CartesianLinearGrid, PB.Grids.CartesianArrayGrid}, internalvalues, indices, coord_source -) - if !isempty(grid.coords) && !isempty(grid.coords_edges) - dims = [ - PB.NamedDimension(grid.dimnames[idx], grid.coords[idx], grid.coords_edges[idx]) - for (idx, ind) in enumerate(indices) if isa(ind, Colon) - ] - elseif !isempty(grid.coords) - dims = [ - PB.NamedDimension(grid.dimnames[idx], grid.coords[idx]) - for (idx, ind) in enumerate(indices) if isa(ind, Colon) - ] - else - dims = [ - PB.NamedDimension(grid.dimnames[idx]) - for (idx, ind) in enumerate(indices) if isa(ind, Colon) - ] - end - - values = internal_to_cartesian(grid, internalvalues) - if !all(isequal(Colon()), indices) - values = values[indices...] - end - - return values, Tuple(dims) -end - -# field and attributes from a dataset that implements PB.get_field(dataset, name)::FieldRecord -function _get_field_attributes(coord_source::Tuple{Any, Int}, name) - dataset, record_idx = coord_source - fr = PB.get_field(dataset, name) - return fr[record_idx], fr.attributes -end - -function _get_field_attributes(coord_source::Tuple{PB.ModelData, AbstractString}, varname) - modeldata, domainname = coord_source - # PB.get_field(model, modeldata, varnamefull) doesn't provide attributes, so do this ourselves - varnamefull = domainname*"."*varname - var = PB.get_variable(modeldata.model, varnamefull) - !isnothing(var) || - throw(ArgumentError("Variable $varnamefull not found")) - f = PB.get_field(var, modeldata) - attributes = copy(var.attributes) - attributes[:var_name] = var.name - attributes[:domain_name] = var.domain.name - - return f, attributes -end \ No newline at end of file