From a2cc110b4b28e3a524953098a6b92dbb7520e235 Mon Sep 17 00:00:00 2001 From: Daines Date: Fri, 5 May 2023 11:44:22 +0100 Subject: [PATCH] Datamodel update test code See https://github.com/PALEOtoolkit/PALEOmodel.jl/issues/48 https://github.com/PALEOtoolkit/PALEOboxes.jl/pull/85 - Generalize FieldArray - move get_region and coordinate manipulation code from PALEOboxes.jl - rename FixedCoord -> RecordCoord as now only used by FieldRecord (and should be used by OutputWriter) TODO this currently allows a model (only the black sea config tested) to run and store output, but the FieldRecord etc (and surely FieldArray and plotting) does not work. Stop at this point as this illustrates a review of the data model for internal data storage (at least) is needed. Which requires implementing a standardised stable output format on disk first. --- src/FieldArray.jl | 349 ++++++++++++++++++++++++++++++++++++++----- src/FieldRecord.jl | 55 ++++++- src/OutputWriters.jl | 2 +- src/PlotRecipes.jl | 98 +++++++++++- 4 files changed, 452 insertions(+), 52 deletions(-) diff --git a/src/FieldArray.jl b/src/FieldArray.jl index 1edfec0..ce0b426 100644 --- a/src/FieldArray.jl +++ b/src/FieldArray.jl @@ -1,5 +1,10 @@ """ - FieldArray + FieldArray( + name::String, values::Array; + dims::Union{Nothing, Vector{<:AbstractString}, Vector{<:Pair}}=nothing, + coords::Vector{<:FieldArray}=FieldArray[], + attributes=nothing, + ) A generic [xarray](https://xarray.pydata.org/en/stable/index.html)-like or [IRIS](https://scitools-iris.readthedocs.io/en/latest/)-like @@ -7,25 +12,111 @@ 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. + +# Examples + +2-dimensional FieldArray with no supplied dimension names or coordinates: + + julia> fa = PALEOmodel.FieldArray("vec", [1 2 3; 4 5 6]) + FieldArray(name="vec", eltype=Int64) + dims: (NamedDimension(name=dim_1, size=2, coords=String[]), NamedDimension(name=dim_2, size=3, coords=String[])) + +2-dimensional FieldArray with supplied dimension names, no coordinates: + + julia> fa = PALEOmodel.FieldArray("vec", [1 2 3; 4 5 6]; dims=["x", "y"]) + FieldArray(name="vec", eltype=Int64) + dims: (NamedDimension(name=x, size=2, coords=String[]), NamedDimension(name=y, size=3, coords=String[])) + +2-dimensional FieldArray with supplied dimension names and coordinates: + + julia> fa = PALEOmodel.FieldArray( + "vec", [1 2 3; 4 5 6]; + dims=[ + "y"=>("ymid",), + "x"=>("xmid",) + ], + coords=[ + PALEOmodel.FieldArray("xmid", [0.5, 1.5, 2.5], dims=["x"]), + PALEOmodel.FieldArray("ymid", [0.5, 1.5], dims=["y"]) + ], + ) + + FieldArray(name="vec", eltype=Int64) + dims: (NamedDimension(name=y, size=2, coords=["ymid"]), NamedDimension(name=x, size=3, coords=["xmid"])) + coords: + FieldArray(name="xmid", eltype=Float64, dims=(NamedDimension(name=x, size=3, coords=String[]),)) + FieldArray(name="ymid", eltype=Float64, dims=(NamedDimension(name=y, size=2, coords=String[]),)) + """ struct FieldArray{T, N} name::String values::T dims::NTuple{N, PB.NamedDimension} + coords::Vector{FieldArray} attributes::Union{Dict, Nothing} end +function FieldArray( + name::String, values::Array; + dims::Union{Nothing, Vector{<:AbstractString}, Vector{<:Pair}}=nothing, + coords::Vector{<:FieldArray}=FieldArray[], + attributes=nothing, +) + if isnothing(dims) + dims = [PB.NamedDimension("dim_$i", s, []) for (i, s) in enumerate(size(values))] + + elseif dims isa Vector{<:AbstractString} + tmpdims = [] + for (d, s) in PB.IteratorUtils.zipstrict(dims, size(values); + errmsg="values has $(length(size(values))) dimensions but $(length(dims)) dims supplied") + push!(tmpdims, PB.NamedDimension(d, s, [])) + end + dims=tmpdims + elseif dims isa Vector{<:Pair} + tmpdims = [] + for ((dname, dcoords), s) in PB.IteratorUtils.zipstrict(dims, size(values); + errmsg="values has $(length(size(values))) dimensions but $(length(dims)) dims supplied") + push!(tmpdims, PB.NamedDimension(dname, s, [dc for dc in dcoords])) + end + dims=tmpdims + else + error("invalid dims: ", dims) + end + + # check supplied coords have same named dimensions of same size + tmpcoords = FieldArray[] # narrow type + for co in coords + for codim in co.dims + arraydimidx = findfirst(d->d.name == codim.name, dims) + !isnothing(arraydimidx) || error("coord '$(co.name)' dimension '$(codim.name)' not present in array dimensions $([d.name for d in dims])") + arraydim = dims[arraydimidx] + codim.size == arraydim.size || error("coord '$(co.name)' dimension '$(codim.name)': size $(codim.size) != array dimension size $(arraydim.size)") + end + push!(tmpcoords, co) + end + + return FieldArray(name, values, Tuple(dims), tmpcoords, attributes) +end + + function Base.show(io::IO, fa::FieldArray) print(io, - "FieldArray(name=\"", fa.name, "\", eltype=", eltype(fa.values),", size=", size(fa.values), ", dims=", fa.dims, ")" + "FieldArray(name=\"", fa.name, "\", eltype=", eltype(fa.values), ", dims=", fa.dims, ")" ) end function Base.show(io::IO, ::MIME"text/plain", fa::FieldArray) println(io, "FieldArray(name=\"", fa.name, "\", eltype=", eltype(fa.values),")") println(io, " dims: ", fa.dims) - println(io, " attributes: ", fa.attributes) + + if !isempty(fa.coords) + println(io, " coords:") + for co in fa.coords + println(io, " ", co) + end + end + # println(io, " attributes: ", fa.attributes) end # basic arithmetic operations @@ -52,6 +143,84 @@ function default_fieldarray_name(attributes::Dict) return name end +""" + get_dimension(fa::FieldArray, dimname::AbstractString) -> nd::NamedDimension + +Get dimension `dimname` (error if not present) +""" +function get_dimension(fa::FieldArray, dimname::AbstractString) + arraydimidx = findfirst(d->d.name == dimname, fa.dims) + !isnothing(arraydimidx) || error("dimension '$dimname' not present") + return dims[arraydimidx] +end + +""" + get_coord_arrays(fa::FieldArray, dimname::AbstractString) -> coord_arrays::Vector{FieldArray} + +Get coordinates for dimension `dimname` (empty Vector if no coordinates defined) +""" +function get_coord_arrays(fa::FieldArray, dimname::AbstractString) + + nd = get_dimension(fa, dimname) + coord_names = nd.coords + coord_arrays = FieldArray[] + for cname in coord_names + cnidx = findfirst(c->c.name == cname, fa.coords) + !isnothing(cnidx) || error("dimension '$dimname' coordinate '$cname' not present") + push!(coord_arrays, fa.coords[cnidx]) + end + + return coord_arrays +end + +""" + 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 + 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 + +# 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}) + ) + ) + ############################################################# # Create from PALEO objects ############################################################# @@ -184,50 +353,154 @@ function get_array( 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}) - ) - ) + + + + +######################################### +# get_region moved from PALEOboxes.jl +######################################### """ - update_coordinates(varray::FieldArray, vec_coord_arrays::AbstractVector) -> FieldArray + get_region(grid::Union{PB.AbstractMesh, Nothing}, values; selectargs...) -> values_subset, (dim_subset::NamedDimension, ...) -Replace or add coordinates `vec_coord_arrays` to `varray`. +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 -`new_coord_arrays` is a Vector of Pairs of "coord_name"=>(var1::FieldArray, var2::FieldArray, ...) +""" + get_region(grid::Nothing, values) -> values[] -Example: to replace a 1D column default pressure coordinate with a z coordinate: - - coords=["z"=>(zmid::FieldArray, zlower::FieldArray, atm.zupper::FieldArray)] +Fallback for Domain with no grid, assumed 1 cell """ -function update_coordinates(varray::FieldArray, vec_coord_arrays::AbstractVector) +function get_region(grid::Nothing, values) + length(values) == 1 || + throw(ArgumentError("grid==Nothing and length(values) != 1")) + return values[], () +end + +""" + get_region(grid::UnstructuredVectorGrid, values; cell) -> + values_subset, (dim_subset::NamedDimension, ...) - 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)]") +# Keywords for region selection: +- `cell::Union{Int, Symbol}`: an Int, or a Symbol to look up in `cellnames` +""" +function get_region(grid::PB.Grids.UnstructuredVectorGrid, values; cell::Union{Int, Symbol}) + if 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 - # 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, - ) + return ( + values[idx], + (), # no dimensions (ie squeeze out a dimension length 1 for single cell) + ) +end + +""" + get_region(grid::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) + + 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] + return ( + values[indices], + (PB.NamedDimension("z", length(indices), PB.get_region(grid.z_coords, indices)), ), + ) + else + # squeeze out z dimension + idx = grid.Icolumns[colidx][cell] + return ( + values[idx], + (), # no dimensions (ie squeeze out a dimension length 1 for single cell) ) end + +end - # replace coordinates - varray_newcoords = FieldArray(varray.name, varray.values, Tuple(named_dimensions), varray.attributes) +""" + get_region(grid::Union{CartesianLinearGrid{2}, CartesianArrayGrid{2}} , internalvalues; [i=i_idx], [j=j_idx]) -> + arrayvalues_subset, (dim_subset::NamedDimension, ...) - return varray_newcoords -end \ No newline at end of file +# 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() +) + return _get_region(grid, internalvalues, [i, j]) +end + +""" + get_region(grid::Union{CartesianLinearGrid{3}, 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() +) + return _get_region(grid, internalvalues, [i, j, k]) +end + +function _get_region( + grid::Union{PB.Grids.CartesianLinearGrid, PB.Grids.CartesianArrayGrid}, internalvalues, indices +) + 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 diff --git a/src/FieldRecord.jl b/src/FieldRecord.jl index 16d5e77..bdb6d84 100644 --- a/src/FieldRecord.jl +++ b/src/FieldRecord.jl @@ -1,3 +1,20 @@ +mutable struct RecordCoord + name::String + values::Vector{Float64} + attributes::Dict{Symbol, Any} +end + + + +function get_region(fc::RecordCoord, indices::AbstractVector) + return FixedCoord(fc.name, fc.values[indices], fc.attributes) +end + +function get_region(fcv::Vector{RecordCoord}, indices::AbstractVector) + return [RecordCoord(fc.name, fc.values[indices], fc.attributes) for fc in fcv] +end + + """ FieldRecord{D <: AbstractData, S <: AbstractSpace, V, N, M, R} FieldRecord( @@ -17,15 +34,17 @@ Fields with single `values` (`field_single_element` true) are stored as a Vector struct FieldRecord{D <: PB.AbstractData, S <: PB.AbstractSpace, V, N, M, R} records::Vector{R} data_dims::NTuple{N, PB.NamedDimension} + record_dim::PB.NamedDimension mesh::M + coords::Vector{FieldArray} attributes::Dict{Symbol, Any} - coords_record::Vector{PB.FixedCoord} # coordinates attached to record dimension + # coords_record::Vector{RecordCoord} # coordinates attached to record dimension end function Base.show(io::IO, fr::FieldRecord) print(io, - "FieldRecord(eltype=", eltype(fr),", length=", length(fr), - ", attributes=", fr.attributes, ", coords_record=", fr.coords_record, ")" + "FieldRecord(eltype=", eltype(fr),", length=", length(fr), ", record_dim=", fr.record_dim, + ", attributes=", fr.attributes, ")" ) end @@ -223,7 +242,7 @@ function _get_array( # find ridx corresponding to a coordinate for cr in fr.coords_record if String(k) == cr.name - ridx, cvalue = PB.find_indices(cr.values, v) + ridx, cvalue = find_indices(cr.values, v) selectargs_records=NamedTuple((k=>cvalue,)) end end @@ -231,7 +250,8 @@ function _get_array( selectargs_region = merge(selectargs_region, NamedTuple((k=>v,))) end end - records_dim = PB.NamedDimension("records", length(ridx), PB.get_region(fr.coords_record, ridx)) + records_coord = + records_dim = PB.NamedDimension("records", length(ridx), get_region(fr.coords_record, ridx)) # add attributes for selection used attributes = copy(fr.attributes) @@ -312,3 +332,28 @@ function _get_array( end +"find indices of coord from first before range[1] to first after range[2]" +function find_indices(coord::AbstractVector, range) + length(range) == 2 || + throw(ArgumentError("find_indices: length(range) != 2 $range")) + + idxstart = findlast(t -> t<=range[1], coord) + isnothing(idxstart) && (idxstart = 1) + + idxend = findfirst(t -> t>=range[2], coord) + isnothing(idxend) && (idxend = length(coord)) + + return idxstart:idxend, (coord[idxstart], coord[idxend]) +end + +"find indices of coord nearest val" +function find_indices(coord::AbstractVector, val::Real) + idx = 1 + for i in 1:length(coord) + if abs(coord[i] - val) < abs(coord[idx] - val) + idx = i + end + end + + return [idx], coord[idx] +end \ No newline at end of file diff --git a/src/OutputWriters.jl b/src/OutputWriters.jl index 9302709..50fcb04 100644 --- a/src/OutputWriters.jl +++ b/src/OutputWriters.jl @@ -440,7 +440,7 @@ function PB.get_field(odom::OutputMemoryDomain, varname) fr = PALEOmodel.wrap_fieldrecord( vdata, field_data, Tuple(data_dims), missing, space, grid, attributes, coords_record=[ - PB.FixedCoord( + PALEOmodel.RecordCoord( String(odom.coords_record), # df[!, odom.coords_record], use_all_records ? df[!, odom.coords_record] : df[1:odom._nrecs, odom.coords_record], diff --git a/src/PlotRecipes.jl b/src/PlotRecipes.jl index e3cc46c..b2958e7 100644 --- a/src/PlotRecipes.jl +++ b/src/PlotRecipes.jl @@ -364,6 +364,12 @@ RecipesBase.@recipe function f( end 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 values = fa.values name = fa.name @@ -386,20 +392,20 @@ RecipesBase.@recipe function f( else label --> labelprefix*name end - f_label = PB.append_units(get_attribute(fa, :var_name, ""), fa.attributes) + f_label = append_units(get_attribute(fa, :var_name, ""), fa.attributes) coords_vec = fa.dims[1].coords if length(coords_vec) == 1 || length(coords_vec) > 3 co = first(coords_vec) return swap_coords_values( co.values, values, - PB.append_units(co.name, co.attributes), f_label, + append_units(co.name, co.attributes), f_label, ) elseif length(coords_vec) in (2, 3) co_lower = coords_vec[end-1] co_upper = coords_vec[end] - c_label = PB.append_units(co_lower.name*", "*co_upper.name, co_lower.attributes) + c_label = append_units(co_lower.name*", "*co_upper.name, co_lower.attributes) return swap_coords_values( create_stepped(co_lower.values, co_upper.values, values)..., c_label, f_label, @@ -412,19 +418,19 @@ RecipesBase.@recipe function f( end elseif length(fa.dims) == 2 - title --> PB.append_units(fa.name, fa.attributes) + title --> append_units(fa.name, fa.attributes) - i_values, i_label = PB.build_coords_edges(fa.dims[1]) - j_values, j_label = PB.build_coords_edges(fa.dims[2]) + i_values, i_label = build_coords_edges(fa, fa.dims[1].name) + j_values, j_label = build_coords_edges(fa, fa.dims[2].name) # workaround for Plots.jl heatmap: needs x, y to either both be midpoints, or both be edges if length(i_values) == size(values, 1) && length(j_values) == size(values, 2)+1 @warn "$(fa.name) guessing '$i_label' coordinate edges from midpoints assuming uniform spacing" - i_values = PB.guess_coords_edges(i_values) + i_values = guess_coords_edges(i_values) end if length(i_values) == size(values, 1)+1 && length(j_values) == size(values, 2) @warn "$(fa.name) guessing '$j_label' coordinate edges from midpoints assuming uniform spacing" - j_values = PB.guess_coords_edges(j_values) + j_values = guess_coords_edges(j_values) end x_values, y_values, f_values = swap_coords_xy( @@ -527,3 +533,79 @@ function broadcast_dict(dv::Vector{<:Dict}) return dvout end 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::PB.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 + +""" + 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(fa::FieldArray, dimname::AbstractString) + + coords_vec = get_coord_arrays(fa, dimname) + + if isempty(coords_vec) + nd = get_dimension(fa, dimname) + @warn "no coords for dimension $(nd), returning indices" + return collect(1:nd.size), nd.name*" (indices)" + 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) 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 + +""" + 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