Field and Grid include dimension

src/TMI.jl

@@ -59,7 +59,7 @@ export config, config_from_mat, config_from_nc,
     wetlocation, iswet,
     control2state, control2state!,
     surfacecontrol2field, surfacecontrol2field!,
-    sparsedatamap, config2nc, gridprops,
+    sparsedatamap, config2nc, axislabels,
     matrix_zyx2xyz,
     surface_oxygensaturation, oxygen, location_obs,
     zerosurfaceboundary,

src/config.jl

@@ -211,14 +211,14 @@ function cartesianindex(file::String)
 end
 
 """
-    function gridprops(file)
+    function axislabels(file)
     Read and assemble the grid properties.
 # Arguments
 - `file`: TMI NetCDF file name
 # Output
 - `grid`: TMI grid coordinates
 """
-function gridprops(file)
+function axislabels(file)
     if file[end-1:end] == "nc"
 
         lon = convert(Vector{Float64},ncread(file,"lon"))

src/field.jl

@@ -14,29 +14,14 @@
     longname::String
     units::String
 """
-struct Field{T}
-    tracer::Array{T,3}
-    γ::Grid
+struct Field{T,N}
+    tracer::Array{T,N}
+    γ::Grid{T,N}
     name::Symbol
     longname::String
     units::String
 end
 
-# """
-#     function Field(tracer::Array{T,3},γ::Grid) where T <: Real
-
-#     Outer constructor for Field if there's no worry about
-#     tracer type, long name, or units.
-# # Arguments
-# - `tracer::Array{T,3}`
-# - `γ::Grid`
-# # Output
-# - `field::Field`
-# """
-#function Field(tracer::Array{T,3},γ::Grid) where T <: Real
-#   return Field(tracer,γ,:none,"unknown","unknown")
-#end
-
 function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, x::Field)
     summary(io, x); println(io)
     print(io, "Field size ")
@@ -121,13 +106,12 @@ end
 
     read MATLAB field and transfer zyx format to xyz
 """
-function readfield(file,tracername,γ::Grid) 
+function readfield(file,tracername,γ::Grid{T,N}) where {T,N} 
 
     # The mode "r" stands for read-only. The mode "r" is the default mode and the parameter can be omitted.
     tracer, units, longname = _read3d(file,tracername)
     checkgrid!(tracer,γ.wet)
-
-    c = Field(tracer,γ,tracerdict()[tracername],longname,units)
+    c = Field{T,N}(tracer,γ,tracerdict()[tracername],longname,units)
     return c
 end
 function readfield(matfile,mattracername,γ::Grid,Izyx) 

src/grid.jl

@@ -1,14 +1,16 @@
+using Base: index_ndims, methods_including_ambiguous
 """
     struct Grid
 
     TMI grid with accounting for wet/dry points
 """
-struct Grid{T <: Real}
-    lon::Vector{T}
-    lat::Vector{T}
-    depth::Vector{T}
-    wet::BitArray{3}
-    interior::BitArray{3}
+struct Grid{T,N}
+    # lon::Vector{T}
+    # lat::Vector{T}
+    # depth::Vector{T}
+    axes::NTuple{N,Vector{T}}
+    wet::BitArray{N}
+    interior::BitArray{N}
 end
 
 """
@@ -23,16 +25,23 @@ function Grid(TMIfile)
 - `γ::Grid`: TMI grid struct
 """
 function Grid(TMIfile::String; A = watermassmatrix(TMIfile))
+
     # get properties of grid
-    lon,lat,depth = gridprops(TMIfile)
+    labels = axislabels(TMIfile)
+
+    ndims = length(labels)
+    #dimsize = Tuple{Int64}(undef,ndims)
+    dimsize = Tuple([length(labels[d]) for d in 1:ndims])
+    #end
 
     # make ocean mask
-    wet = wetmask(TMIfile,length(lon),length(lat),length(depth))
+    #wet = wetmask(TMIfile,length(lon),length(lat),length(depth))
+    wet = wetmask(TMIfile,dimsize)
 
     # make interior mask
-    interior = interiormask(A,wet,length(lon),length(lat),length(depth))
-
-    return Grid(lon,lat,depth,wet,interior)
+    #interior = interiormask(A,wet,length(lon),length(lat),length(depth))
+    interior = interiormask(A,wet,dimsize)
+    return Grid(labels,wet,interior)
 end
 
 """
@@ -70,7 +79,7 @@ function Grid(foreign_file::S, maskname::S, lonname::S, latname::S, depthname::S
     interior = deepcopy(wet)
     interior[:,:,1] .= false
 
-    return Grid(lon,lat,depth,wet,interior)
+    return Grid((lon,lat,depth),wet,interior)
 end
 
 """
@@ -79,12 +88,18 @@ end
     Do not store Cartesian and linear indices.
     Compute them on demand.
 """ 
-Base.propertynames(γ::Grid) = (:I,:R,fieldnames(typeof(γ))...)
+Base.propertynames(γ::Grid) = (:I,:R,:lon,:lat,:depth,fieldnames(typeof(γ))...)
 function Base.getproperty(γ::Grid, d::Symbol)
     if d === :I
         return cartesianindex(γ.wet)
     elseif d === :R
         return linearindex(γ.wet)
+    elseif d === :lon
+        return γ.axes[1]
+    elseif d === :lat
+        return γ.axes[2]
+    elseif d === :depth
+        return γ.axes[3]
     else
         return getfield(γ, d)
     end
@@ -137,12 +152,17 @@ function checkgrid!(tracer,wet)
     end
 end
 
+function wetmask(TMIfile::String,dimsize::NTuple)
+    # read Cartesian Index from file.
+    I = cartesianindex(TMIfile)
+    wet = falses(dimsize)
+    wet[I] .= 1
+    return wet
+end
 function wetmask(TMIfile::String,nx,ny,nz)
+    # older version that assumes N = 3
     # read Cartesian Index from file.
     I = cartesianindex(TMIfile)
-    # make a mask
-    # first choice: smaller but inconsistent with input grid
-    #wet = falses(maximum(I)[1],maximum(I)[2],maximum(I)[3])
     wet = falses(nx,ny,nz)
     wet[I] .= 1
     return wet
@@ -158,3 +178,10 @@ function interiormask(A,wet,nx,ny,nz)
     interior[I[list]] .= true 
     return interior
 end
+function interiormask(A,wet,dimsize::NTuple)
+    interior = falses(dimsize)
+    I = cartesianindex(wet)
+    list = findall(.!(isone.(sum(abs.(A),dims=2))))
+    interior[I[list]] .= true 
+    return interior
+end