read/write axis and standard_name for netcdf (#331)

* detect axis name for netcdf

* also detect standard_name

* bugfix

Author: rafaqz

src/sources/ncdatasets.jl

@@ -6,7 +6,7 @@ const NCDAllowedType = Union{Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64,F
 
 # CF standards don't enforce dimension names.
 # But these are common, and should take care of most dims.
-const NCD_DIMMAP = Dict(
+const NCD_DIM_MAP = Dict(
     "lat" => Y,
     "latitude" => Y,
     "lon" => X,
@@ -23,6 +23,20 @@ const NCD_DIMMAP = Dict(
     "band" => Band,
 )
 
+const NCD_AXIS_MAP = Dict(
+    "X" => X,
+    "Y" => Y,
+    "Z" => Z,
+    "T" => Ti,
+)
+
+const NCD_STANDARD_NAME_MAP = Dict(
+    "longitude" => X,
+    "latitude" => Y,
+    "depth" => Z,
+    "time" => Ti,
+)
+
 haslayers(::Type{NCDfile}) = true
 defaultcrs(::Type{NCDfile}) = EPSG(4326)
 defaultmappedcrs(::Type{NCDfile}) = EPSG(4326)
@@ -60,7 +74,7 @@ Write an NCDarray to a NetCDF file using NCDatasets.jl
 
 Returns `filename`.
 """
-function Base.write(filename::AbstractString, ::Type{NCDfile}, A::AbstractRaster; append = false, kw...)
+function Base.write(filename::AbstractString, ::Type{NCDfile}, A::AbstractRaster; append=false, kw...)
     mode  = !isfile(filename) || !append ? "c" : "a";
     ds = NCD.Dataset(filename, mode; attrib=_attribdict(metadata(A)))
     try
@@ -161,7 +175,7 @@ function DD.layerdims(ds::NCD.Dataset)
 end
 function DD.layerdims(var::NCD.Variable)
     map(NCD.dimnames(var)) do dimname
-        _ncddimtype(dimname)()
+        _ncddimtype(var.attrib, dimname)()
     end
 end
 
@@ -228,15 +242,17 @@ cleanreturn(A::NCD.CFVariable) = Array(A)
 
 function _ncddim(ds, dimname::Key, crs=nothing, mappedcrs=nothing)
     if haskey(ds, dimname)
-        D = _ncddimtype(dimname)
+        var = NCD.variable(ds, dimname)
+        D = _ncddimtype(var.attrib, dimname)
         lookup = _ncdlookup(ds, dimname, D, crs, mappedcrs)
         return D(lookup)
     else
         # The var doesn't exist. Maybe its `complex` or some other marker,
         # so make it a custom `Dim` with `NoLookup`
         len = _ncfinddimlen(ds, dimname)
         len === nothing && _unuseddimerror()
-        return Dim{Symbol(dimname)}(NoLookup(Base.OneTo(len)))
+        lookup = NoLookup(Base.OneTo(len))
+        return Dim{Symbol(dimname)}(lookup)
     end
 end
 
@@ -253,25 +269,42 @@ end
 
 # Find the matching dimension constructor. If its an unknown name
 # use the generic Dim with the dim name as type parameter
-_ncddimtype(dimname) = haskey(NCD_DIMMAP, dimname) ? NCD_DIMMAP[dimname] : DD.basetypeof(DD.key2dim(Symbol(dimname)))
+function _ncddimtype(attrib, dimname)
+    if haskey(attrib, "axis") 
+        k = attrib["axis"] 
+        if haskey(attrib, k) 
+            return NCD_AXIS_MAP[k] 
+        end
+    end
+    if haskey(attrib, "standard_name")
+        k = attrib["standard_name"]
+        if haskey(NCD_STANDARD_NAME_MAP, k) 
+            return NCD_STANDARD_NAME_MAP[k]
+        end
+    end
+    if haskey(NCD_DIM_MAP, dimname) 
+        return NCD_DIM_MAP[dimname] 
+    end
+    return DD.basetypeof(DD.key2dim(Symbol(dimname)))
+end
 
 # _ncdlookup
 # Generate a `LookupArray` from a netcdf dim.
-function _ncdlookup(ds::NCD.Dataset, dimname, D, crs, mappedcrs)
+function _ncdlookup(ds::NCD.Dataset, dimname, D::Type, crs, mappedcrs)
     dvar = ds[dimname]
     index = dvar[:]
     metadata = _metadatadict(NCDfile, dvar.attrib)
     return _ncdlookup(ds, dimname, D, index, metadata, crs, mappedcrs)
 end
 # For unknown types we just make a Categorical lookup
-function _ncdlookup(ds::NCD.Dataset, dimname, D, index::AbstractArray, metadata, crs, mappedcrs)
+function _ncdlookup(ds::NCD.Dataset, dimname, D::Type, index::AbstractArray, metadata, crs, mappedcrs)
     Categorical(index; order=Unordered(), metadata=metadata)
 end
 # For Number and AbstractTime we generate order/span/sampling
 # We need to include `Missing` in unions in case `_FillValue` is used
 # on coordinate variables in a file and propagates here.
 function _ncdlookup(
-    ds::NCD.Dataset, dimname, D, index::AbstractArray{<:Union{Missing,Number,Dates.AbstractTime}},
+    ds::NCD.Dataset, dimname, D::Type, index::AbstractArray{<:Union{Missing,Number,Dates.AbstractTime}},
     metadata, crs, mappedcrs
 )
     # Assume the locus is at the center of the cell if boundaries aren't provided.
@@ -426,7 +459,10 @@ function _def_dim_var!(ds::NCD.Dataset, dim::Dimension)
     if dim isa Y || dim isa X
         dim = convertlookup(Mapped, dim)
     end
+    # Attributes
     attrib = _attribdict(metadata(dim))
+    _ncd_set_axis_attrib!(attrib, dim)
+    # Bounds variables
     if span(dim) isa Explicit
         bounds = val(span(dim))
         boundskey = get(metadata(dim), :bounds, string(dimkey, "_bnds"))
@@ -437,6 +473,15 @@ function _def_dim_var!(ds::NCD.Dataset, dim::Dimension)
     return nothing
 end
 
+# Add axis and standard name attributes to dimension variabls
+# We need to get better at guaranteeing if X/Y is actually measured in `longitude/latitude`
+# CF standards requires that we specify "units" if we use these standard names
+_ncd_set_axis_attrib!(atr, dim::X) = atr["axis"] = "X" # at["standard_name"] = "longitude";
+_ncd_set_axis_attrib!(atr, dim::Y) = atr["axis"] = "Y" # at["standard_name"] = "latitude"; 
+_ncd_set_axis_attrib!(atr, dim::Z) = (atr["axis"] = "Z"; atr["standard_name"] = "depth")
+_ncd_set_axis_attrib!(atr, dim::Ti) = (atr["axis"] = "T"; atr["standard_name"] = "time")
+_ncd_set_axis_attrib!(atr, dim) = nothing
+
 _ncdshiftlocus(dim::Dimension) = _ncdshiftlocus(lookup(dim), dim)
 _ncdshiftlocus(::LookupArray, dim::Dimension) = dim
 function _ncdshiftlocus(lookup::AbstractSampled, dim::Dimension)

test/sources/ncdatasets.jl

@@ -83,8 +83,6 @@ stackkeys = (
         @test typeof(lookup(ncarray)) <: Tuple{<:Mapped,<:Mapped,<:Sampled}
         @test bounds(ncarray) == ((0.0, 360.0), (-80.0, 90.0), (DateTime360Day(2001, 1, 1), DateTime360Day(2003, 1, 1)))
     end
-    tempfile = tempname() * ".nc"
-
 
     @testset "other fields" begin
         @test ismissing(missingval(ncarray))
@@ -211,6 +209,11 @@ stackkeys = (
             @test size(geoA) == size(ncarray)
             filename = tempname() * ".nc"
             write(filename, geoA)
+            @testset "CF attributes" begin
+                @test NCDatasets.Dataset(filename)[:x].attrib["axis"] == "X"
+                @test NCDatasets.Dataset(filename)[:x].attrib["bounds"] == "x_bnds"
+                # TODO  better units and standard name handling
+            end
             saved = read(Raster(filename))
             @test size(saved) == size(geoA)
             @test refdims(saved) == refdims(geoA)