Add Events from External Data to Timetable

This example shows how to use events with timestamped data stored in timetables. An event consists of a timestamp (when something happened), often a description (what happened), and sometimes additional information about the event. The timestamp for an instantaneous event is a single point in time, either a datetime or duration value. The timestamp for an interval event consists of two points in time that define the start and end of the event.

You can use events to find or mark data of interest for plotting and analysis. This example shows how you can define events in data using information that is external to your data set. To find events that are within your timestamped data, see Find Events in Timestamped Data.

You can store events separately from the main data set. A timetable is a convenient way to store events because it can hold the event timestamps along with their descriptions and any other information. You can use the event timestamps as subscripts into the main data set to select data at, before, or after an event, or between two events. You can also store events in the same timetable as the data set they describe by using one or more additional timetable variables. This example shows how to use both representations. The representation that is more useful depends on the data analysis you plan to perform.

Define Data Set Using Variations in Length of Day

By definition a day is 86,400 seconds long, where the second has a precise definition in the International System of Units (SI). However, the length of a day actually varies due to several physical causes. It varies with the seasons by as much as 30 seconds over and 21 seconds under the SI definition because of the eccentricity of Earth's orbit and the tilt of its axis. Averaging these seasonal effects enables the definition of the mean solar day, which does not vary in length over a year.

Also, there is a very long-term slowing in the rotational speed of the Earth due to tidal interaction with the moon; a smaller, opposite, shorter-term component believed to be due to melting of continental ice sheets; very short-term cycles on the order of decades; and unpredictable fluctuations due to geological events and other causes. Because of those effects, the length of a mean solar day might increase or decrease. In recent decades, it has fluctuated up and down, but has mostly been 1–3 milliseconds longer than 86,400 seconds. That difference is known as the excess Length of Day, or excess LOD.

For this example, create a timetable that contains the excess LOD for every day from January 1, 1962, to the present. The International Earth Rotation and Reference Systems Service (IERS) collects and publishes this data. However, this data needs preprocessing before storing in a MATLAB timetable because the dates are modified Julian dates. To read the IERS data into a table, use the readtable function. Rename the two variables of interest to MJD and ExcessLOD.

file = "https://datacenter.iers.org/data/latestVersion/223_EOP_C04_14.62-NOW.IAU1980223.txt";
IERSdata = readtable(file,"NumHeaderLines",14);
IERSdata.Properties.VariableNames([4 8]) = ["MJD","ExcessLOD"];

To store the excess LOD values in a timetable, convert the modified Julian dates to datetime values. Use the datetime function with the "ConvertFrom","mjd" name-value argument. Convert the excess LOD values to duration values by using the seconds function. Then convert IERSdata from a table to a timetable using the table2timetable function.

IERSdata.Date = datetime(IERSdata.MJD,"ConvertFrom","mjd");
IERSdata.ExcessLOD = seconds(IERSdata.ExcessLOD);
IERSdata = table2timetable(IERSdata(:,["Date","ExcessLOD"]))

IERSdata=22017×1 timetable
       Date         ExcessLOD  
    ___________    ____________

    01-Jan-1962    0.001723 sec
    02-Jan-1962    0.001669 sec
    03-Jan-1962    0.001582 sec
    04-Jan-1962    0.001496 sec
    05-Jan-1962    0.001416 sec
    06-Jan-1962    0.001382 sec
    07-Jan-1962    0.001413 sec
    08-Jan-1962    0.001505 sec
    09-Jan-1962    0.001628 sec
    10-Jan-1962    0.001738 sec
    11-Jan-1962    0.001794 sec
    12-Jan-1962    0.001774 sec
    13-Jan-1962    0.001667 sec
    14-Jan-1962     0.00151 sec
    15-Jan-1962    0.001312 sec
    16-Jan-1962    0.001112 sec
      ⋮

Plot the excess LOD as a function of time. The excess LOD is currently decreasing on average but has remained positive except during very brief periods.

plot(IERSdata.Date,IERSdata.ExcessLOD,"b-");
ylabel("Excess LOD");

Calculate Cumulative Excess Length of Day

A clock that defines a day as 86,400 SI seconds is effectively running fast with respect to the Earth's rotation, gaining time with respect to the sun each day. A few extra milliseconds per day might seem unimportant, but the excess LOD since 1962 accumulates.

Add the cumulative excess LOD to IERSdata as another table variable. To calculate the cumulative sum of the excess LOD, use the cumsum function. Then use the stackedplot function to plot excess LOD and cumulative excess LOD together. This plot shows that the shift due to the accumulated excess LOD since the 1960s has been less than one minute.

IERSdata.CumulativeELOD = [0; cumsum(IERSdata.ExcessLOD(1:end-1))];
h = stackedplot(IERSdata,["ExcessLOD","CumulativeELOD"])

h = 
  StackedLineChart with properties:

         SourceTable: [22017×2 timetable]
    DisplayVariables: {'ExcessLOD'  'CumulativeELOD'}
               Color: [0 0.4470 0.7410]
           LineStyle: '-'
           LineWidth: 0.5000
              Marker: 'none'
          MarkerSize: 6

  Show all properties

Define Leap Seconds as Events in Separate TImetable

At midnight on January 1, 1972, the current version of the system of time known as Coordinated Universal Time (UTC) was enacted, which uses 86,400 SI seconds per day. On that date, UTC was defined to be roughly in sync with solar time at the Greenwich Meridian, or more precisely, with the time known as UT1. However, timekeeping based on exactly 86,400 SI seconds per day would drift away from our physical experience of solar time because of accumulated excess LOD. So when needed, an extra leap second adjustment is inserted to keep UTC approximately in sync with solar time. Without these leap second adjustments after 1972, accumulated excess LOD would have caused UTC to drift away from UT1 over the years. Within decades, the difference would have grown to tens of seconds.

Show the difference since 1972. First, select the post-1971 excess LOD data by subscripting into IERSdata with a time range starting on January 1, 1972. To create that time range, use the timerange function.

IERSdata1972 = IERSdata(timerange("1972-01-01",Inf),"ExcessLOD")

IERSdata1972=18365×1 timetable
       Date         ExcessLOD  
    ___________    ____________

    01-Jan-1972    0.002539 sec
    02-Jan-1972    0.002708 sec
    03-Jan-1972    0.002897 sec
    04-Jan-1972    0.003065 sec
    05-Jan-1972    0.003177 sec
    06-Jan-1972     0.00322 sec
    07-Jan-1972    0.003201 sec
    08-Jan-1972    0.003137 sec
    09-Jan-1972    0.003033 sec
    10-Jan-1972    0.002898 sec
    11-Jan-1972    0.002772 sec
    12-Jan-1972    0.002672 sec
    13-Jan-1972    0.002621 sec
    14-Jan-1972    0.002642 sec
    15-Jan-1972     0.00274 sec
    16-Jan-1972    0.002937 sec
      ⋮

Add another variable with the unadjusted difference since 1972. Start with the cumulative excess LOD on January 1, 1972, which was 0.0454859 second. For every date that follows, calculate the unadjusted difference. By the beginning of 2017, the unadjusted difference would have grown to about 26 seconds. Display IERSdata1972 using a time range that includes the start of 2017.

DiffUT1_1972 = seconds(0.0454859);
IERSdata1972.UnadjustedDiff = DiffUT1_1972 + [0; cumsum(IERSdata1972.ExcessLOD(1:end-1))];
IERSdata1972(timerange("2016-12-29","2017-01-04"),:)

ans=6×2 timetable
       Date          ExcessLOD      UnadjustedDiff
    ___________    _____________    ______________

    29-Dec-2016    0.0008055 sec      26.402 sec  
    30-Dec-2016    0.0008525 sec      26.402 sec  
    31-Dec-2016    0.0009173 sec      26.403 sec  
    01-Jan-2017     0.001016 sec      26.404 sec  
    02-Jan-2017    0.0011845 sec      26.405 sec  
    03-Jan-2017    0.0013554 sec      26.406 sec  

This drift is what leap seconds are designed to mitigate. You can think of each leap second that is inserted into the UTC timeline as an event. These events are not part of the LOD data. However, the leapseconds function lists each of the leap second events that has occurred since 1972. (The leap second events listed by leapseconds come from another data set provided by the IERS.) The timetable returned by leapseconds contains a timestamp, a description of each event (+ for leap second insertions, - for removals), and the cumulative number of leap seconds up to and including the event. To date, all leap second events have been insertions.

lsEvents = leapseconds()

lsEvents=27×2 timetable
       Date        Type    CumulativeAdjustment
    ___________    ____    ____________________

    30-Jun-1972     +              1 sec       
    31-Dec-1972     +              2 sec       
    31-Dec-1973     +              3 sec       
    31-Dec-1974     +              4 sec       
    31-Dec-1975     +              5 sec       
    31-Dec-1976     +              6 sec       
    31-Dec-1977     +              7 sec       
    31-Dec-1978     +              8 sec       
    31-Dec-1979     +              9 sec       
    30-Jun-1981     +             10 sec       
    30-Jun-1982     +             11 sec       
    30-Jun-1983     +             12 sec       
    30-Jun-1985     +             13 sec       
    31-Dec-1987     +             14 sec       
    31-Dec-1989     +             15 sec       
    31-Dec-1990     +             16 sec       
      ⋮

At this point, the main LOD data is in the IERSdata1972 timetable. The leap second events are in the lsEvents timetable.

Use Event Timestamps to Select Data

Add the unadjusted difference, or cumulative excess LOD, at each date in lsEvents. To find the unadjusted differences on these dates, index into the UnadjustedDiff variable of IERSdata1972 and select the unadjusted difference at the timestamps from lsEvents. Then append those differences to each event in lsEvents as additional information about the event.

lsEvents.UnadjustedDiff = IERSdata1972.UnadjustedDiff(lsEvents.Date)

lsEvents=27×3 timetable
       Date        Type    CumulativeAdjustment    UnadjustedDiff
    ___________    ____    ____________________    ______________

    30-Jun-1972     +              1 sec            0.63494 sec  
    31-Dec-1972     +              2 sec             1.1864 sec  
    31-Dec-1973     +              3 sec             2.2976 sec  
    31-Dec-1974     +              4 sec              3.289 sec  
    31-Dec-1975     +              5 sec             4.2718 sec  
    31-Dec-1976     +              6 sec             5.3334 sec  
    31-Dec-1977     +              7 sec             6.3469 sec  
    31-Dec-1978     +              8 sec              7.398 sec  
    31-Dec-1979     +              9 sec             8.3526 sec  
    30-Jun-1981     +             10 sec             9.6281 sec  
    30-Jun-1982     +             11 sec             10.389 sec  
    30-Jun-1983     +             12 sec             11.249 sec  
    30-Jun-1985     +             13 sec             12.451 sec  
    31-Dec-1987     +             14 sec             13.634 sec  
    31-Dec-1989     +             15 sec             14.669 sec  
    31-Dec-1990     +             16 sec             15.379 sec  
      ⋮

Plot Events Against Data

The lsEvents timetable shows that at each event, the IERS inserted a leap second whenever the excess LOD accumulated by roughly one additional second. To confirm this observation visually, plot the instantaneous leap second events as points overlaid on the LOD data over time.

plot(IERSdata1972.Date,IERSdata1972.UnadjustedDiff);
hold on
plot(lsEvents.Date,IERSdata1972.UnadjustedDiff(lsEvents.Date),"r+");
hold off

Represent Events in Additional Timetable Variable

Each leap second is an instantaneous event that occurs at a specific time. Another way to represent them is as a state variable appended to the original LOD data. A state variable describes the "state" of the process being measured at each time point in the data rather than tagging a specific instant. One possible state variable for this leap second data is an indicator variable that records an adjustment on each event date, with missing values on all the other dates. Assign the leap second event types to a state variable in IERSdata1972, using timestamps from lsEvents to assign event types to the correct rows. This approach works because all timestamps from lsEvents are also timestamps in IERSdata1972. The assignment automatically fills the other elements of the new Adjustment variable with missing values.

IERSdata1972.Adjustment(lsEvents.Date) = lsEvents.Type;
IERSdata1972(timerange("2016-12-29","2017-01-04"),:)

ans=6×3 timetable
       Date          ExcessLOD      UnadjustedDiff    Adjustment 
    ___________    _____________    ______________    ___________

    29-Dec-2016    0.0008055 sec      26.402 sec      <undefined>
    30-Dec-2016    0.0008525 sec      26.402 sec      <undefined>
    31-Dec-2016    0.0009173 sec      26.403 sec      +          
    01-Jan-2017     0.001016 sec      26.404 sec      <undefined>
    02-Jan-2017    0.0011845 sec      26.405 sec      <undefined>
    03-Jan-2017    0.0013554 sec      26.406 sec      <undefined>

Another possibility, more useful in this example, is to add a state variable that indicates the cumulative sum of leap seconds at any given date in the LOD data. Leap seconds occur at the end of a day, so first assign a value of one second to the day after the date for each event.

IERSdata1972.Adjustment = []; % remove the previous state variable
IERSdata1972.Adjustment(lsEvents.Date+caldays(1)) = seconds(1);
IERSdata1972(timerange("2016-12-29","2017-01-04"),:)

ans=6×3 timetable
       Date          ExcessLOD      UnadjustedDiff    Adjustment
    ___________    _____________    ______________    __________

    29-Dec-2016    0.0008055 sec      26.402 sec        0 sec   
    30-Dec-2016    0.0008525 sec      26.402 sec        0 sec   
    31-Dec-2016    0.0009173 sec      26.403 sec        0 sec   
    01-Jan-2017     0.001016 sec      26.404 sec        1 sec   
    02-Jan-2017    0.0011845 sec      26.405 sec        0 sec   
    03-Jan-2017    0.0013554 sec      26.406 sec        0 sec   

Then use a cumulative sum to fill in the Adjustment state variable at all other times in the data. Instead of containing values of + and <undefined>, the state variable Adjustment now lists the number of leap seconds added as of the date listed in each row of IERSdata1972.

IERSdata1972.Adjustment = cumsum(IERSdata1972.Adjustment);
IERSdata1972(timerange("2016-12-29","2017-01-04"),:)

ans=6×3 timetable
       Date          ExcessLOD      UnadjustedDiff    Adjustment
    ___________    _____________    ______________    __________

    29-Dec-2016    0.0008055 sec      26.402 sec        26 sec  
    30-Dec-2016    0.0008525 sec      26.402 sec        26 sec  
    31-Dec-2016    0.0009173 sec      26.403 sec        26 sec  
    01-Jan-2017     0.001016 sec      26.404 sec        27 sec  
    02-Jan-2017    0.0011845 sec      26.405 sec        27 sec  
    03-Jan-2017    0.0013554 sec      26.406 sec        27 sec  

Although not used here, the retime function is often a useful tool for converting instantaneous events to a state variable.

While the new state variable represents the instantaneous events, it is just like all the other variables in IERSdata1972, with a value that is defined at each time. You can use it to compute the actual differences between UTC and UT1, given all the leap second adjustments, by subtracting it from the unadjusted difference. The convention is to compute the actual difference with the opposite sign as UT1 – UTC and denote it as DUT1.

IERSdata1972.DUT1 = IERSdata1972.UnadjustedDiff - IERSdata1972.Adjustment;
IERSdata1972(timerange("2016-12-29","2017-01-04"),:)

ans=6×4 timetable
       Date          ExcessLOD      UnadjustedDiff    Adjustment        DUT1    
    ___________    _____________    ______________    __________    ____________

    29-Dec-2016    0.0008055 sec      26.402 sec        26 sec       0.40159 sec
    30-Dec-2016    0.0008525 sec      26.402 sec        26 sec       0.40239 sec
    31-Dec-2016    0.0009173 sec      26.403 sec        26 sec       0.40324 sec
    01-Jan-2017     0.001016 sec      26.404 sec        27 sec      -0.59584 sec
    02-Jan-2017    0.0011845 sec      26.405 sec        27 sec      -0.59482 sec
    03-Jan-2017    0.0013554 sec      26.406 sec        27 sec      -0.59364 sec

Plot State Variable Against Data

The IERS makes leap second adjustments to keep UTC roughly in sync with UT1. To show how this adjustment works, plot the Adjustment state variable as a piecewise-constant step function over time. The step function approximates the accumulated excess LOD, so subtracting it keeps DUT1 small. For visual confirmation, plot DUT1 along the same time axis by using the stackedplot function. Plot UnadjustedDiff and Adjustment together along one y-axis, with a legend, and plot DUT1 along a second y-axis.

h = stackedplot(IERSdata1972,{["UnadjustedDiff","Adjustment"],"DUT1"})

h = 
  StackedLineChart with properties:

         SourceTable: [18365×4 timetable]
    DisplayVariables: {{1×2 cell}  'DUT1'}
               Color: [0 0.4470 0.7410]
           LineStyle: '-'
           LineWidth: 0.5000
              Marker: 'none'
          MarkerSize: 6

  Show all properties

Add the value of DUT1 before and after each leap second to lsEvents. To get these values, use the instantaneous form of the leap second events. Index into the DUT1 variable of IERSdata1972 using the timestamps from lsEvents. Get the value of DUT1 on the date each leap second was added and on the date of the day after. In general, DUT1 has a small negative value the day after a leap second is added.

lsEvents.UnadjustedDiff = []; % remove the previous additional event information
lsEvents.DUT1_before = IERSdata1972.DUT1(lsEvents.Date);
lsEvents.DUT1_after = IERSdata1972.DUT1(lsEvents.Date+caldays(1))

lsEvents=27×4 timetable
       Date        Type    CumulativeAdjustment    DUT1_before     DUT1_after 
    ___________    ____    ____________________    ___________    ____________

    30-Jun-1972     +              1 sec           0.63494 sec    -0.36236 sec
    31-Dec-1972     +              2 sec           0.18641 sec     -0.8107 sec
    31-Dec-1973     +              3 sec           0.29763 sec    -0.69952 sec
    31-Dec-1974     +              4 sec           0.28898 sec    -0.70817 sec
    31-Dec-1975     +              5 sec            0.2718 sec    -0.72578 sec
    31-Dec-1976     +              6 sec           0.33338 sec    -0.66411 sec
    31-Dec-1977     +              7 sec           0.34693 sec    -0.64986 sec
    31-Dec-1978     +              8 sec           0.39802 sec    -0.59892 sec
    31-Dec-1979     +              9 sec           0.35261 sec    -0.64508 sec
    30-Jun-1981     +             10 sec            0.6281 sec    -0.37053 sec
    30-Jun-1982     +             11 sec           0.38922 sec    -0.60911 sec
    30-Jun-1983     +             12 sec           0.24852 sec    -0.75008 sec
    30-Jun-1985     +             13 sec           0.45106 sec    -0.54797 sec
    31-Dec-1987     +             14 sec           0.63443 sec    -0.36417 sec
    31-Dec-1989     +             15 sec           0.66947 sec    -0.32879 sec
    31-Dec-1990     +             16 sec           0.37936 sec    -0.61885 sec
      ⋮

Choose Events or State Variables

The two representations, a separate list of instantaneous events and a state variable defined at all times, are conceptually different. In some cases, there might not be a useful definition of a "state" that corresponds to the periods between instantaneous events. But when there is, the two representations are equivalent and useful in similar ways. For example, to select all the data after the 22nd leap second but before the 23rd, you can use event times and timerange.

from22to23 = timerange(lsEvents.Date(22),lsEvents.Date(23),"openleft");
data22to23 = IERSdata1972(from22to23,:)

data22to23=2557×4 timetable
       Date          ExcessLOD      UnadjustedDiff    Adjustment        DUT1    
    ___________    _____________    ______________    __________    ____________

    01-Jan-1999    0.0009738 sec      21.283 sec        22 sec      -0.71702 sec
    02-Jan-1999     0.000888 sec      21.284 sec        22 sec      -0.71604 sec
    03-Jan-1999    0.0008605 sec      21.285 sec        22 sec      -0.71516 sec
    04-Jan-1999    0.0008798 sec      21.286 sec        22 sec      -0.71429 sec
    05-Jan-1999    0.0008935 sec      21.287 sec        22 sec      -0.71342 sec
    06-Jan-1999    0.0009693 sec      21.287 sec        22 sec      -0.71252 sec
    07-Jan-1999    0.0010658 sec      21.288 sec        22 sec      -0.71155 sec
    08-Jan-1999    0.0010857 sec       21.29 sec        22 sec      -0.71049 sec
    09-Jan-1999    0.0010715 sec      21.291 sec        22 sec       -0.7094 sec
    10-Jan-1999    0.0010519 sec      21.292 sec        22 sec      -0.70833 sec
    11-Jan-1999    0.0010021 sec      21.293 sec        22 sec      -0.70728 sec
    12-Jan-1999    0.0008986 sec      21.294 sec        22 sec      -0.70628 sec
    13-Jan-1999    0.0007891 sec      21.295 sec        22 sec      -0.70538 sec
    14-Jan-1999    0.0007236 sec      21.295 sec        22 sec      -0.70459 sec
    15-Jan-1999    0.0006842 sec      21.296 sec        22 sec      -0.70386 sec
    16-Jan-1999    0.0006091 sec      21.297 sec        22 sec      -0.70318 sec
      ⋮

As an alternative, you can create a logical subscript from the state variable that selects the same data.

from22to23 = (IERSdata1972.Adjustment == seconds(22));
data22to23 = IERSdata1972(from22to23,:)

data22to23=2557×4 timetable
       Date          ExcessLOD      UnadjustedDiff    Adjustment        DUT1    
    ___________    _____________    ______________    __________    ____________

    01-Jan-1999    0.0009738 sec      21.283 sec        22 sec      -0.71702 sec
    02-Jan-1999     0.000888 sec      21.284 sec        22 sec      -0.71604 sec
    03-Jan-1999    0.0008605 sec      21.285 sec        22 sec      -0.71516 sec
    04-Jan-1999    0.0008798 sec      21.286 sec        22 sec      -0.71429 sec
    05-Jan-1999    0.0008935 sec      21.287 sec        22 sec      -0.71342 sec
    06-Jan-1999    0.0009693 sec      21.287 sec        22 sec      -0.71252 sec
    07-Jan-1999    0.0010658 sec      21.288 sec        22 sec      -0.71155 sec
    08-Jan-1999    0.0010857 sec       21.29 sec        22 sec      -0.71049 sec
    09-Jan-1999    0.0010715 sec      21.291 sec        22 sec       -0.7094 sec
    10-Jan-1999    0.0010519 sec      21.292 sec        22 sec      -0.70833 sec
    11-Jan-1999    0.0010021 sec      21.293 sec        22 sec      -0.70728 sec
    12-Jan-1999    0.0008986 sec      21.294 sec        22 sec      -0.70628 sec
    13-Jan-1999    0.0007891 sec      21.295 sec        22 sec      -0.70538 sec
    14-Jan-1999    0.0007236 sec      21.295 sec        22 sec      -0.70459 sec
    15-Jan-1999    0.0006842 sec      21.296 sec        22 sec      -0.70386 sec
    16-Jan-1999    0.0006091 sec      21.297 sec        22 sec      -0.70318 sec
      ⋮

Both representations have their uses. For example, while each event can be plotted as a point, the state variable is the more convenient form for highlighting regions between events in a plot. To highlight the region between the 22nd and 23rd leap second, use the isbetween function with the Adjustment state variable from IERSdata1972.

plot(IERSdata1972.Date,IERSdata1972.UnadjustedDiff);
hold on
from22to23 = isbetween(IERSdata1972.Adjustment,seconds(22),seconds(23));
plot(IERSdata1972.Date(from22to23),IERSdata1972.UnadjustedDiff(from22to23),"r-","LineWidth",4);
hold off

You can switch between the two representations. In this case, to get the leap second event dates from the adjustment state variable, find the locations where the adjustment changes, and subtract one day. The eventTimes output represents the dates on which leap seconds were added, which are instantaneous events.

eventTimes = IERSdata1972.Date(diff(IERSdata1972.Adjustment) ~= 0) - caldays(1)

eventTimes = 27×1 datetime
   29-Jun-1972
   30-Dec-1972
   30-Dec-1973
   30-Dec-1974
   30-Dec-1975
   30-Dec-1976
   30-Dec-1977
   30-Dec-1978
   30-Dec-1979
   29-Jun-1981
   29-Jun-1982
   29-Jun-1983
   29-Jun-1985
   30-Dec-1987
   30-Dec-1989
   30-Dec-1990
   29-Jun-1992
   29-Jun-1993
   29-Jun-1994
   30-Dec-1995
   29-Jun-1997
   30-Dec-1998
   30-Dec-2005
   30-Dec-2008
   29-Jun-2012
   29-Jun-2015
   30-Dec-2016

Whether you use instantaneous events, state variables, or switch between them depends on which representation is more convenient and useful for the data analysis that you plan to perform. Both representations are useful ways to add information about events to your timestamped data in a timetable.