File: gnssProcessingStepWriteResiduals.h

package info (click to toggle)
groops 0%2Bgit20240830%2Bds-1
  • links: PTS, VCS
  • area: non-free
  • in suites: trixie
  • size: 11,052 kB
  • sloc: cpp: 134,939; fortran: 1,569; makefile: 20
file content (155 lines) | stat: -rw-r--r-- 6,592 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
/***********************************************/
/**
* @file gnssProcessingStepWriteResiduals.h
*
* @brief GNSS processing step: WriteResiduals.
*
* @author Torsten Mayer-Guerr
* @date 2021-09-05
*
*/
/***********************************************/

#ifndef __GROOPS_GNSSPROCESSINGSTEPWRITERESIDUALS__
#define __GROOPS_GNSSPROCESSINGSTEPWRITERESIDUALS__

// Latex documentation
#ifdef DOCSTRING_GnssProcessingStep
static const char *docstringGnssProcessingStepWriteResiduals = R"(
\subsection{WriteResiduals}\label{gnssProcessingStepType:writeResiduals}
Writes the \file{observation residuals}{instrument} for all
\configClass{selectReceivers}{platformSelectorType}.
For for each station a file is written. The file name is interpreted as
a template with the variable \verb|{station}| being replaced by the station name.
)";
#endif

/***********************************************/

#include "config/config.h"
#include "classes/platformSelector/platformSelector.h"
#include "gnss/gnssProcessingStep/gnssProcessingStep.h"

/***** CLASS ***********************************/

/** @brief GNSS processing step: WriteResiduals.
* @ingroup gnssProcessingStepGroup
* @see GnssProcessingStep */
class GnssProcessingStepWriteResiduals : public GnssProcessingStepBase
{
  PlatformSelectorPtr selectReceivers;
  FileName            fileNameResiduals;

public:
  GnssProcessingStepWriteResiduals(Config &config);
  void process(GnssProcessingStep::State &state) override;
};

/***********************************************/

inline GnssProcessingStepWriteResiduals::GnssProcessingStepWriteResiduals(Config &config)
{
  try
  {
    readConfig(config, "selectReceivers",     selectReceivers,   Config::MUSTSET, "", "subset of used stations");
    readConfig(config, "outputfileResiduals", fileNameResiduals, Config::MUSTSET, "output/residuals_{loopTime:%D}.{station}.dat", "variable {station} available");
  }
  catch(std::exception &e)
  {
    GROOPS_RETHROW(e)
  }
}

/***********************************************/

inline void GnssProcessingStepWriteResiduals::process(GnssProcessingStep::State &state)
{
  try
  {
    auto selectedReceivers = state.gnss->selectReceivers(selectReceivers);
    VariableList fileNameVariableList;
    fileNameVariableList.setVariable("station", "****");
    logStatus<<"write residuals to file <"<<fileNameResiduals(fileNameVariableList)<<">"<<Log::endl;
    for(auto recv : state.gnss->receivers)
      if(selectedReceivers.at(recv->idRecv()) && state.normalEquationInfo.estimateReceiver.at(recv->idRecv()) && recv->isMyRank())
      {
        GnssReceiverArc arc;
        for(UInt idEpoch : state.normalEquationInfo.idEpochs)
          if(recv->useable(idEpoch))
          {
            // get types
            std::vector<GnssType> types;
            for(UInt idTrans=0; idTrans<recv->idTransmitterSize(idEpoch); idTrans++)
              if(recv->observation(idTrans, idEpoch) && state.gnss->transmitters.at(idTrans)->useable(idEpoch))
                for(UInt idType=0; idType<recv->observation(idTrans, idEpoch)->size(); idType++)
                  if(!recv->observation(idTrans, idEpoch)->at(idType).type.isInList(types))
                    types.push_back(recv->observation(idTrans, idEpoch)->at(idType).type & ~(GnssType::PRN+GnssType::FREQ_NO));
            std::sort(types.begin(), types.end());
            if(!types.size())
              continue;

            GnssReceiverEpoch epoch;
            GnssType system = GnssType::SYSTEM;
            for(UInt idType=0; idType<types.size(); idType++)
            {
              if(types.at(idType) != system)
              {
                system = types.at(idType) & GnssType::SYSTEM;
                epoch.obsType.push_back( GnssType::AZIMUT    + GnssType::L1 + system );
                epoch.obsType.push_back( GnssType::ELEVATION + GnssType::L1 + system );
                epoch.obsType.push_back( GnssType::AZIMUT    + GnssType::L2 + system );
                epoch.obsType.push_back( GnssType::ELEVATION + GnssType::L2 + system );
                epoch.obsType.push_back( GnssType::IONODELAY + system );
              }
              // residuals, redundancy, sigma
              epoch.obsType.insert(epoch.obsType.end(), {types.at(idType), types.at(idType), types.at(idType)});
            }

            for(UInt idTrans=0; idTrans<recv->idTransmitterSize(idEpoch); idTrans++)
              if(recv->observation(idTrans, idEpoch) && state.gnss->transmitters.at(idTrans)->useable(idEpoch))
              {
                const GnssObservation &obs = *recv->observation(idTrans, idEpoch);
                const GnssObservationEquation eqn(obs, *recv, *state.gnss->transmitters.at(idTrans),
                                                  state.gnss->funcRotationCrf2Trf, state.gnss->funcReduceModels, idEpoch, FALSE, {});
                const GnssType prn = obs.at(0).type & (GnssType::SYSTEM + GnssType::PRN + GnssType::FREQ_NO);
                UInt idType = std::distance(epoch.obsType.begin(), std::find(epoch.obsType.begin(), epoch.obsType.end(), prn));
                if(idType >= epoch.obsType.size())
                  continue;

                epoch.time = eqn.timeRecv;
                epoch.satellite.push_back(prn);
                epoch.observation.insert(epoch.observation.end(), {eqn.azimutRecvAnt, eqn.elevationRecvAnt, eqn.azimutTrans, eqn.elevationTrans, obs.STEC});
                idType += 5;

                for(; (idType<epoch.obsType.size()) && (epoch.obsType.at(idType) == prn); idType+=3)
                {
                  epoch.observation.insert(epoch.observation.end(), {0., 0., 1.});
                  for(UInt i=0; i<obs.size(); i++)
                    if(obs.at(i).type == epoch.obsType.at(idType))
                    {
                      epoch.observation.at(epoch.observation.size()-3) = obs.at(i).residuals;
                      epoch.observation.at(epoch.observation.size()-2) = obs.at(i).redundancy;
                      epoch.observation.at(epoch.observation.size()-1) = obs.at(i).sigma/obs.at(i).sigma0;
                      break;
                    }
                }
              } // for(idTrans)

            if(epoch.satellite.size())
              arc.push_back(epoch);
          } // for(idEpoch)

        VariableList fileNameVariableList;
        fileNameVariableList.setVariable("station", recv->name());
        InstrumentFile::write(fileNameResiduals(fileNameVariableList), arc);
      } // for(recv)
  }
  catch(std::exception &e)
  {
    GROOPS_RETHROW(e)
  }
}

/***********************************************/

#endif