Merge pull request #107 from PALEOtoolkit/dims_coords

Simplify and align dimensions and coordinates with netcdf Common Data Model

Project.toml

@@ -1,20 +1,21 @@
 name = "PALEOmodel"
 uuid = "bf7b4fbe-ccb1-42c5-83c2-e6e9378b660c"
 authors = ["Stuart Daines <[email protected]>"]
-version = "0.15.49"
+version = "0.16.0"
 
 [deps]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DiffRules = "b552c78f-8df3-52c6-915a-8e097449b14b"
+DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b"
-JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 MultiFloats = "bdf0d083-296b-4888-a5b6-7498122e68a5"
 NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
+OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 PALEOboxes = "804b410e-d900-4b2a-9ecd-f5a06d4c1fd4"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
@@ -34,14 +35,15 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 [compat]
 DataFrames = "1.1"
 DiffRules = "1.0"
+DocStringExtensions = "0.8, 0.9"
 FileIO = "1.0"
 ForwardDiff = "0.10"
 Infiltrator = "1.0"
-JLD2 = "0.4, 0.5"
 MultiFloats = "1.0, 2.0"
 NCDatasets = "0.12, 0.14.2"
 NLsolve = "4.5"
-PALEOboxes = "0.21.37"
+OrderedCollections = "1.7.0"
+PALEOboxes = "0.22"
 RecipesBase = "1.2"
 Requires = "1.0"
 Revise = "3.1"

docs/src/PALEOmodelOutput.md

@@ -69,20 +69,19 @@ OutputMemoryDomain
 ```@docs
 save_netcdf
 load_netcdf!
-save_jld2
-load_jld2!
 ```
 
 ## Field Array
 
 ```@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)
+FixedCoord
 ```
 
 ## Plotting output

src/CoordsDims.jl

@@ -0,0 +1,105 @@
+
+################################
+# Coordinates
+#################################
+
+"""
+    FixedCoord(name, values::Array{Float64, N}, attributes::Dict)
+
+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)
+"""
+mutable struct FixedCoord{N}
+    name::String
+    values::Array{Float64, N}
+    attributes::Dict{Symbol, Any}
+end
+
+is_boundary_coordinate(fc::FixedCoord) = ndims(fc.values) == 2 && size(fc.values, 1) == 2
+
+
+
+##################################################################
+# Coordinate filtering and selection
+################################################################
+
+"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
+
+"""
+    dimscoord_subset(dim::PB.NamedDimension, coords::Vector{FixedCoord}, select_dimvals::AbstractString, select_filter)
+        -> cidx, dim_subset, coords_subset, coords_used
+
+Filter dimension `dim` according to key `select_dimvals` and `select_filter` (typically a single value or a range)
+
+Filtering may be applied either to dimension indices from `dim`, or to coordinates from `coords`:
+- If `select_dimvals` is of form "<dimname>_isel", use dimension indices and return `coords_used=nothing`
+- Otherwise use coordinate values and return actual coordinate values used in `coords_used`
+
+`cidx` are the filtered indices to use:
+- if `cidx` is a scalar (an Int), `dim_subset=nothing` and `coords_subset=nothing` indicating this dimension should be squeezed out
+- otherwise `cidx` is a Vector and `dim_subset` and `coords_subset` are the filtered subset of `dim` and `coords`
+"""
+function dimscoord_subset(dim::PB.NamedDimension, coords::Vector{FixedCoord}, select_dimvals::AbstractString, select_filter)
+    if length(select_dimvals) > 5 && select_dimvals[end-4:end] == "_isel"
+        @assert select_dimvals[1:end-5] == dim.name
+        cidx = select_filter
+        coords_used=nothing
+    else
+        # find cidx corresponding to a coordinate
+        # find coordinate to use in coords
+        ccidx = findfirst(c -> c.name == select_dimvals, coords)
+        @assert !isnothing(ccidx)
+        cc = coords[ccidx]
+        cidx, cvalue = find_indices(cc.values, select_filter)
+        # reset to the value actually used
+        coords_used = cvalue
+    end
+
+    if cidx isa AbstractVector
+        dim_subset = PB.NamedDimension(dim.name, length(cidx))
+        coords_subset = FixedCoord[]
+        for c in coords
+            if is_boundary_coordinate(c)
+                cs = FixedCoord(c.name, c.values[:, cidx], c.attributes)
+            else
+                cs = FixedCoord(c.name, c.values[cidx], c.attributes)
+            end
+            push!(coords_subset, cs)
+        end
+    else
+        # squeeze out dimensions
+        dim_subset = nothing
+        coords_subset = nothing
+    end
+
+    return cidx, dim_subset, coords_subset, coords_used
+end
+