teach reporter to do ephemeris calculations (still using the WRONG af0!) and feed them into influxdb

pull/125/head
bert hubert 2020-07-03 23:23:24 +02:00
parent b62e01197d
commit 2eb900c433
2 changed files with 364 additions and 32 deletions

View File

@ -79,7 +79,7 @@ decrypt: decrypt.o bits.o ext/fmt-6.1.2/src/format.o
navparse: navparse.o ext/fmt-6.1.2/src/format.o $(H2OPP) $(SIMPLESOCKETS) minicurl.o ubx.o bits.o navmon.pb.o gps.o ephemeris.o beidou.o glonass.o $(patsubst %.cc,%.o,$(wildcard ext/sgp4/libsgp4/*.cc)) tle.o navmon.o coverage.o osen.o trkmeas.o influxpush.o ${EXTRADEP} githash.o sbas.o rtcm.o
$(CXX) -std=gnu++17 $^ -o $@ -pthread -L/usr/local/lib -L/usr/local/opt/openssl/lib/ -lh2o-evloop -lssl -lcrypto -lz -lcurl -lprotobuf $(WSLAY)
reporter: reporter.o ext/fmt-6.1.2/src/format.o $(SIMPLESOCKETS) minicurl.o ubx.o bits.o navmon.pb.o gps.o ephemeris.o beidou.o glonass.o $(patsubst %.cc,%.o,$(wildcard ext/sgp4/libsgp4/*.cc)) tle.o navmon.o coverage.o osen.o githash.o
reporter: reporter.o ext/fmt-6.1.2/src/format.o $(SIMPLESOCKETS) minicurl.o ubx.o bits.o navmon.pb.o gps.o ephemeris.o beidou.o glonass.o $(patsubst %.cc,%.o,$(wildcard ext/sgp4/libsgp4/*.cc)) tle.o navmon.o coverage.o osen.o githash.o influxpush.o
$(CXX) -std=gnu++17 $^ -o $@ -pthread -L/usr/local/lib -lprotobuf -lcurl
sp3feed: sp3feed.o ext/fmt-6.1.2/src/format.o $(SIMPLESOCKETS) minicurl.o ubx.o bits.o navmon.pb.o gps.o ephemeris.o beidou.o glonass.o $(patsubst %.cc,%.o,$(wildcard ext/sgp4/libsgp4/*.cc)) tle.o navmon.o coverage.o osen.o influxpush.o githash.o sp3.o

View File

@ -4,9 +4,13 @@
#include "navmon.hh"
#include "fmt/format.h"
#include "fmt/printf.h"
#include "galileo.hh"
#include "gps.hh"
#include "CLI/CLI.hpp"
#include "version.hh"
#include "ephemeris.hh"
#include "influxpush.hh"
#include "sp3.hh"
static char program[]="reporter";
@ -25,7 +29,7 @@ public:
struct Results
{
vector<double> d;
bool dirty{false};
mutable bool dirty{false};
double median() const
{
return quantile(0.5);
@ -41,7 +45,7 @@ public:
}
};
const Results& done()
const Results& done() const
{
if(results.dirty) {
sort(results.d.begin(), results.d.end());
@ -49,9 +53,12 @@ public:
}
return results;
}
bool empty() const
{
return results.d.empty();
}
private:
Results results;
mutable Results results;
};
/*
@ -81,6 +88,7 @@ struct IntervalStat
bool ripe{false};
bool expired{false};
double rtcmDist{-1};
std::optional<double> rtcmDClock;
};
@ -100,11 +108,15 @@ int main(int argc, char **argv)
CLI::App app(program);
string periodarg("1d");
string beginarg, endarg;
string sp3src("default");
int gnssid=2;
app.add_flag("--version", doVERSION, "show program version and copyright");
app.add_option("--period,-p", periodarg, "period over which to report (1h, 1w)");
app.add_option("--begin,-b", beginarg, "Beginning");
app.add_option("--end,-e", endarg, "End");
app.add_option("--sp3src", sp3src, "Identifier of SP3 source");
app.add_option("--sigid,-s", sigid, "Signal identifier. 1 or 5 for Galileo.");
app.add_option("--gnssid,-g", gnssid, "gnssid, 0 GPS, 2 Galileo");
app.add_option("--influxdb", influxDBName, "Name of influxdb database");
try {
app.parse(argc, argv);
@ -131,7 +143,8 @@ int main(int argc, char **argv)
string url="http://127.0.0.1:8086/query?db="+influxDBName+"&epoch=s&q=";
string query="select distinct(value) from sisa where "+period+" and sigid='"+to_string(sigid)+"' group by gnssid,sv,sigid,time(10m)";
string sisaname = (gnssid==2) ? "sisa" : "gpsura";
string query="select distinct(value) from "+sisaname+" where "+period+" and sigid='"+to_string(sigid)+"' group by gnssid,sv,sigid,time(10m)";
cout<<"query: "<<query<<endl;
auto res = mc.getURL(url + mc.urlEncode(query));
@ -150,8 +163,9 @@ int main(int argc, char **argv)
}
}
res = mc.getURL(url + mc.urlEncode("select distinct(e1bhs) from galhealth where "+period+" and sigid='"+to_string(sigid)+"' group by gnssid,sv,sigid,time(10m)"));
string healthname = (gnssid == 2) ? "galhealth" : "gpshealth";
string healthfieldname = (gnssid==2) ? "e1bhs" : "value";
res = mc.getURL(url + mc.urlEncode("select distinct("+healthfieldname+") from "+healthname+" where "+period+" and sigid='"+to_string(sigid)+"' group by gnssid,sv,sigid,time(10m)"));
j = nlohmann::json::parse(res);
for(const auto& sv : j["results"][0]["series"]) {
@ -161,7 +175,7 @@ int main(int argc, char **argv)
for(const auto& v : sv["values"]) {
auto healthy = (int)v[1];
g_stats[id][(int)v[0]].unhealthy = healthy;
g_stats[id][(int)v[0]].unhealthy = healthy; // hngg
}
}
@ -170,7 +184,7 @@ int main(int argc, char **argv)
for(const auto& sv : j["results"][0]["series"]) {
const auto& tags=sv["tags"];
SatID id{(unsigned int)std::stoi((string)tags["gnssid"]), (unsigned int)std::stoi((string)tags["sv"]), (unsigned int)std::stoi((string)tags["sigid"])};
for(const auto& v : sv["values"]) {
if(v.size() > 1 && v[1] != nullptr) {
@ -187,17 +201,50 @@ int main(int argc, char **argv)
}
}
/////////////////////
///////////////////// rtcm-eph
res = mc.getURL(url + mc.urlEncode("select mean(\"total-dist\") from \"rtcm-eph-correction\" where "+period+" and sigid='"+to_string(sigid)+"' group by gnssid,sv,sigid,time(10m)"));
string rtcmQuery = "select mean(\"radial\") from \"rtcm-eph-correction\" where "+period+" and sigid='"+to_string(sigid)+"' and gnssid='"+to_string(gnssid)+"' group by gnssid,sv,sigid,time(10m)";
cout<<"rtcmquery: "<<rtcmQuery<<endl;
res = mc.getURL(url + mc.urlEncode(rtcmQuery));
j = nlohmann::json::parse(res);
for(const auto& sv : j["results"][0]["series"]) {
try {
const auto& tags=sv["tags"];
SatID id{(unsigned int)std::stoi((string)tags["gnssid"]), (unsigned int)std::stoi((string)tags["sv"]), (unsigned int)std::stoi((string)tags["sigid"])};
for(const auto& v : sv["values"]) {
g_stats[id][(int)v[0]].rtcmDist = v[1];
try {
auto val = (double)v[1]; // might trigger exception
g_stats[id][(int)v[0]].rtcmDist = val;
}
catch(...){ continue; }
}
}
catch(...) {
continue;
}
}
/////////////////////
///////////////////// rtcm-clock
string rtcmClockQuery = "select mean(\"dclock0\") from \"rtcm-clock-correction\" where "+period+" and sigid='"+to_string(sigid)+"' and gnssid='"+to_string(gnssid)+"' group by gnssid,sv,sigid,time(10m)";
cout<<"rtcmquery: "<<rtcmClockQuery<<endl;
res = mc.getURL(url + mc.urlEncode(rtcmClockQuery));
j = nlohmann::json::parse(res);
for(const auto& sv : j["results"][0]["series"]) {
try {
const auto& tags=sv["tags"];
SatID id{(unsigned int)std::stoi((string)tags["gnssid"]), (unsigned int)std::stoi((string)tags["sv"]), (unsigned int)std::stoi((string)tags["sigid"])};
for(const auto& v : sv["values"]) {
try {
auto val = (double) v[1]; // might trigger an exception
g_stats[id][(int)v[0]].rtcmDClock = val;
}
catch(...){ continue; }
}
}
catch(...) {
@ -209,6 +256,253 @@ int main(int argc, char **argv)
/////////////////////
map<SatID, map<time_t, GalileoMessage>> galephs;
map<SatID, map<time_t, GPSState>> gpsephs;
res = mc.getURL(url + mc.urlEncode("select * from \"ephemeris-actual\" where "+period+" and sigid='"+to_string(sigid)+"' and gnssid='"+to_string(gnssid)+"' group by gnssid,sv,sigid"));
j = nlohmann::json::parse(res);
for(const auto& sv : j["results"][0]["series"]) {
try {
const auto& tags=sv["tags"];
SatID id{(unsigned int)std::stoi((string)tags["gnssid"]), (unsigned int)std::stoi((string)tags["sv"]), (unsigned int)std::stoi((string)tags["sigid"])};
// cout << makeSatIDName(id) <<": "<<sv["columns"] << endl;
map<string, int> cmap;
for(const auto& c : sv["columns"]) {
cmap[c]=cmap.size();
}
for(const auto& v : sv["values"]) {
// cout << makeSatIDName(id)<<": crc "<<v[cmap["crc"]]<<" e " <<v[cmap["e"]] << endl;
if(id.gnss==2) {
GalileoMessage gm;
gm.e = v[cmap["e"]];
gm.crc = v[cmap["crc"]];
gm.crs = v[cmap["crs"]];
gm.cuc = v[cmap["cuc"]];
gm.cus = v[cmap["cus"]];
gm.cic = v[cmap["cic"]];
gm.cis = v[cmap["cis"]];
gm.sqrtA = v[cmap["sqrta"]];
gm.t0e = v[cmap["t0e"]];
gm.m0 = v[cmap["m0"]];
gm.deltan = v[cmap["deltan"]];
gm.omega0 = v[cmap["omega0"]];
gm.omegadot = v[cmap["omegadot"]];
gm.idot = v[cmap["idot"]];
gm.omega = v[cmap["omega"]];
gm.i0 = v[cmap["i0"]];
gm.t0e = v[cmap["t0e"]];
gm.iodnav = v[cmap["iod"]];
if(cmap.count("af0")) {
gm.af0 = v[cmap["af0"]];
gm.af1 = v[cmap["af1"]];
gm.af2 = v[cmap["af2"]];
gm.t0c = v[cmap["t0c"]];
}
Point pos;
galephs[id][v[cmap["time"]]]= gm;
}
else if(id.gnss==0) {
GPSState gm{};
gm.e = v[cmap["e"]];
gm.crc = v[cmap["crc"]];
gm.crs = v[cmap["crs"]];
gm.cuc = v[cmap["cuc"]];
gm.cus = v[cmap["cus"]];
gm.cic = v[cmap["cic"]];
gm.cis = v[cmap["cis"]];
gm.sqrtA = v[cmap["sqrta"]];
gm.t0e = v[cmap["t0e"]];
gm.m0 = v[cmap["m0"]];
gm.deltan = v[cmap["deltan"]];
gm.omega0 = v[cmap["omega0"]];
gm.omegadot = v[cmap["omegadot"]];
gm.idot = v[cmap["idot"]];
gm.omega = v[cmap["omega"]];
gm.i0 = v[cmap["i0"]];
gm.t0e = v[cmap["t0e"]];
gm.gpsiod = v[cmap["iod"]];
if(cmap.count("af0")) {
gm.af0 = v[cmap["af0"]];
gm.af1 = v[cmap["af1"]];
gm.af2 = v[cmap["af2"]];
gm.t0c = v[cmap["t0c"]];
}
Point pos;
gpsephs[id][v[cmap["time"]]]= gm;
}
}
}
catch(...) {
continue;
}
}
cout<<"Gathered ephemerides for "<<galephs.size()<<" galileo + "<<gpsephs.size()<<" GPS signals"<<endl;
/////////////////////
map<SatID, map<time_t, SP3Entry>> sp3s;
string spq3="select * from \"sp3\" where "+period+" and sp3src =~ /"+sp3src+"/ and gnssid='"+to_string(gnssid)+"' group by gnssid,sv,sigid";
cout<<"sp3 query: "<<spq3<<endl;
res = mc.getURL(url + mc.urlEncode(spq3));
j = nlohmann::json::parse(res);
cout<<"Gathered sp, got "<< j["results"][0]["series"].size()<< " tags"<<endl;
for(const auto& sv : j["results"][0]["series"]) {
try {
const auto& tags=sv["tags"];
SatID id{(unsigned int)std::stoi((string)tags["gnssid"]), (unsigned int)std::stoi((string)tags["sv"]), (unsigned int)sigid};
// SP3 data does not have a sigid, it refers to the center of mass, not the antenna phase center
// cout << makeSatIDName(id) <<": "<<sv["columns"] << endl;
map<string, int> cmap;
for(const auto& c : sv["columns"]) {
cmap[c]=cmap.size();
}
for(const auto& v : sv["values"]) {
// cout << makeSatIDName(id)<<": time "<<v[cmap["time"]] <<" x "<<v[cmap["x"]]<<" y " <<v[cmap["y"]] << " z " << v[cmap["z"]] << endl;
SP3Entry e;
e.t = v[cmap["time"]]; // UTC!!
e.x = v[cmap["x"]];
e.y = v[cmap["y"]];
e.z = v[cmap["z"]];
e.clockBias = v[cmap["clock-bias"]];
sp3s[id][e.t]=e;
}
}
catch(std::exception& e)
{
cerr<<"Error: "<<e.what()<<endl;
}
}
ofstream csv("sp3.csv");
csv<<"timestamp gnss sv sigid zerror"<<endl;
InfluxPusher idb(influxDBName);
map<SatID, Stats> sp3zerrors, sp3clockerrors;
/////////////////////
for(const auto& svsp3 : sp3s) {
const auto& id = svsp3.first;
const auto& es = svsp3.second;
// cout<<"Looking at SP3 for " << makeSatIDName(id)<<", have "<<es.size()<< " entries"<<endl;
for(const auto& e : es) {
// cout << humanTimeShort(e.first)<<": ";
if(id.gnss == 2) {
const auto& svephs = galephs[id];
auto iter = svephs.lower_bound(e.first);
// this logic is actually sort of wrong & ignores the last ephemeris
if(iter != svephs.end() && iter != svephs.begin()) {
--iter;
// cout << "found ephemeris from "<< humanTimeShort(iter->first)<<" iod "<<iter->second.iodnav;
// our UTC timestamp from SP3 need to be converted into a tow
int offset = id.gnss ? 935280000 : 315964800;
int sp3tow = (e.first - offset) % (7*86400);
Point epos;
getCoordinates(sp3tow, iter->second, &epos);
double clkoffset= iter->second.getAtomicOffset(sp3tow).first - e.second.clockBias;
// cout<<" "<<iter->second.getAtomicOffset(sp3tow).first <<" v " << e.second.clockBias<<" ns ";
// cout << " ("<<epos.x<<", "<<epos.y<<", "<<epos.z<<") - > ("<<e.second.x<<", "<<e.second.y<<", "<<e.second.z<<") " ;
// cout <<" ("<<epos.x - e.second.x<<", "<<epos.y - e.second.y <<", "<<epos.z - e.second.z<<")";
Point sp3pos(e.second.x, e.second.y, e.second.z);
Vector v(epos, sp3pos);
// cout<< " -> " << v.length();
Vector dir(Point(0,0,0), sp3pos);
dir.norm();
Point cv=sp3pos;
cv.x -= 0.519 * dir.x;
cv.y -= 0.519 * dir.y;
cv.z -= 0.519 * dir.z;
Vector v2(epos, cv);
// cout<< " -> " << v2.length();
// cout<<" z-error: "<<dir.inner(v);
csv << e.first << " " << id.gnss <<" " << id.sv << " " << id.sigid <<" " << dir.inner(v) << endl;
idb.addValue({{"gnssid", id.gnss}, {"sv", id.sv}, {"sp3src", sp3src}},
"sp3delta",
{{"ecef-dx", v.x}, {"ecef-dy", v.y}, {"ecef-dz", v.z}, {"sv-dz", dir.inner(v)}, {"dclock", clkoffset},
{"iod-nav",iter->second.iodnav}},
e.first);
sp3clockerrors[id].add(clkoffset); // nanoseconds
sp3zerrors[id].add(100*dir.inner(v) - 80); // meters -> cm
}
}
else if(id.gnss==0) {
const auto& svephs = gpsephs[id];
// this is keyed on the moment of _issue_
auto iter = svephs.lower_bound(e.first);
if(iter != svephs.end() && iter != svephs.begin()) {
--iter;
// cout << "found ephemeris from "<< humanTimeShort(iter->first)<<" iod "<<iter->second.iodnav;
// our UTC timestamp from SP3 need to be converted into a tow
int offset = 315964800;
int sp3tow = (e.first - offset) % (7*86400);
Point epos;
getCoordinates(sp3tow, iter->second, &epos);
double clkoffset= getGPSAtomicOffset(sp3tow, iter->second).first - e.second.clockBias;
// cout<<" "<<iter->second.getAtomicOffset(sp3tow).first <<" v " << e.second.clockBias<<" ns ";
// cout << " ("<<epos.x<<", "<<epos.y<<", "<<epos.z<<") - > ("<<e.second.x<<", "<<e.second.y<<", "<<e.second.z<<") " ;
// cout <<" ("<<epos.x - e.second.x<<", "<<epos.y - e.second.y <<", "<<epos.z - e.second.z<<")";
Point sp3pos(e.second.x, e.second.y, e.second.z);
Vector v(epos, sp3pos);
// cout<< " -> " << v.length();
Vector dir(Point(0,0,0), sp3pos);
dir.norm();
Point cv=sp3pos;
cv.x -= 0.519 * dir.x;
cv.y -= 0.519 * dir.y;
cv.z -= 0.519 * dir.z;
Vector v2(epos, cv);
// cout<< " -> " << v2.length();
// cout<<" z-error: "<<dir.inner(v);
csv << e.first << " " << id.gnss <<" " << id.sv << " " << id.sigid <<" " << dir.inner(v) << endl;
idb.addValue({{"gnssid", id.gnss}, {"sv", id.sv}, {"sp3src", sp3src}},
"sp3delta",
{{"ecef-dx", v.x}, {"ecef-dy", v.y}, {"ecef-dz", v.z}, {"sv-dz", dir.inner(v)}, {"dclock", clkoffset},
{"iod-nav",iter->second.gpsiod}},
e.first);
sp3clockerrors[id].add(clkoffset); // nanoseconds
sp3zerrors[id].add(100*dir.inner(v)); // meters -> cm
}
}
}
}
/////////////////////
g_stats.erase({2,14,1});
g_stats.erase({2,18,1});
g_stats.erase({2,14,5});
@ -225,9 +519,17 @@ int main(int argc, char **argv)
if(sv.second.rbegin()->first > stop)
stop = sv.second.rbegin()->first;
}
if(sv.second.size() > maxintervals)
maxintervals = sv.second.size();
unsigned int liveInterval=0;
for(const auto i : sv.second) {
if(i.second.sisa.has_value())
liveInterval++;
else {
// cout<<makeSatIDName(sv.first)<<": no Sisa, "<< i.second.rtcmDClock.has_value()<<" " << i.second.unhealthy.has_value()<<" " <<i.second.rtcmDist<<" ripe "<<i.second.ripe<<endl;
}
}
if(liveInterval > maxintervals)
maxintervals = liveInterval;
}
cout<<"Report on "<<g_stats.size()<<" SVs from "<<humanTime(start) <<" to " <<humanTime(stop) << endl;
@ -239,13 +541,17 @@ int main(int argc, char **argv)
int napa=0, unhealthy=0, healthy=0, testing=0, ripe=0, expired=0;
Stats rtcm;
for(const auto& i : sv.second) {
Stats rtcm, clockRtcm;
for(const auto& i : sv.second) {
// cout<<humanTimeShort(i.first)<<" "<<((i.second.sisa.has_value()) ? "S" : " ")<<" ";
if(i.second.rtcmDist >= 0) {
rtcm.add(i.second.rtcmDist);
totRTCM.add(i.second.rtcmDist);
}
if(i.second.rtcmDClock)
clockRtcm.add(*i.second.rtcmDClock * 100);
if(i.second.ripe)
ripe++;
@ -273,37 +579,63 @@ int main(int argc, char **argv)
napa++;
}
}
// cout<<endl;
totnapa += napa;
totunhealthy += unhealthy;
tottesting += testing;
tothealthy += healthy;
totripe += ripe;
totexpired += expired;
totunobserved += maxintervals-sv.second.size();
cout<<fmt::sprintf("E%02d: %6.2f%% unobserved, %6.2f%% unhealthy, %6.2f%% healthy, %6.2f%% testing, %6.2f%% napa, %6.2f%% ripe, %6.2f%% expired, %.1f - %.1f - %.1f cm",
sv.first.sv,
100.0*(maxintervals-sv.second.size())/maxintervals,
int liveInterval=0;
for(const auto i : sv.second)
if(i.second.sisa.has_value())
liveInterval++;
totunobserved += maxintervals - liveInterval;
cout<<fmt::sprintf("%s: %6.2f%% unobserved, %6.2f%% unhealthy, %6.2f%% healthy, %6.2f%% testing, %6.2f%% napa, %6.2f%% ripe, %6.2f%% expired",
makeSatPartialName(sv.first),
100.0*(maxintervals-liveInterval)/maxintervals,
100.0*unhealthy/maxintervals,
100.0*healthy/maxintervals,
100.0*testing/maxintervals,
100.0*napa/maxintervals,
100.0*ripe/maxintervals,
100.0*expired/maxintervals,
rtcm.done().quantile(0.10)/10, rtcm.done().median()/10, rtcm.done().quantile(0.9)/10
)<<endl;
100.0*expired/maxintervals);
if(!rtcm.empty())
cout<<fmt::sprintf(", %.1f - %.1f - %.1f cm",
rtcm.done().quantile(0.10)/10, rtcm.done().median()/10, rtcm.done().quantile(0.9)/10);
if(!clockRtcm.empty())
cout<<fmt::sprintf(", c %.1f - %.1f - %.1f cm",
clockRtcm.done().quantile(0.10), clockRtcm.done().median(), clockRtcm.done().quantile(0.9));
if(!sp3zerrors[sv.first].empty()) {
const auto& z = sp3zerrors[sv.first];
cout<<fmt::sprintf(", z %.1f - %.1f - %.1f cm",
z.done().quantile(0.10), z.done().median(), z.done().quantile(0.9));
}
cout<<endl;
}
cout<<"------------------------------------------------------------------------------------------"<<endl;
cout<<fmt::sprintf("Tot: %6.2f%% unobserved, %6.2f%% unhealthy, %6.2f%% healthy, %6.2f%% testing, %6.2f%% napa, %6.2f%% ripe, %6.2f%% expired, %.1f - %.1f - %.1f cm",
cout<<fmt::sprintf("Tot: %6.2f%% unobserved, %6.2f%% unhealthy, %6.2f%% healthy, %6.2f%% testing, %6.2f%% napa, %6.2f%% ripe, %6.2f%% expired",
100.0*(totunobserved)/maxintervals/g_stats.size(),
100.0*totunhealthy/maxintervals/g_stats.size(),
100.0*tothealthy/maxintervals/g_stats.size(),
100.0*tottesting/maxintervals/g_stats.size(),
100.0*totnapa/maxintervals/g_stats.size(),
100.0*totripe/maxintervals/g_stats.size(),
100.0*totexpired/maxintervals/g_stats.size(),
totRTCM.done().quantile(0.10)/10, totRTCM.done().median()/10, totRTCM.done().quantile(0.9)/10
)<<endl;
100.0*totexpired/maxintervals/g_stats.size());
if(!totRTCM.empty())
cout<<fmt::sprintf(", %.1f - %.1f - %.1f cm",
totRTCM.done().quantile(0.10)/10, totRTCM.done().median()/10, totRTCM.done().quantile(0.9)/10);
cout<<endl;
}