gempa_public/libcapsclient
7
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Files
44609d367f6a6de37cb92a1d3565b0c0cccbb36f
libcapsclient/libs/gempa/caps/mseedpacket.cpp
Enrico Ellguth 44609d367f Update to version 3.2
2026-03-18 15:18:05 +01:00

533 lines
14 KiB
C++
Raw Blame History

/***************************************************************************
* Copyright (C) 2009 by gempa GmbH *
* *
* All Rights Reserved. *
* *
* NOTICE: All information contained herein is, and remains *
* the property of gempa GmbH and its suppliers, if any. The intellectual *
* and technical concepts contained herein are proprietary to gempa GmbH *
* and its suppliers. *
* Dissemination of this information or reproduction of this material *
* is strictly forbidden unless prior written permission is obtained *
* from gempa GmbH. *
***************************************************************************/
#include <gempa/caps/endianess.h>
#include <gempa/caps/mseedpacket.h>
#include <gempa/caps/log.h>
#include <gempa/caps/utils.h>
#include <cmath>
#include <climits>
#include <cstdio>
#include <numeric>
#include <streambuf>
#include <iostream>
// Force local include
#include "./3rd-party/libmseed/libmseed.h"
#include "./3rd-party/libmseed/mseedformat.h"
#define MS_ISVALIDYEARDAY(Y,D) (Y >= 1900 && Y <= 2100 && D >= 1 && D <= 366)
#define LOG_SID(FUNC, format,...)\
do {\
char n[6];\
char s[6];\
char c[6];\
char l[6];\
ms_strncpclean(n, pMS2FSDH_NETWORK(head), 2);\
ms_strncpclean(s, pMS2FSDH_STATION(head), 5);\
ms_strncpclean(l, pMS2FSDH_LOCATION(head), 2);\
ms_strncpclean(c, pMS2FSDH_CHANNEL(head), 3);\
FUNC("[%s.%s.%s.%s] " format, n, s, l, c, ##__VA_ARGS__);\
} while(0)
namespace Gempa {
namespace CAPS {
namespace {
uint16_t ms2_blktlen(uint16_t blkttype, const char *blkt, int8_t swapflag) {
uint16_t blktlen = 0;
switch ( blkttype ) {
case 100: // Sample Rate
blktlen = 12;
break;
case 200: // Generic Event Detection
blktlen = 28;
break;
case 201: // Murdock Event Detection
blktlen = 36;
break;
case 300: // Step Calibration
blktlen = 32;
break;
case 310: // Sine Calibration
blktlen = 32;
break;
case 320: // Pseudo-random Calibration
blktlen = 28;
break;
case 390: // Generic Calibration
blktlen = 28;
break;
case 395: // Calibration Abort
blktlen = 16;
break;
case 400: // Beam
blktlen = 16;
break;
case 500: // Timing
blktlen = 8;
break;
case 1000: // Data Only SEED
blktlen = 8;
break;
case 1001: // Data Extension
blktlen = 8;
break;
case 2000: // Opaque Data
// First 2-byte field after the blockette header is the length
if ( blkt ) {
blktlen = *reinterpret_cast<const uint16_t*>(blkt + 4);
if ( swapflag ) {
ms_gswap2(&blktlen);
}
}
break;
}
return blktlen;
}
bool double2frac(uint16_t &numerator, uint16_t &denominator, double value) {
using namespace std;
// Check numeric limits
if ( value > numeric_limits<uint16_t>::max() ) {
numerator = numeric_limits<uint16_t>::max();
denominator = 1;
return false;
}
if ( value < numeric_limits<uint16_t>::min() ) {
numerator = numeric_limits<uint16_t>::min();
denominator = 1;
return false;
}
// Operate on positive numbers
int sign = 1;
if ( value < 0 ) {
sign = -1;
value = -value;
}
// Calculatate the largest possible power of 10 giving numeric integer
// limits and the current input number
static auto max_exp = floor(log10(numeric_limits<uint16_t>::max())) - 1.0;
auto exp = max_exp;
if ( value >= 10 ) {
exp -= floor(log10(value));
}
// Expand input number with power of 10
denominator = static_cast<int>(pow(10, exp));
numerator = static_cast<int>(round(value * denominator));
// Simplify the fraction by calculating the greatest common divisor
int gcd_ = gcd(numerator, denominator);
numerator = sign * numerator / gcd_;
denominator = denominator / gcd_;
return true;
}
}
MSEEDDataRecord::MSEEDDataRecord() {}
DataRecord *MSEEDDataRecord::clone() const {
return new MSEEDDataRecord(*this);
}
const char *MSEEDDataRecord::formatName() const {
return "MSEED";
}
bool MSEEDDataRecord::readMetaData(std::streambuf &buf, int size,
Header &header, Time &startTime, Time &endTime) {
char head[48];
if ( size <= 0 ) {
CAPS_WARNING("read metadata: invalid size of record: %d", size);
return false;
}
if ( (size < MINRECLEN) || (size > MAXRECLEN) ) {
CAPS_WARNING("read metadata: invalid MSEED record size: %d", size);
return false;
}
// Read first 40 byte
size_t read = buf.sgetn((char*)&head, sizeof(head));
if ( read < sizeof(head) ) {
CAPS_WARNING("read metadata: input buffer underflow: only %d/%d bytes read",
(int)read, (int)sizeof(head));
return false;
}
if ( MS2_ISVALIDHEADER(((char*)&head)) ) {
int headerswapflag = 0;
// Swap byte order?
if ( !MS_ISVALIDYEARDAY(*pMS2FSDH_YEAR(head), *pMS2FSDH_DAY(head)) ) {
headerswapflag = 1;
ms_gswap2(pMS2FSDH_YEAR(head));
ms_gswap2(pMS2FSDH_DAY(head));
ms_gswap2(pMS2FSDH_FSEC(head));
ms_gswap2(pMS2FSDH_NUMSAMPLES(head));
ms_gswap2(pMS2FSDH_SAMPLERATEFACT(head));
ms_gswap2(pMS2FSDH_SAMPLERATEMULT(head));
ms_gswap4(pMS2FSDH_TIMECORRECT(head));
ms_gswap2(pMS2FSDH_DATAOFFSET(head));
ms_gswap2(pMS2FSDH_BLOCKETTEOFFSET(head));
}
auto year = *pMS2FSDH_YEAR(head);
auto yday = *pMS2FSDH_DAY(head);
auto hours = *pMS2FSDH_HOUR(head);
auto minutes = *pMS2FSDH_MIN(head);
auto seconds = *pMS2FSDH_SEC(head);
auto fsec = *pMS2FSDH_FSEC(head) * 100;
header.dataType = DT_Unknown;
startTime = Time::FromYearDay(year, yday);
startTime += TimeSpan(hours * 3600 + minutes * 60 + seconds, fsec);
header.quality.ID = 0;
header.quality.str[0] = *pMS2FSDH_DATAQUALITY(head);
if ( *pMS2FSDH_TIMECORRECT(head) != 0 && !(*pMS2FSDH_ACTFLAGS(head) & 0x02) ) {
startTime += TimeSpan(0, *pMS2FSDH_TIMECORRECT(head) * 100);
}
// Parse blockettes
uint32_t blkt_offset = *pMS2FSDH_BLOCKETTEOFFSET(head);
uint32_t blkt_length;
uint16_t blkt_type;
uint16_t next_blkt;
if ( blkt_offset < sizeof(head) ) {
LOG_SID(CAPS_DEBUG, "%s", "read metadata: blockette "
"offset points into header");
return false;
}
uint32_t coffs = 0;
while ( (blkt_offset != 0) && ((int)blkt_offset < size) &&
(blkt_offset < static_cast<uint32_t>(size)) ) {
char bhead[6];
int seek_ofs = blkt_offset - sizeof(head) - coffs;
buf.pubseekoff(seek_ofs, std::ios_base::cur, std::ios_base::in);
coffs += seek_ofs;
if ( buf.sgetn(bhead, 6) != 6 ) {
LOG_SID(CAPS_DEBUG, "%s", "read metadata: "
"failed to read blockette header");
break;
}
coffs += 6;
memcpy(&blkt_type, bhead, 2);
memcpy(&next_blkt, bhead + 2, 2);
if ( headerswapflag ) {
ms_gswap2(&blkt_type);
ms_gswap2(&next_blkt);
}
blkt_length = ms2_blktlen(blkt_type, bhead, headerswapflag);
if ( blkt_length == 0 ) {
LOG_SID(CAPS_DEBUG, "read metadata: "
"unknown blockette length for type %d",
blkt_type);
break;
}
/* Make sure blockette is contained within the msrecord buffer */
if ( (int)(blkt_offset - 4 + blkt_length) > size ) {
LOG_SID(CAPS_DEBUG, "read metadata: blockette "
"%d extends beyond record size, truncated?",
blkt_type);
break;
}
if ( blkt_type == 1000 ) {
switch ( (int)bhead[4] ) {
case DE_ASCII:
header.dataType = DT_INT8;
break;
case DE_INT16:
header.dataType = DT_INT16;
break;
case DE_INT32:
case DE_STEIM1:
case DE_STEIM2:
case DE_CDSN:
case DE_DWWSSN:
case DE_SRO:
header.dataType = DT_INT32;
break;
case DE_FLOAT32:
header.dataType = DT_FLOAT;
break;
case DE_FLOAT64:
header.dataType = DT_DOUBLE;
break;
case DE_GEOSCOPE24:
case DE_GEOSCOPE163:
case DE_GEOSCOPE164:
header.dataType = DT_FLOAT;
break;
default:
break;
}
}
else if ( blkt_type == 1001 ) {
// Add usec correction
startTime += TimeSpan(0, *reinterpret_cast<int8_t*>(bhead + 5));
}
/* Check that the next blockette offset is beyond the current blockette */
if ( next_blkt && next_blkt < (blkt_offset + blkt_length - 4) ) {
LOG_SID(CAPS_DEBUG, "read metadata: offset to "
"next blockette (%d) is within current blockette "
"ending at byte %d",
blkt_type, (blkt_offset + blkt_length - 4));
break;
}
/* Check that the offset is within record length */
else if ( next_blkt && next_blkt > size ) {
LOG_SID(CAPS_DEBUG, "read metadata: offset to "
"next blockette (%d) from type %d is beyond record "
"length", next_blkt, blkt_type);
break;
}
else
blkt_offset = next_blkt;
}
endTime = startTime;
if ( *pMS2FSDH_SAMPLERATEFACT(head) > 0 ) {
header.samplingFrequencyNumerator = *pMS2FSDH_SAMPLERATEFACT(head);
header.samplingFrequencyDenominator = 1;
}
else {
header.samplingFrequencyNumerator = 1;
header.samplingFrequencyDenominator = -*pMS2FSDH_SAMPLERATEFACT(head);
}
if ( *pMS2FSDH_SAMPLERATEMULT(head) > 0 ) {
header.samplingFrequencyNumerator *= *pMS2FSDH_SAMPLERATEMULT(head);
}
else {
header.samplingFrequencyDenominator *= -*pMS2FSDH_SAMPLERATEMULT(head);
}
if ( header.samplingFrequencyNumerator > 0.0 && *pMS2FSDH_NUMSAMPLES(head) > 0 ) {
int64_t dt = static_cast<int64_t>(*pMS2FSDH_NUMSAMPLES(head)) * 1000000 * header.samplingFrequencyDenominator / header.samplingFrequencyNumerator;
endTime += TimeSpan(0, dt);
}
timeToTimestamp(header.samplingTime, startTime);
}
else if ( MS3_ISVALIDHEADER(((char*)&head)) ) {
startTime = Time::FromYearDay(
Endianess::Converter::ToLittleEndian(*pMS3FSDH_YEAR(head)),
Endianess::Converter::ToLittleEndian(*pMS3FSDH_DAY(head))
);
startTime += TimeSpan(
Endianess::Converter::ToLittleEndian(*pMS3FSDH_HOUR(head)) * 3600 +
Endianess::Converter::ToLittleEndian(*pMS3FSDH_MIN(head)) * 60 +
Endianess::Converter::ToLittleEndian(*pMS3FSDH_SEC(head)),
Endianess::Converter::ToLittleEndian(*pMS3FSDH_NSEC(head)) / 1000
);
auto nsamp = Endianess::Converter::ToLittleEndian(*pMS3FSDH_NUMSAMPLES(head));
auto fsamp = Endianess::Converter::ToLittleEndian(*pMS3FSDH_SAMPLERATE(head));
if ( fsamp < 0 ) {
fsamp = -1.0 / fsamp;
}
auto encoding = Endianess::Converter::ToLittleEndian(*pMS3FSDH_ENCODING(head));
header.dataType = DT_Unknown;
switch ( encoding ) {
case DE_ASCII:
header.dataType = DT_INT8;
break;
case DE_INT16:
header.dataType = DT_INT16;
break;
case DE_INT32:
case DE_STEIM1:
case DE_STEIM2:
case DE_CDSN:
case DE_DWWSSN:
case DE_SRO:
header.dataType = DT_INT32;
break;
case DE_FLOAT32:
header.dataType = DT_FLOAT;
break;
case DE_FLOAT64:
header.dataType = DT_DOUBLE;
break;
case DE_GEOSCOPE24:
case DE_GEOSCOPE163:
case DE_GEOSCOPE164:
header.dataType = DT_FLOAT;
break;
default:
break;
}
header.quality.ID = 0;
if ( !double2frac(header.samplingFrequencyNumerator, header.samplingFrequencyDenominator, fsamp) ) {
CAPS_WARNING("read metadata: invalid sampling rate: %f", fsamp);
return false;
}
endTime = startTime;
if ( header.samplingFrequencyNumerator > 0.0 && *pMS3FSDH_NUMSAMPLES(head) > 0 ) {
int64_t dt = static_cast<int64_t>(nsamp) * 1000000 * header.samplingFrequencyDenominator / header.samplingFrequencyNumerator;
endTime += TimeSpan(0, dt);
}
timeToTimestamp(header.samplingTime, startTime);
return false;
}
else {
CAPS_WARNING("read metadata: invalid MSEED header");
return false;
}
return true;
}
const DataRecord::Header *MSEEDDataRecord::header() const {
return &_header;
}
Time MSEEDDataRecord::startTime() const {
return _startTime;
}
Time MSEEDDataRecord::endTime() const {
return _endTime;
}
bool MSEEDDataRecord::canTrim() const {
return false;
}
bool MSEEDDataRecord::canMerge() const {
return false;
}
bool MSEEDDataRecord::trim(const Time &start,
const Time &end) const {
return false;
}
size_t MSEEDDataRecord::dataSize(bool /*withHeader*/) const {
return _data.size();
}
DataRecord::ReadStatus MSEEDDataRecord::get(std::streambuf &buf, int size,
const Time &start,
const Time &end,
int) {
//MSRecord *ms_rec = NULL;
if ( size <= 0 ) {
CAPS_WARNING("get: invalid size of record: %d", size);
return RS_Error;
}
_data.resize(size);
size_t read = buf.sgetn(&_data[0], _data.size());
if ( read != _data.size() ) {
CAPS_WARNING("get: input buffer underflow: only %d/%d bytes read",
(int)read, (int)_data.size());
return RS_Error;
}
arraybuf tmp(&_data[0], size);
readMetaData(tmp, size, _header, _startTime, _endTime);
if ( start.valid() || end.valid() ) {
// Out of scope?
if ( end.valid() && (end <= _startTime) )
return RS_AfterTimeWindow;
if ( start.valid() && (start >= _endTime) )
return RS_BeforeTimeWindow;
}
return RS_Complete;
}
bool MSEEDDataRecord::put(std::streambuf &buf, bool /*withHeader*/) const {
return buf.sputn(&_data[0], _data.size()) == (int)_data.size();
}
void MSEEDDataRecord::setData(const void *data, size_t size) {
_data.resize(size);
memcpy(_data.data(), data, size);
unpackHeader();
}
void MSEEDDataRecord::unpackHeader() {
arraybuf tmp(_data.data(), _data.size());
if ( !readMetaData(tmp, _data.size(), _header, _startTime, _endTime) ) {
CAPS_WARNING("invalid MSEED header");
}
}
}
}
Reference in New Issue View Git Blame Copy Permalink