Add time and date bounds

Author: ph-kev

NEWS.md

@@ -1,4 +1,18 @@
 # NEWS
+v0.2.13
+-------
+
+## Features
+
+### Add time and date bounds for netCDF files
+In a netCDF file produced by ClimaDiagnostics, there are now the dimensions
+`time_bnds` and `date_bnds`. Each time or date value is a representative of the
+corresponding time or date bound. For example, if the time is `10.0` and
+corresponding time bound is `[0.0, 10.0]`, then the time of `10.0` represents
+the interval `[0.0, 10.0]`. If one knows that the data represents a time
+average, then the time of `10.0` is the time average over the interval
+`[0.0, 10.0]`.
+
 v0.2.12
 -------
 

src/netcdf_writer.jl

@@ -355,6 +355,7 @@ function write_field!(writer::NetCDFWriter, field, diagnostic, u, p, t)
         axis = "T",
         standard_name = "time",
         long_name = "Time",
+        bounds = "time_bnds",
     )
 
     dim_names = add_space_coordinates_maybe!(
@@ -373,12 +374,25 @@ function write_field!(writer::NetCDFWriter, field, diagnostic, u, p, t)
     else
         start_date = writer.start_date
     end
+
+    add_time_bounds_maybe!(
+        nc,
+        TT;
+        comments = "time bounds for each time value",
+        units = "s",
+    )
+
     if !isnothing(start_date)
         add_date_maybe!(
             nc;
             units = "seconds since $start_date",
             bounds = "date_bnds",
         )
+        add_date_bounds_maybe!(
+            nc;
+            comments = "date bounds for each date value",
+            units = "seconds since $start_date",
+        )
     end
 
     if haskey(nc, "$(var.short_name)")
@@ -411,6 +425,9 @@ function write_field!(writer::NetCDFWriter, field, diagnostic, u, p, t)
 
     # TODO: Use ITime here
     nc["time"][time_index] = float(t)
+    nc["time_bnds"][:, time_index] =
+        time_index == 1 ? [zero(float(t)); float(t)] :
+        [nc["time"][time_index - 1]; float(t)]
 
     # FIXME: We are hardcoding p.start_date !
     # FIXME: We are rounding t
@@ -419,6 +436,9 @@ function write_field!(writer::NetCDFWriter, field, diagnostic, u, p, t)
         curr_date = start_date + Dates.Millisecond(round(1000 * float(t)))
         date_type = typeof(curr_date) # not necessarily a Dates.DateTime
         nc["date"][time_index] = curr_date
+        nc["date_bnds"][:, time_index] =
+            time_index == 1 ? [start_date; curr_date] :
+            [date_type(nc["date"][time_index - 1]); curr_date]
     end
 
     # TODO: It would be nice to find a cleaner way to do this

src/netcdf_writer_coordinates.jl

@@ -132,6 +132,49 @@ function add_date_maybe!(nc::NCDatasets.NCDataset; kwargs...)
         dim.attrib[String(k)] = v
     end
 end
+
+"""
+    add_time_bounds_maybe!(nc::NCDatasets.NCDataset,
+                           float_type::Type{FT};
+                           kwargs...)
+
+Add the `time_bnds` dimension (with infinite size) to the given NetCDF file if
+not already there. Optionally, add all the keyword arguments as attributes.
+
+This function is meant to be called after `add_time_maybe!`.
+"""
+function add_time_bounds_maybe!(
+    nc::NCDatasets.NCDataset,
+    float_type::Type{FT};
+    kwargs...,
+) where {FT}
+    haskey(nc, "time_bnds") && return nothing
+    !haskey(nc, "time") && error(
+        "Time dimension does not exist. Call add_time_maybe! before calling this function",
+    )
+    !("nv" in NCDatasets.dimnames(nc)) && NCDatasets.defDim(nc, "nv", 2) # number of vertices
+    dim = NCDatasets.defVar(nc, "time_bnds", FT, ("nv", "time"))
+    for (k, v) in kwargs
+        dim.attrib[String(k)] = v
+    end
+    return nothing
+end
+
+"""
+    add_date_bounds_maybe!(nc::NCDatasets.NCDataset; kwargs...)
+
+Add the `date_bnds` dimension (with infinite size) to the given NetCDF file if
+not already there. Optionally, add all the keyword arguments as attributes.
+
+This function is meant to be called after `add_time_maybe!`.
+"""
+function add_date_bounds_maybe!(nc::NCDatasets.NCDataset; kwargs...)
+    haskey(nc, "date_bnds") && return nothing
+    !haskey(nc, "date") && error(
+        "Time dimension does not exist. Call add_time_maybe! before calling this function",
+    )
+    !("nv" in NCDatasets.dimnames(nc)) && NCDatasets.defDim(nc, "nv", 2) # number of vertices
+    dim = NCDatasets.defVar(nc, "date_bnds", Float64, ("nv", "time"))
     for (k, v) in kwargs
         dim.attrib[String(k)] = v
     end

test/writers.jl

@@ -107,6 +107,13 @@ end
         @test nc["lat"].attrib["standard_name"] == "latitude"
         @test nc["lat"].attrib["long_name"] == "Latitude"
         @test nc["lat"].attrib["axis"] == "Y"
+
+        # Test bounds
+        @test nc["time_bnds"][:, 1] == [0.0; 10.0]
+        @test nc["date_bnds"][:, 1] == [
+            Dates.DateTime(1453, 5, 29)
+            Dates.DateTime(1453, 5, 29) + Dates.Second(10.0)
+        ]
     end
 
     # Disable vertical interpolation