fix timestamp to address phase offset (#1332)

kmdeck · web-flow · commit 1a412ea07644 · 2025-08-20T15:02:58.000-07:00
diff --git a/docs/src/tutorials/calibration/obs_site_level_calibration.jl b/docs/src/tutorials/calibration/obs_site_level_calibration.jl
@@ -70,8 +70,9 @@ site_ID_val = FluxnetSimulations.replace_hyphen(site_ID);
     FluxnetSimulations.get_location(FT, Val(site_ID_val))
 (; atmos_h) = FluxnetSimulations.get_fluxtower_height(FT, Val(site_ID_val));
 
-# Get maximum simulation start and end dates in UTC; these must be included in the forcing data range
-start_date = DateTime(2010, 3, 1)
+# Get simulation start and stop dates in UTC; these must be included in the forcing data range
+# Here we calibrate with only two months of data.
+start_date = DateTime(2010, 5, 1)
 stop_date = DateTime(2010, 7, 1)
 Δt = 450.0; # seconds
 
@@ -185,14 +186,16 @@ function model(Vcmax25, g1)
         land_model,
     )
 
-    #md # Configure diagnostics to output sensible and latent heat fluxes hourly
+    #md # Configure diagnostics to output sensible and latent heat fluxes half-hourly
+    #md # Since fluxnet data is recorded every half hour, with the average of the half hour,
+    #md # we need to be consistent here.
     output_vars = ["lhf"]
     diagnostics = ClimaLand.default_diagnostics(
         land_model,
         start_date;
         output_writer = ClimaDiagnostics.Writers.DictWriter(),
         output_vars,
-        average_period = :hourly,
+        average_period = :halfhourly,
     )
 
     #md # Create and run the simulation
@@ -238,7 +241,7 @@ end;
 function get_lhf(simulation)
     return ClimaLand.Diagnostics.diagnostic_as_vectors(
         simulation.diagnostics[1].output_writer,
-        "lhf_1h_average",
+        "lhf_30m_average",
     )
 end;
 
@@ -402,3 +405,9 @@ axislegend(
 CairoMakie.resize_to_layout!(fig)
 CairoMakie.save("G_first_and_last.png", fig)
 fig
+# Here we can see that the calibration has improved the fit to the data;
+# but this it not the fully story. If other parameters have been set to
+# unrealistic values, the parameters here may be compensating to
+# achieve a lower loss. Moreover, we have only calibrated with two months
+# of data. Production-level calibration runs require careful choice
+# of free parameters and target observations.
diff --git a/docs/src/tutorials/integrated/snowy_land_fluxnet_tutorial.jl b/docs/src/tutorials/integrated/snowy_land_fluxnet_tutorial.jl
@@ -103,6 +103,8 @@ diagnostics = ClimaLand.default_diagnostics(
 );
 
 # Choose how often we want to update the forcing.
+# Choosing a frequency > the data frequency results in linear
+# interpolation in time to the intermediate times.
 data_dt = Second(FluxnetSimulations.get_data_dt(site_ID));
 updateat = Array(start_date:data_dt:stop_date);
 
diff --git a/docs/src/tutorials/standalone/Canopy/canopy_tutorial.jl b/docs/src/tutorials/standalone/Canopy/canopy_tutorial.jl
@@ -65,27 +65,26 @@ import ClimaLand.LandSimVis as LandSimVis
 const FT = Float32;
 earth_param_set = LP.LandParameters(FT);
 
-# First provide some information about the site:
-# Timezone (offset from UTC in hrs)
-time_offset = 7
-start_date = DateTime(2010) + Hour(time_offset)
-
-# Select a time range to perform time stepping over, and a dt. As usual,
-# the timestep depends on the problem you are solving, the accuracy of the
-# solution required, and the timestepping algorithm you are using.
-N_days = 364
-stop_date = start_date + Day(N_days)
-dt = 225.0
+# We will use prescribed atmospheric and radiative forcing from the
+# US-MOz tower.
+site_ID = "US-MOz";
+site_ID_val = FluxnetSimulations.replace_hyphen(site_ID);
+# Get the latitude and longitude in degrees, as well as the
+# time offset in hours of local time from UTC
+(; time_offset, lat, long) =
+    FluxnetSimulations.get_location(FT, Val(site_ID_val));
+# Get the height of the sensors in m
+(; atmos_h) = FluxnetSimulations.get_fluxtower_height(FT, Val(site_ID_val));
+# Set a start and stop date of the simulation in UTC, as well as
+# a timestep in seconds
+start_date = DateTime("2010-05-01", "yyyy-mm-dd")
+stop_date = DateTime("2010-09-01", "yyyy-mm-dd")
+dt = 450.0
 
 # Site latitude and longitude
 lat = FT(38.7441) # degree
 long = FT(-92.2000) # degree
 
-# Height of the sensor at the site
-atmos_h = FT(32)
-# Site ID
-site_ID = "US-MOz";
-
 # # Setup the Canopy Model
 
 # We want to simulate a vegetative canopy in standalone mode, without coupling
@@ -103,8 +102,6 @@ site_ID = "US-MOz";
 # would change.
 domain = Point(; z_sfc = FT(0.0), longlat = (long, lat));
 
-
-
 # We will be using prescribed atmospheric and radiative drivers from the
 # US-MOz tower, which we read in here. We are using prescribed
 # atmospheric and radiative flux conditions, but it is also possible to couple
diff --git a/experiments/integrated/fluxnet/ozark_pft.jl b/experiments/integrated/fluxnet/ozark_pft.jl
@@ -302,13 +302,12 @@ diags = ClimaLand.default_diagnostics(
     start_date;
     output_writer = output_writer,
     output_vars,
-    average_period = :hourly,
+    average_period = :halfhourly,
 )
 
 ## How often we want to update the drivers
 ## defined in the simulatons file
-data_dt = Float64(FluxnetSimulations.get_data_dt(site_ID))
-updateat = Array(start_date:Second(data_dt):stop_date)
+updateat = Array(start_date:Second(dt):stop_date)
 simulation = LandSimulation(
     start_date,
     stop_date,
diff --git a/experiments/integrated/fluxnet/ozark_pmodel.jl b/experiments/integrated/fluxnet/ozark_pmodel.jl
@@ -246,12 +246,11 @@ diags = ClimaLand.default_diagnostics(
     start_date;
     output_writer = ClimaDiagnostics.Writers.DictWriter(),
     output_vars,
-    average_period = :hourly,
+    average_period = :halfhourly,
 );
 
 ## How often we want to update the drivers.
-data_dt = Second(FluxnetSimulations.get_data_dt(site_ID))
-updateat = Array(start_date:data_dt:stop_date)
+updateat = Array(start_date:Second(dt):stop_date)
 pmodel_cb = ClimaLand.make_PModel_callback(FT, start_date, dt, land.canopy)
 simulation = LandSimulation(
     start_date,
diff --git a/experiments/integrated/fluxnet/ozark_soilsnow.jl b/experiments/integrated/fluxnet/ozark_soilsnow.jl
@@ -186,6 +186,7 @@ simulation = LandSimulation(
     set_ic! = set_ic!,
     updateat,
     solver_kwargs = (; saveat = saveat),
+    diagnostics = nothing,
 )
 sol = solve!(simulation)
 
diff --git a/experiments/integrated/fluxnet/run_fluxnet.jl b/experiments/integrated/fluxnet/run_fluxnet.jl
@@ -270,13 +270,11 @@ diags = ClimaLand.default_diagnostics(
     start_date;
     output_writer = ClimaDiagnostics.Writers.DictWriter(),
     output_vars,
-    average_period = :hourly,
+    average_period = :halfhourly,
 );
 
 ## How often we want to update the drivers.
-data_dt = Second(FluxnetSimulations.get_data_dt(site_ID))
-updateat = Array(start_date:data_dt:stop_date)
-
+updateat = Array(start_date:Second(dt):stop_date)
 simulation = LandSimulation(
     start_date,
     stop_date,
diff --git a/experiments/standalone/Vegetation/no_vegetation.jl b/experiments/standalone/Vegetation/no_vegetation.jl
@@ -29,7 +29,7 @@ long = FT(-92.2000) # degree
 land_domain = Point(; z_sfc = FT(0.0), longlat = (long, lat))
 atmos_h = FT(32)
 site_ID = "US-MOz"
-start_date = DateTime(2010) + Hour(time_offset)
+start_date = DateTime(2010, 1, 2)
 N_days = 10
 stop_date = start_date + Day(N_days)
 dt = 225.0;
diff --git a/experiments/standalone/Vegetation/varying_lai.jl b/experiments/standalone/Vegetation/varying_lai.jl
@@ -28,7 +28,7 @@ long = FT(-92.2000) # degree
 land_domain = Point(; z_sfc = FT(0.0), longlat = (long, lat))
 atmos_h = FT(32)
 site_ID = "US-MOz"
-start_date = DateTime(2010) + Hour(time_offset)
+start_date = DateTime(2010, 1, 2)
 N_days = 364
 stop_date = start_date + Day(N_days)
 dt = 225.0;
diff --git a/experiments/standalone/Vegetation/varying_lai_with_stem.jl b/experiments/standalone/Vegetation/varying_lai_with_stem.jl
@@ -28,7 +28,7 @@ long = FT(-92.2000) # degree
 land_domain = Point(; z_sfc = FT(0.0), longlat = (long, lat))
 atmos_h = FT(32)
 site_ID = "US-MOz"
-start_date = DateTime(2010) + Hour(time_offset)
+start_date = DateTime(2010, 1, 2)
 N_days = 60
 stop_date = start_date + Day(N_days)
 dt = FT(225);
diff --git a/ext/fluxnet_simulations/forcing.jl b/ext/fluxnet_simulations/forcing.jl
@@ -23,7 +23,9 @@ precipitation that is in snow (if split_precip ==true), and (5)
 making the TimeVaryingInput objects.
 
 Note that the TimeVaryingInput objects can be used to interpolate in time,
-which is why we drop missing data.
+which is why we drop missing data. For the timestamp of the observation, we use
+the halfway point between the timestamp start and timestamp end:
+see https://fluxnet.org/data/fluxnet2015-dataset/fullset-data-product/ for details.
 
 This assumes that the first row of the CSV file is the list of of column names,
 and that these names are:
@@ -54,7 +56,17 @@ function FluxnetSimulations.prescribed_forcing_fluxnet(
     (data, columns) = readdlm(fluxnet_csv_path, ','; header = true)
 
     # Determine which column index corresponds to which varname
-    varnames = ("TA_F", "VPD_F", "PA_F", "P_F", "WS_F", "LW_IN_F", "SW_IN_F")
+    varnames = (
+        "TIMESTAMP_START",
+        "TIMESTAMP_END",
+        "TA_F",
+        "VPD_F",
+        "PA_F",
+        "P_F",
+        "WS_F",
+        "LW_IN_F",
+        "SW_IN_F",
+    )
     column_name_map = Dict(
         varname => findfirst(columns[:] .== varname) for varname in varnames
     )
@@ -68,7 +80,19 @@ function FluxnetSimulations.prescribed_forcing_fluxnet(
     # Convert the local timestamp to UTC
     # Since it was read in as Float64 type, convert to a string before
     # converting to a DateTime
-    local_datetime = DateTime.(string.(Int.(data[:, 1])), "yyyymmddHHMM")
+    local_datetime_start =
+        DateTime.(
+            string.(Int.(data[:, column_name_map["TIMESTAMP_START"]])),
+            "yyyymmddHHMM",
+        )
+    local_datetime_end =
+        DateTime.(
+            string.(Int.(data[:, column_name_map["TIMESTAMP_END"]])),
+            "yyyymmddHHMM",
+        )
+    local_datetime =
+        local_datetime_start .+
+        (local_datetime_end .- local_datetime_start) ./ 2
     UTC_datetime = local_datetime .+ Dates.Hour(hour_offset_from_UTC)
 
     # The TimeVaryingInput interface for columns expects the time in seconds
@@ -268,6 +292,10 @@ Fluxnet data at `site_ID`, given the offset in hours of local time
 from UTC. If `duration` is provided, it is used to determine the end date,
 otherwise the end date is the last date in the data. The `start_offset` is
 added to the start date, and must be non-negative.
+
+Please note that we use date halfway between the TIMESTAMP_START
+and TIMESTAMP_END date of the Fluxnet averaging window as the timestamp of
+the data for forcing.
 """
 function FluxnetSimulations.get_data_dates(
     site_ID,
@@ -277,7 +305,33 @@ function FluxnetSimulations.get_data_dates(
 )
     fluxnet_csv_path = ClimaLand.Artifacts.experiment_fluxnet_data_path(site_ID)
     (data, columns) = readdlm(fluxnet_csv_path, ','; header = true)
-    local_datetime = DateTime.(string.(Int.(data[:, 1])), "yyyymmddHHMM")
+    # Determine which column index corresponds to which varname
+    varnames = ("TIMESTAMP_START", "TIMESTAMP_END")
+    column_name_map = Dict(
+        varname => findfirst(columns[:] .== varname) for varname in varnames
+    )
+    # If any of these are missing, error, because we need all of them
+    nothing_id = findall(collect(values(column_name_map)) .== nothing)
+    if !isempty(nothing_id)
+        @error("$(labels[nothing_id]) is missing in the data, but required.")
+    end
+
+    # Convert the local timestamp to UTC
+    # Since it was read in as Float64 type, convert to a string before
+    # converting to a DateTime
+    local_datetime_start =
+        DateTime.(
+            string.(Int.(data[:, column_name_map["TIMESTAMP_START"]])),
+            "yyyymmddHHMM",
+        )
+    local_datetime_end =
+        DateTime.(
+            string.(Int.(data[:, column_name_map["TIMESTAMP_END"]])),
+            "yyyymmddHHMM",
+        )
+    local_datetime =
+        local_datetime_start .+
+        (local_datetime_end .- local_datetime_start) ./ 2
     UTC_datetime = local_datetime .+ Dates.Hour(hour_offset_from_UTC)
     earliest_date, latest_date = extrema(UTC_datetime)
     Dates.value(start_offset) < 0 && error("start_offset must be non-negative")
@@ -354,23 +408,38 @@ end
      get_comparison_data(
         site_ID,
         hour_offset_from_UTC;
-        val = -9999
+        val = -9999,
+        timestamp_end = true
 )
 
 Gets and returns the a NamedTuple with the comparison
 data for the Fluxnet `site_ID`, given its hour offset from
-UTC, and the value used to indicate missing data (`val`)
+UTC, and the value used to indicate missing data (`val`), and the
+averaging period.
+
+Fluxnet as has two timestamps associated with it: TIMESTAMP_START and TIMESTAMP_END,
+with the first referring to the start of the averaging period, and the second the end.
+ClimaLand diagnostics are averaged, accumulated, or otherwise reduced
+over a time period (e.g. hourly, daily, monthly). They are saved with the first date following
+that average period. For example, the hourly average from 11-noon is saved with a timestamp of
+noon. To make a true comparison to Fluxnet data, therefore, we must use halfhourly diagnostics in ClimaLAnd,
+and return the UTC time that corresponds to the end 
+of the averaging period in fluxnet; this is true by default. If `timestamp_end = false`,
+we return the point halfway between TIMESTAMP_START and TIMESTAMP_END
 """
 function FluxnetSimulations.get_comparison_data(
     site_ID::String,
     hour_offset_from_UTC::Int;
     val = -9999,
+    timestamp_end = true,
 )
     fluxnet_csv_path = ClimaLand.Artifacts.experiment_fluxnet_data_path(site_ID)
     (data, columns) = readdlm(fluxnet_csv_path, ','; header = true)
 
     # Determine which column index corresponds to which varname
     varnames = (
+        "TIMESTAMP_START",
+        "TIMESTAMP_END",
         "GPP_DT_VUT_REF",
         "LE_CORR",
         "H_CORR",
@@ -384,11 +453,26 @@ function FluxnetSimulations.get_comparison_data(
         varname => findfirst(columns[:] .== varname) for varname in varnames
     )
 
-    # Convert the local timestamp to UTC
-    local_datetime = DateTime.(string.(Int.(data[:, 1])), "yyyymmddHHMM")
-    UTC_datetime = local_datetime .+ Dates.Hour(hour_offset_from_UTC)
-    data_dt = Second(local_datetime[2] - local_datetime[1]).value # seconds
 
+    local_datetime_end =
+        DateTime.(
+            string.(Int.(data[:, column_name_map["TIMESTAMP_END"]])),
+            "yyyymmddHHMM",
+        )
+    data_dt = Second(local_datetime_end[2] - local_datetime_end[1]).value # seconds  
+    if timestamp_end
+        UTC_datetime = local_datetime_end .+ Dates.Hour(hour_offset_from_UTC)
+    else
+        local_datetime_start =
+            DateTime.(
+                string.(Int.(data[:, column_name_map["TIMESTAMP_START"]])),
+                "yyyymmddHHMM",
+            )
+        local_datetime =
+            local_datetime_start .+
+            (local_datetime_end .- local_datetime_start) ./ 2
+        UTC_datetime = local_datetime .+ Dates.Hour(hour_offset_from_UTC)
+    end
     gpp = FluxnetSimulations.get_comparison_data(
         data,
         "GPP_DT_VUT_REF",
@@ -397,6 +481,7 @@ function FluxnetSimulations.get_comparison_data(
         preprocess_func = (x) -> x * 1e-6, # converts from μmol/m^2/s to mol/m^2/s
         val,
     )
+
     lhf = FluxnetSimulations.get_comparison_data(
         data,
         "LE_CORR",
diff --git a/ext/land_sim_vis/plotting_utils.jl b/ext/land_sim_vis/plotting_utils.jl
@@ -189,7 +189,7 @@ function LandSimVis.check_conservation(
     end
     titles =
         ["Global mean energy per area", "Global mean water volume per area"]
-    name = ["energy", "water"]
+    quantity_names = ["energy", "water"]
     errors = [energy_error, water_volume_error]
     typical_value =
         [@sprintf("%1.2le", mean_energy), @sprintf("%1.2le", mean_water_volume)]
@@ -204,7 +204,7 @@ function LandSimVis.check_conservation(
         )
         CairoMakie.lines!(ax, times ./ 24 ./ 3600 ./ 365, errors[i])
         CairoMakie.save(
-            joinpath(savedir, "$(names[i])_$(plot_stem_name).pdf"),
+            joinpath(savedir, "$(quantity_names[i])_$(plot_stem_name).pdf"),
             fig_cycle,
         )
     end

Original file line number	Diff line number	Diff line change
`@@ -186,6 +186,7 @@ simulation = LandSimulation(`
`186`	`186`	`set_ic! = set_ic!,`
`187`	`187`	`updateat,`
`188`	`188`	`solver_kwargs = (; saveat = saveat),`
	`189`	`+ diagnostics = nothing,`
`189`	`190`	`)`
`190`	`191`	`sol = solve!(simulation)`
`191`	`192`