diff --git a/docs/src/datahandling.md b/docs/src/datahandling.md index fa6a3dee..80c5b94a 100644 --- a/docs/src/datahandling.md +++ b/docs/src/datahandling.md @@ -189,6 +189,8 @@ ClimaUtilities.DataHandling.available_times ClimaUtilities.DataHandling.available_dates ClimaUtilities.DataHandling.previous_time ClimaUtilities.DataHandling.next_time +ClimaUtilities.DataHandling.previous_date +ClimaUtilities.DataHandling.next_date ClimaUtilities.DataHandling.regridded_snapshot ClimaUtilities.DataHandling.regridded_snapshot! ClimaUtilities.DataHandling.dt diff --git a/ext/DataHandlingExt.jl b/ext/DataHandlingExt.jl index 82f71e09..b05b3d0b 100644 --- a/ext/DataHandlingExt.jl +++ b/ext/DataHandlingExt.jl @@ -511,6 +511,46 @@ function DataHandling.next_time(data_handler::DataHandler, date::Dates.DateTime) return data_handler.available_times[index] end +""" + DataHandling.previous_date(data_handler::DataHandler, time::Dates.TimeType) + +Return the date of the snapshot before the given `date`. +If `date` is one of the snapshots, return itself. + +If `date` is not in the `data_handler`, return an error. +""" +function DataHandling.previous_date( + data_handler::DataHandler, + date::Dates.TimeType, +) + if date in data_handler.available_dates + index = searchsortedfirst(data_handler.available_dates, date) + else + index = searchsortedfirst(data_handler.available_dates, date) - 1 + end + index < 1 && error("Date $date is before available dates") + return data_handler.available_dates[index] +end + +""" + DataHandling.next_date(data_handler::DataHandler, time::Dates.TimeType) + +Return the date of the snapshot after the given `time`. +If `date` is one of the snapshots, return itself. + +If `date` is not in the `data_handler`, return an error. +""" +function DataHandling.next_date(data_handler::DataHandler, date::Dates.TimeType) + if date in data_handler.available_dates + index = searchsortedfirst(data_handler.available_dates, date) + 1 + else + index = searchsortedfirst(data_handler.available_dates, date) + end + index > length(data_handler.available_dates) && + error("Date $date is after available dates") + return data_handler.available_dates[index] +end + """ regridded_snapshot(data_handler::DataHandler, date::Dates.DateTime) regridded_snapshot(data_handler::DataHandler, time::AbstractFloat) diff --git a/ext/TimeVaryingInputsExt.jl b/ext/TimeVaryingInputsExt.jl index d66a3d4b..292d551e 100644 --- a/ext/TimeVaryingInputsExt.jl +++ b/ext/TimeVaryingInputsExt.jl @@ -33,12 +33,13 @@ import ClimaUtilities.TimeVaryingInputs: extrapolation_bc import ClimaUtilities.DataHandling import ClimaUtilities.DataHandling: DataHandler, - previous_time, - next_time, regridded_snapshot!, available_times, available_dates, - date_to_time + date_to_time, + time_to_date, + previous_date, + next_date # Ideally, we should be able to split off the analytic part in a different # extension, but precompilation stops working when we do so @@ -110,9 +111,16 @@ end Check if the given `time` is in the range of definition for `itp`. """ function Base.in(time, itp::InterpolatingTimeVaryingInput23D) - return itp.range[1] <= time <= itp.range[2] + return itp.data_handler.available_dates[begin] <= + time <= + itp.data_handler.available_dates[end] end +function Base.in(time::Number, itp::InterpolatingTimeVaryingInput23D) + return Base.in(time_to_date(itp.data_handler, time), itp) +end + + function TimeVaryingInputs.TimeVaryingInput( data_handler; method = LinearInterpolation(), @@ -198,12 +206,13 @@ function TimeVaryingInputs.evaluate!( time in itp || error("TimeVaryingInput does not cover time $time") end if extrapolation_bc(itp.method) isa Flat - t_init, t_end = itp.range - if time >= t_end - regridded_snapshot!(dest, itp.data_handler, t_end) + date_init, date_end = itp.data_handler.available_dates[begin], + itp.data_handler.available_dates[end] + if time >= date_end + regridded_snapshot!(dest, itp.data_handler, date_end) else - time <= t_init - regridded_snapshot!(dest, itp.data_handler, t_init) + time <= date_init + regridded_snapshot!(dest, itp.data_handler, date_init) end else TimeVaryingInputs.evaluate!(dest, itp, time, itp.method) @@ -211,6 +220,23 @@ function TimeVaryingInputs.evaluate!( return nothing end +function TimeVaryingInputs.evaluate!( + dest, + itp::InterpolatingTimeVaryingInput23D, + time::Number, + args...; + kwargs..., +) + TimeVaryingInputs.evaluate!( + dest, + itp, + Dates.Millisecond(round(1_000 * time)) + itp.data_handler.start_date, + args..., + kwargs..., + ) + return nothing +end + function _time_range_dt_dt_e(itp::InterpolatingTimeVaryingInput23D) return _time_range_dt_dt_e(itp, extrapolation_bc(itp.method)) end @@ -220,7 +246,10 @@ function _time_range_dt_dt_e( extrapolation_bc::PeriodicCalendar{Nothing}, ) dt = DataHandling.dt(itp.data_handler) - return itp.range[begin], itp.range[end], dt, dt / 2 + return itp.data_handler.available_dates[begin], + itp.data_handler.available_dates[end], + Dates.Millisecond(round(1_000 * dt)), + Dates.Millisecond(round((1_000 * dt) / 2)) end function _time_range_dt_dt_e( @@ -240,11 +269,9 @@ function _time_range_dt_dt_e( date_to_time(itp.data_handler, date_end) # We have to add 1 Second because endofperiod(date_end, period) returns the very last # second before the next period - dt_e = - (endofperiod(date_end, period) + Dates.Second(1) - date_end) |> - period_to_seconds_float - dt = (date_init + period - date_end) |> period_to_seconds_float - return t_init, t_end, dt, dt_e + dt_e = (endofperiod(date_end, period) + Dates.Second(1) - date_end) + dt = (date_init + period - date_end) + return date_init, date_end, dt, dt_e end """ @@ -320,15 +347,14 @@ function TimeVaryingInputs.evaluate!( end return nothing end - - t0, t1 = - previous_time(itp.data_handler, time), next_time(itp.data_handler, time) + date0, date1 = + previous_date(itp.data_handler, time), next_date(itp.data_handler, time) # The closest regridded_snapshot could be either the previous or the next one - if (time - t0) <= (t1 - time) - regridded_snapshot!(dest, itp.data_handler, t0) + if (time - date0) <= (date1 - time) + regridded_snapshot!(dest, itp.data_handler, date0) else - regridded_snapshot!(dest, itp.data_handler, t1) + regridded_snapshot!(dest, itp.data_handler, date1) end return nothing end @@ -351,7 +377,7 @@ function TimeVaryingInputs.evaluate!( field_t0, field_t1 = itp.preallocated_regridded_fields[begin:(begin + 1)] if extrapolation_bc(itp.method) isa PeriodicCalendar - time, t_init, t_end, dt, _ = + time, date_init, date_end, dt, _ = _interpolation_times_periodic_calendar(time, itp) # We have to handle separately the edge case where the desired time is past t_end. @@ -361,10 +387,10 @@ function TimeVaryingInputs.evaluate!( # TODO: It would be nice to handle this edge case directly instead of copying the # code - if time > t_end - regridded_snapshot!(field_t0, itp.data_handler, t_end) - regridded_snapshot!(field_t1, itp.data_handler, t_init) - coeff = (time - t_end) / dt + if time > date_end + regridded_snapshot!(field_t0, itp.data_handler, date_end) + regridded_snapshot!(field_t1, itp.data_handler, date_init) + coeff = (time - date_end) / dt dest .= (1 - coeff) .* field_t0 .+ coeff .* field_t1 return nothing end @@ -372,15 +398,15 @@ function TimeVaryingInputs.evaluate!( # We have to consider the edge case where time is precisely the last available_time. # This is relevant also because it can be triggered by LinearPeriodFilling - if time in DataHandling.available_times(itp.data_handler) + if time in DataHandling.available_dates(itp.data_handler) regridded_snapshot!(dest, itp.data_handler, time) else - t0, t1 = previous_time(itp.data_handler, time), - next_time(itp.data_handler, time) - coeff = (time - t0) / (t1 - t0) + date0, date1 = previous_date(itp.data_handler, time), + next_date(itp.data_handler, time) + coeff = (time - date0) / (date1 - date0) - regridded_snapshot!(field_t0, itp.data_handler, t0) - regridded_snapshot!(field_t1, itp.data_handler, t1) + regridded_snapshot!(field_t0, itp.data_handler, date0) + regridded_snapshot!(field_t1, itp.data_handler, date1) dest .= (1 - coeff) .* field_t0 .+ coeff .* field_t1 end diff --git a/ext/time_varying_inputs_linearperiodfilling.jl b/ext/time_varying_inputs_linearperiodfilling.jl index 20f81f94..2fed7097 100644 --- a/ext/time_varying_inputs_linearperiodfilling.jl +++ b/ext/time_varying_inputs_linearperiodfilling.jl @@ -211,7 +211,7 @@ function TimeVaryingInputs.evaluate!( available_periods = unique_periods(available_dates, period) # E.g., Date(1987, 2, 1) - target_date = DataHandling.time_to_date(itp.data_handler, time) + target_date = time # If the target date falls within an available period, we just interpolate it with # LinearInterpolation() @@ -274,8 +274,6 @@ function TimeVaryingInputs.evaluate!( if in_interpolable_region date_pre = _move_date_to_period(target_date, period_left, period) date_post = _move_date_to_period(target_date, period_right, period) - time_pre = date_to_time(itp.data_handler, date_pre) - time_post = date_to_time(itp.data_handler, date_post) # Linear interpolation: y = y0 * (1 - coeff) + coeff * y1 # @@ -288,12 +286,12 @@ function TimeVaryingInputs.evaluate!( period_offset = period_right - period_left coeff = offset_periods_left / period_offset - TimeVaryingInputs.evaluate!(dest, itp, time_pre, LinearInterpolation()) + TimeVaryingInputs.evaluate!(dest, itp, date_pre, LinearInterpolation()) dest .*= (1 - coeff) TimeVaryingInputs.evaluate!( tmp_field1, itp, - time_post, + date_post, LinearInterpolation(), ) tmp_field1 .*= coeff @@ -359,17 +357,32 @@ function TimeVaryingInputs.evaluate!( error("We should not be here!") end - time_pre = date_to_time(itp.data_handler, date_pre) - time_post = date_to_time(itp.data_handler, date_post) - # y = y0 * (1 - coeff) + coff * y1 - TimeVaryingInputs.evaluate!(dest, itp, time_pre, method) - coeff = (time - time_pre) / (time_post - time_pre) + TimeVaryingInputs.evaluate!(dest, itp, date_pre, method) + coeff = (time - date_pre) / (date_post - date_pre) dest .*= (1 - coeff) - TimeVaryingInputs.evaluate!(tmp_field2, itp, time_post, method) + TimeVaryingInputs.evaluate!(tmp_field2, itp, date_post, method) tmp_field2 .*= coeff dest .+= tmp_field2 return nothing end return nothing end + +function TimeVaryingInputs.evaluate!( + dest, + itp::InterpolatingTimeVaryingInput23D, + time::Number, + method::LinearPeriodFillingInterpolation, + args...; + kwargs..., +) + TimeVaryingInputs.evaluate!( + dest, + itp, + Dates.Millisecond(round(1_000 * time)) + itp.data_handler.start_date, + args..., + kwargs..., + ) + return nothing +end diff --git a/src/DataHandling.jl b/src/DataHandling.jl index 0c56a326..842ff87a 100644 --- a/src/DataHandling.jl +++ b/src/DataHandling.jl @@ -51,6 +51,10 @@ function previous_time end function next_time end +function previous_date end + +function next_date end + function regridded_snapshot end function regridded_snapshot! end @@ -68,6 +72,8 @@ extension_fns = [ :available_dates, :previous_time, :next_time, + :previous_date, + :next_date, :regridded_snapshot, :regridded_snapshot!, :dt, @@ -80,6 +86,8 @@ extension_fns = [ :available_dates, :previous_time, :next_time, + :previous_date, + :next_date, :regridded_snapshot, :regridded_snapshot!, :dt, diff --git a/test/data_handling.jl b/test/data_handling.jl index bcbd9ea5..ad288d1d 100644 --- a/test/data_handling.jl +++ b/test/data_handling.jl @@ -361,6 +361,31 @@ end available_dates[10], ) == available_times[10] + # Previous date + @test_throws ErrorException DataHandling.previous_date( + data_handler, + available_dates[1] - Second(1), + ) + + @test DataHandling.previous_date( + data_handler, + available_dates[10] + Second(1), + ) == available_dates[10] + @test DataHandling.previous_date( + data_handler, + available_dates[1] + Second(1), + ) == available_dates[1] + + # Previous time with time, boundaries (return the node) + @test DataHandling.previous_date( + data_handler, + available_dates[1], + ) == available_dates[1] + @test DataHandling.previous_date( + data_handler, + available_dates[end], + ) == available_dates[end] + # Next time with time @test_throws ErrorException DataHandling.next_time( data_handler, @@ -391,6 +416,31 @@ end available_dates[10], ) == available_times[11] + # Next date + @test_throws ErrorException DataHandling.next_date( + data_handler, + available_dates[end] + Second(1), + ) + + @test DataHandling.next_date( + data_handler, + available_dates[10] + Second(1), + ) == available_dates[11] + @test DataHandling.next_date( + data_handler, + available_dates[1] + Second(1), + ) == available_dates[2] + + # On node + @test DataHandling.next_date( + data_handler, + available_dates[10], + ) == available_dates[11] + @test DataHandling.next_date( + data_handler, + available_dates[10], + ) == available_dates[11] + # Asking for a regridded_snapshot without specifying the time @test_throws ErrorException DataHandling.regridded_snapshot( data_handler, diff --git a/test/time_varying_inputs23.jl b/test/time_varying_inputs23.jl index aa99f9b7..c5878880 100644 --- a/test/time_varying_inputs23.jl +++ b/test/time_varying_inputs23.jl @@ -142,170 +142,199 @@ include("TestTools.jl") ), ), ) - @test 0.0 in input_nearest - @test !(1e23 in input_nearest) + @test !(82801.0 in input_nearest) + @test Dates.DateTime(2021, 1, 1, 0) in input_nearest + @test !(Dates.DateTime(2021, 1, 2, 0) in input_nearest) available_times = DataHandling.available_times(data_handler) + available_dates = DataHandling.available_dates(data_handler) dest = Fields.zeros(target_space) # Time outside of range - @test_throws ErrorException TimeVaryingInputs.evaluate!( - dest, - input_nearest, - FT(-40000), - ) + for t in (FT(-40000), Dates.DateTime(2020)) + @test_throws ErrorException TimeVaryingInputs.evaluate!( + dest, + input_nearest, + t, + ) + end # We are testing NearestNeighbor, so we can just have to check if the fields agree # Left nearest point target_time = available_times[10] + 1 - TimeVaryingInputs.evaluate!(dest, input_nearest, target_time) - - # We use isequal to handle NaNs - @test isequal( - Array(parent(dest)), - Array( - parent( - DataHandling.regridded_snapshot( - data_handler, - available_times[10], + target_date = available_dates[10] + Second(1) + for t in (target_time, target_date) + TimeVaryingInputs.evaluate!( + dest, + input_nearest, + target_time, + ) + + # We use isequal to handle NaNs + @test isequal( + Array(parent(dest)), + Array( + parent( + DataHandling.regridded_snapshot( + data_handler, + available_times[10], + ), ), ), - ), - ) + ) + end # Right nearest point target_time = available_times[9] - 1 - TimeVaryingInputs.evaluate!(dest, input_nearest, target_time) - - @test isequal( - Array(parent(dest)), - Array( - parent( - DataHandling.regridded_snapshot( - data_handler, - available_times[9], + target_date = available_dates[9] - Second(1) + for t in (target_time, target_date) + TimeVaryingInputs.evaluate!(dest, input_nearest, t) + + @test isequal( + Array(parent(dest)), + Array( + parent( + DataHandling.regridded_snapshot( + data_handler, + available_times[9], + ), ), ), - ), - ) + ) + end # On node target_time = available_times[11] - TimeVaryingInputs.evaluate!(dest, input_nearest, target_time) - - @test isequal( - Array(parent(dest)), - Array( - parent( - DataHandling.regridded_snapshot( - data_handler, - available_times[11], + target_date = available_dates[11] + for t in (target_time, target_date) + TimeVaryingInputs.evaluate!( + dest, + input_nearest, + target_time, + ) + + @test isequal( + Array(parent(dest)), + Array( + parent( + DataHandling.regridded_snapshot( + data_handler, + available_times[11], + ), ), ), - ), - ) + ) + end # Flat left target_time = available_times[begin] - 1 - TimeVaryingInputs.evaluate!( - dest, - input_nearest_flat, - target_time, - ) - - @test isequal( - Array(parent(dest)), - Array( - parent( - DataHandling.regridded_snapshot( - data_handler, - available_times[begin], + target_date = available_dates[begin] - Second(1) + for t in (target_time, target_date) + TimeVaryingInputs.evaluate!(dest, input_nearest_flat, t) + + @test isequal( + Array(parent(dest)), + Array( + parent( + DataHandling.regridded_snapshot( + data_handler, + available_times[begin], + ), ), ), - ), - ) + ) + end # Flat right target_time = available_times[end] + 1 - TimeVaryingInputs.evaluate!( - dest, - input_nearest_flat, - target_time, - ) - - @test isequal( - Array(parent(dest)), - Array( - parent( - DataHandling.regridded_snapshot( - data_handler, - available_times[end], + target_date = available_dates[end] + Second(1) + for t in (target_time, target_date) + TimeVaryingInputs.evaluate!(dest, input_nearest_flat, t) + + @test isequal( + Array(parent(dest)), + Array( + parent( + DataHandling.regridded_snapshot( + data_handler, + available_times[end], + ), ), ), - ), - ) + ) + end # Nearest periodic calendar dt = available_times[2] - available_times[1] target_time = available_times[end] + 0.1dt - TimeVaryingInputs.evaluate!( - dest, - input_nearest_periodic_calendar, - target_time, - ) + d_date = available_dates[2] - available_dates[1] + target_date = available_dates[end] + 0.1d_date + for t in (target_time, target_date) + TimeVaryingInputs.evaluate!( + dest, + input_nearest_periodic_calendar, + target_time, + ) - @test isequal( - Array(parent(dest)), - Array( - parent( - DataHandling.regridded_snapshot( - data_handler, - available_times[end], + @test isequal( + Array(parent(dest)), + Array( + parent( + DataHandling.regridded_snapshot( + data_handler, + available_times[end], + ), ), ), - ), - ) + ) + end # With date - TimeVaryingInputs.evaluate!( - dest, - input_nearest_periodic_calendar_date, - target_time, - ) + for t in (target_time, target_date) + TimeVaryingInputs.evaluate!( + dest, + input_nearest_periodic_calendar_date, + target_time, + ) - @test isequal( - Array(parent(dest)), - Array( - parent( - DataHandling.regridded_snapshot( - data_handler, - available_times[end], + @test isequal( + Array(parent(dest)), + Array( + parent( + DataHandling.regridded_snapshot( + data_handler, + available_times[end], + ), ), ), - ), - ) + ) + end dt = available_times[2] - available_times[1] target_time = available_times[end] + 0.6dt - TimeVaryingInputs.evaluate!( - dest, - input_nearest_periodic_calendar, - target_time, - ) + d_date = available_dates[2] - available_dates[1] + target_date = available_dates[end] + 0.6d_date + for t in (target_time, target_date) + TimeVaryingInputs.evaluate!( + dest, + input_nearest_periodic_calendar, + target_time, + ) - @test isequal( - Array(parent(dest)), - Array( - parent( - DataHandling.regridded_snapshot( - data_handler, - available_times[begin], + @test isequal( + Array(parent(dest)), + Array( + parent( + DataHandling.regridded_snapshot( + data_handler, + available_times[begin], + ), ), ), - ), - ) + ) + end # Now testing LinearInterpolation input_linear = TimeVaryingInputs.TimeVaryingInput(data_handler) @@ -316,6 +345,11 @@ include("TestTools.jl") input_linear, FT(-40000), ) + @test_throws ErrorException TimeVaryingInputs.evaluate!( + dest, + input_linear, + Dates.DateTime(2021, 1, 1) + Second(-40000), + ) left_value = DataHandling.regridded_snapshot( data_handler, @@ -330,7 +364,7 @@ include("TestTools.jl") left_time = available_times[10] right_time = available_times[11] - TimeVaryingInputs.evaluate!(dest, input_linear, target_time) + target_date = available_dates[10] + Second(30) expected = Fields.zeros(target_space) expected .= @@ -338,7 +372,10 @@ include("TestTools.jl") (target_time - left_time) / (right_time - left_time) .* (right_value .- left_value) - @test parent(dest) ≈ parent(expected) + for t in (target_time, target_date) + TimeVaryingInputs.evaluate!(dest, input_linear, t) + @test parent(dest) ≈ parent(expected) + end # LinearInterpolation with PeriodicCalendar time_delta = 0.1dt @@ -358,26 +395,32 @@ include("TestTools.jl") left_time = available_times[10] right_time = available_times[11] - TimeVaryingInputs.evaluate!( - dest, - input_linear_periodic_calendar, - target_time, - ) - expected = Fields.zeros(target_space) expected .= left_value .+ time_delta / dt .* (right_value .- left_value) - @test parent(dest) ≈ parent(expected) + target_date = available_dates[end] + 0.1d_date + for t in (target_time, target_date) + TimeVaryingInputs.evaluate!( + dest, + input_linear_periodic_calendar, + target_time, + ) + @test parent(dest) ≈ parent(expected) + end # With offset of one period - TimeVaryingInputs.evaluate!( - dest, - input_linear_periodic_calendar_date, - target_time + 86400, - ) + target_time += 86400 + target_date += Second(86400) + for t in (target_time, target_date) + TimeVaryingInputs.evaluate!( + dest, + input_linear_periodic_calendar_date, + t, + ) - @test parent(dest) ≈ parent(expected) + @test parent(dest) ≈ parent(expected) + end close(input_multiple_vars) close(input_nearest) diff --git a/test/utils.jl b/test/utils.jl index 670a64d2..26aab8d9 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -55,6 +55,18 @@ end # Wrapping when time equals t_init or t_end @test wrap_time(10, t_init, t_end) == 10 @test wrap_time(20, t_init, t_end) == 10 + + date_init = Dates.DateTime(2010) + Dates.Second(10) + date_end = Dates.DateTime(2010) + Dates.Second(20) + dt = Dates.Second(1) + + @test wrap_time(date_init + 5dt, date_init, date_end) == date_init + 5dt + @test wrap_time(date_end + 5dt, date_init, date_end) == date_init + 5dt + @test wrap_time(date_init - 5dt, date_init, date_end) == date_init + 5dt + + # Wrapping when time equals date_init or date_end + @test wrap_time(date_init, date_init, date_end) == date_init + @test wrap_time(date_end, date_init, date_end) == date_init end @testset "bounding_dates" begin