Miniseed factory's endtime is not inclusive

To Do:

• adjust starttime/endtime by subtracting/adding a half-interval to the call to underlying clients/factories that retrieve data series, then use ObsPy class member functions to trim to the originally requested interval. If this can be generalized to work with all factories, that would be ideal, but at least do this for MiniSeedFactory and EdgeFactory.

Details:

If you use the MiniSeedFactory.get_timeseries() to pull 1-minute data from starttime​ to endtime​, you will get the correct length trace(s) returned (i.e., (endtime - starttime)/delta + 1​), but the last data point always seems to be a NaN. The MiniSeedFactory class wrapper to ObsPy's neic (miniseed) client, which in turn, simply opens a socket to a specified Edge host to request miniseed data using -s, -b, and -d parameters much like the java-based CWBQuery tool. If 500 microseconds (half-millisecond) is added to endtime, the expected top-of-the-minute value is returned correctly. In fact, if you request 1-second data, but the last requested second falls on a top-of-the-minute, everything occurs identically.

Here's an example using 1-minute data (assumes Edge running locally):

from geomagio.edge import MiniSeedFactory
from obspy.core import UTCDateTime

miniseed_factory = MiniSeedFactory(
    host = 'localhost',
    port = 2061,
    interval = 'minute',
    type = 'variation',
)

miniseed_stream = miniseed_factory.get_timeseries(
   starttime=UTCDateTime(2022, 10, 5, 0, 15),
   endtime=UTCDateTime(2022, 10, 5, 0, 20),
   observatory='BOU',
   channels=['U', 'V', 'W']
)

print(miniseed_stream[0].data)
[20768.2265625, 20768.23828125, 20768.26953125, 20768.45117188, 20768.4765625,            nan]

In contrast, if you use the EdgeFactory to pull data from starttime​ to endtime​, you will get the correct length trace(s) returned (i.e., (endtime - starttime)/delta + 1​), and the last data point is as expected. The EdgeFactory is a geomag-algorithms wrapper to ObsPy's earthworm (waveserver) client. Here's an example:

from geomagio.edge import EdgeFactory
from obspy.core import UTCDateTime

edge_factory = EdgeFactory(
    host = 'localhost',
    port = 2060,
    interval = 'minute',
    type = 'variation',
)

edge_stream = edge_factory.get_timeseries(
    starttime=UTCDateTime(2022, 10, 5, 0, 15),
    endtime=UTCDateTime(2022, 10, 5, 0, 20),
    observatory='BOU',
    channels=['H', 'E', 'Z']
)

print(edge_stream[0].data)
[20768.226, 20768.238, 20768.269, 20768.451, 20768.476, 20768.445]

Using Python's debugger, I was able to trace down a few additional potentially helpful facts:

• when requesting 1-minute data using the MiniSeedFactory (neic client) as above, 5 512-byte miniseed blocks are returned; in other words, it seems that there is 1 miniseed block per 1-minute data point.
• when requesting 1-minute data using the MiniSeedFactory (neic client) as above, but with endtime + 500 microseconds, 6 512-byte miniseed blocks are returned.
• when requesting 1-second data using the MiniSeedFactory (neic client), and an endtime that falls on the top-of-the-minute, the final miniseed block returned does NOT include the top-of-the-minute; however...
• when requesting 1-second data using the MiniSeedFactory (neic client), and 500 microseconds is added to an endtime that falls on the top-of-the-minute, the final miniseed block returne3d DOES include the top-of-the-minute.

This is a somewhat pressing problem because the Geomag program has at least one customer that regularly pulls 1-minute realtime data, and low-latency is a priority. The issue with the seconds data and top-of-the-minute gaps is not really all that bad, but when 1-minute data is delayed by 1-minute, this can be problematic. There is also reason to believe that our back-filling Python routines are not working quite as expected due to these top-of-the-minute gaps.

changed the description

changed the description

Disregard what I said about 1-second data above. I thought I had tested this, but when I tried a 1-second data request just now, it worked as expected. This is only a problem with 1-minute data.

Ugh, sorry for the clutter. It is a problem with 1-second data, but only if the last requested value falls on the top of a minute. Updating the issue description again.

changed the description

@ketchum : fyi, @erigler is seeing different results from the earthworm vs miniseed querymom endpoints. Can you comment on the expected behavior (whether endtime should be included in a response)?

The CWB EW returns 6 samples :

mage2::ketchum:~>cwbws -c GETSCNLRAW -b 2022/10/05-00:15 -e 2022/10/05-00:20 -dbg -s NTBOU--MVER0
17:45:43.183 cmd=GETSCNLRAW to localhost/2060 scnl=NTBOU  MVER0 start=2022/10/05-00:15 end=2022/10/05-00:20 fill=2147000000
cmd=GETSCNLRAW
GETSCNLRAW: GS BOU MVE NT R0 1664928900.000 1664929200.000
.
.
.
00:15:00.000 00:20:00.000 n=6 rt=0.016* [0]=-63048 [1]=-63056 [2]=-63115 [3]=-62953

      0:  -63048  -63056  -63115  -62953  -62849  -62996
End of execution

So I think if the sample matches exactly to the tenths of milliseconds, then the sample will be included using the get GETSCNLRAW command. It looks to me like the array is the right length, but that an NAN is in the last position. Perhaps the Python is copying the data and computing the number of samples incorrectly leaving out the last sample.

Hi Dave. Thanks for looking at this.

To clarify, it is not the waveserver client that is returning one-too-few samples, but rather obspy's NEIC client.

I'm not sure what the equivalent command for CWBQuery would be to match the output you copied above, but when I try the following in order to eliminate obspy as the culprit:

query -t sac -ascii -b 2022/10/05-00:15 -d 300 -s NTBOU--MVER0

I get a text file with the following:

5
2022,278 00:15:00.000
0.016666666666666666

I'm pretty sure the "5" is the number of returned values in the sac file (which I don't have software to read), since it increases to "6" when I request a duration of 360 seconds. I don't know why no actual data values are printed.

In any case, I've now debugged obspy down to the actual socket connection in its get_waveforms_nscl method. I am fairly certain it is the CWB that is dropping the last top-of-the-minute data point.

changed the description

I'm looping @bgeels and @awernle into this discussion in case they can add anything. I am also re-pinging @ketchum to ask him to read my updated description of the problem, and see if it rings any bells.

@erigler were you able to reproduce this issue using just obspy? I was unable to get the same problem to appear. I'm using obspy 1.2.2.

With:

from obspy import UTCDateTime
from obspy.clients.neic import Client

client = Client("cwbp2.cr.usgs.gov", port=22061)
t0=UTCDateTime(2022, 10, 5, 0, 15)
t1=UTCDateTime(2022, 10, 5, 0, 20)

st = client.get_waveforms("NT", "BOU", "R0", "MVH", t0, t1)

print(st[0].data)

I get:

[20768226 20768238 20768269 20768451 20768476 20768445]

Seconds and 10Hz data also appear to load properly for me. Same is true if I change the query server to mage1 directly.

@bgeels, you are requesting the "MVH" channel above, which is a copy of "UFU" stored the old pre-ObsRIO way (i.e., converted to scaled integers and uploaded via Raw Input Client). I don't have any trouble with those channels either.

But you got me to thinking...try yours with the following dates:

from obspy import UTCDateTime
from obspy.clients.neic import Client

client = Client("cwbp2.cr.usgs.gov", port=22061)
t0=UTCDateTime(2022, 10, 5, 4, 23, 0, 0)
t1=endtime=UTCDateTime(2022, 10, 5, 4, 28, 0, 0)
st = client.get_waveforms("NT", "BOU", "R0", "MVH", t0, t1)

print(st[0].data)

Do you get 5 instead of the expected 6 returned values? If you add 500 microseconds to t1, do you get the expected 6 returned values? Assuming you answer "yes", here's why:

• using the debugger, I was able to see the entire 512 byte miniseed packet returned by Edge (it is ObsPy that trims these to the desired length for subsequent processing);
• the final requested time stamp, in this case 2022-10-05T04:28:00Z, was the starttime for the next miniseed packet, which we don't receive unless we add 500 microseconds to t1, which somehow forces Edge to look for the next miniseed packet;

So, this is not an issue of ObsPy's waveserver vs. neic clients (or Geomag's thin wrappers to these clients), since they all use miniseed in the end anyway. However, it is more pronounced for the data we store the new ObsRIO way because, for some reason I don't fully understand, the miniseed packets used to store these data are aligned with times we frequently request, which happen to be tops-of-the-minute more often than not.

(on edit)

I will point out that if I use the waveserver client, it is not necessary to add the 500 microseconds. I get 6 returned values, as expected, when I use the dates above. Presumably, these are coming from the same miniseed packets.

I also see the behavior that you mentioned. [Edit]: I originally typed that up on 12/27 and now on 12/29 I'm getting 6 points! Try it again yourself.

I decided to play around with cwbquery on mage2 directly:

query -ascii -b 2022-10-05 04:23:00.000000 -d 300 -s NTBOU--MVER0

Returns 5 datapoints.

query -ascii -b 2022-10-05 04:23:00.000000 -d 329 -s NTBOU--MVER0

Also return 5 datapoints.

query -ascii -b 2022-10-05 04:23:00.000000 -d 330 -s NTBOU--MVER0

Returns 6 datapoints.

Similar behavior occurs with seconds data if you use a duration of 3.49 vs a duration of 3.5, and for 10Hz data with 0.349 vs 0.35, i.e.: You get another datapoint if the duration is increased by half of the sampling period or more. This seems to indicate that CWBquery is:

1. Rounding the duration to the nearest sampling period.
2. Is 'exclusive' rather than an 'inclusive' when it comes to the last point in the requested time period.

The second point above makes sense for a program that uses 'starttime & duration' since if you request 300 seconds of 1Hz data then you would probably expect to receive 300 points. Dave's example of cwbws above (and docs for obspy.clients.earthworm.client) seem to indicate that the earthworm & waveserver calls use the 'starttime & endtime' approach and (I presume) take an 'inclusive' approach with the last datapoint in the interval.

Perhaps as workaround we could automatically add an extra sampling_period/2 to endtime?

I am consistently getting the last datapoint included whenever I use ObsPy NEIC Client now, so maybe the above suggestion is unnecessary. I don't know why the behavior appears to be different when the data is queried via the socket connection.

Just a follow up to, and brief summary of, a recent Teams chat between @bgeels and @erigler:

• "top of the minute" issues are a red herring that arise sometimes because the data are uploaded to Edge in a particular way;
• using a Python debugger, we confirmed that Edge always returns full 512-byte MiniSEED packets that are subsequently trimmed to requested intervals using ObsPy functions;
• using a Python debugger, we confirmed that Edge returns no MiniSEED packet if a single time stamp is requested that corresponds to the starttime of a MiniSEED packet;
• using a Python debugger, I confirmed that Edge's earthworm/waveserver client returns no tracebuffs if a single time step is requested; this seems to have no relation to starttimes of MiniSEED packets;
• using a Python debugger, I confirmed that Edge's earthworm/waveserver client returns tracebuffs that are inclusive of non-zero-length requested intervals;

Finally, I (@erigler) will self-assign this issue and implement a version of the workaround suggested by @bgeels above in the Geomag-Algorithms project, since it seems like "inclusive" starttimes and endtimes are only desired by the Geomagnetism Program, and not by most users of Edge or ObsPy.

A few comments that might enlighten this. In most cases time series data does NOT align exactly on minute boundaries in seismology even if the data samples rate aligns, the group delay for the filtering would shift it back by the FIR filter delays, so the data is never aligned. Consider the much more common case (at least in seismology, you ask for 1 day of data. The requester does not know the alignment of the data, but if -s 2022,001:00:00:00.0000 -d 86400 the correct number of samples is returned - that is all of the data on that day. Normally you would want to get back all of the data for one 24 hour period or 86400*rate. If the data happened to be exactly aligned with the minute mark, would you want to get one more sample representing the first sample of the next day which in the CWB case is actually in a different file. In seismology that answer is no, never.

For the geomag case - if you ask for an exact period of time of 300 seconds from a minute boundary or any other boundary for 1 minute data, you would expect 5 samples, you asked for 5 minutes of data. For real time 1 minute data from geomag the data come in 5 sample packets. So if the last packet was for 01:01:00 the last sample is 01:05:00. If you query for 01:01:00 for 300 seconds you get back those 5 samples, if it worked such that you want to get 01:06:00 (inclusive interval) you would then have to wait for 5 more minutes for that data to show up. The software does the right thing and returns the last 5 mintute packet of data. If it worked the other way then asking for that same 300 seconds after the next real time packets shows up, gets 6 samples. Should the same query return different amounts of data?

The non-inclusive behavior from the CWB is NOT going to change - it is relied upon in seismology to insure the case above. Brendan's comment is correct in that you KNOW that the data aligns on the minute boundary so adding 1/2 of a sample interval will get you the next sample, but so would adding a whole sample interval. The actual adjustments the code makes to insure the interval is complete with random alighnment is actually 1/2 sample interval less 2 milliseconds.

assigned to @erigler

changed the description

closed with merge request !214 (merged)

mentioned in commit d1998c64

mentioned in commit c2c6c644