File: gnssGriddedDataTimeSeries2Ionex.cpp

package info (click to toggle)
groops 0%2Bgit20250907%2Bds-1
  • links: PTS, VCS
  • area: non-free
  • in suites: forky, sid
  • size: 11,140 kB
  • sloc: cpp: 135,607; fortran: 1,603; makefile: 20
file content (173 lines) | stat: -rw-r--r-- 8,527 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
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
/***********************************************/
/**
* @file programTemplate.cpp
*
* @brief Converts TEC maps from GROOPS GriddedDataTimeSeries format to IGS IONEX file format.
*
* @author Sebastian Strasser
* @date 2021-09-15
*/
/***********************************************/

// Latex documentation
#define DOCSTRING docstring
static const char *docstring = R"(
Converts TEC maps from GROOPS \file{gridded data time series}{griddedDataTimeSeries} format
to IGS \href{https://files.igs.org/pub/data/format/ionex1.pdf}{IONEX file format}.

Currently only supports 2D TEC maps.

See also \program{GnssIonex2GriddedDataTimeSeries}, \configClass{IonosphereMap}{gnssParametrizationType:ionosphereMap}.
)";

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

#include "programs/program.h"
#include "inputOutput/system.h"
#include "classes/timeSeries/timeSeries.h"
#include "files/fileGriddedDataTimeSeries.h"
#include "parser/dataVariables.h"

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

/** @brief Converts TEC maps from GROOPS GriddedDataTimeSeries format to IGS IONEX file format.
* @ingroup programsGroup */
class GnssGriddedDataTimeSeries2Ionex
{
public:
  void run(Config &config, Parallel::CommunicatorPtr comm);
};

GROOPS_REGISTER_PROGRAM(GnssGriddedDataTimeSeries2Ionex, SINGLEPROCESS, "Converts TEC maps from GROOPS GriddedDataTimeSeries format to IGS IONEX file format.", Conversion, Gnss, Grid, TimeSplines)

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

void GnssGriddedDataTimeSeries2Ionex::run(Config &config, Parallel::CommunicatorPtr /*comm*/)
{
  try
  {
    FileName fileNameIn, fileNameOut;
    TimeSeriesPtr timeSeries;
    std::string program, institution, mappingFunction, observablesUsed;
    std::vector<std::string> comment, description;
    Double elevationCutoff;
    Int exponent;
    ExpressionVariablePtr expression;

    readConfig(config, "outputfileIonex",                fileNameOut,     Config::MUSTSET,  "", "");
    readConfig(config, "inputfileGriddedDataTimeSeries", fileNameIn,      Config::MUSTSET,  "", "must contain regular grid");
    readConfig(config, "value",                          expression,      Config::MUSTSET,  "data0", "expression (e.g. data column)");
    readConfig(config, "timeSeries",                     timeSeries,      Config::OPTIONAL, "", "(empty = use input file time series)");
    readConfig(config, "program",                        program,         Config::MUSTSET,  "GROOPS", "name of program (for first line)");
    readConfig(config, "institution",                    institution,     Config::MUSTSET,  "TUG (TU Graz)", "name of agency (for first line)");
    readConfig(config, "description",                    description,     Config::OPTIONAL, "", "description in header");
    readConfig(config, "comment",                        comment,         Config::OPTIONAL, R"(["TEC values in 0.1 TECU, 9999 if no value available"])", "comment in header");
    readConfig(config, "mappingFunction",                mappingFunction, Config::DEFAULT,  "NONE", "see IONEX documentation");
    readConfig(config, "elevationCutoff",                elevationCutoff, Config::DEFAULT,  "0", "see IONEX documentation (0 if unknown)");
    readConfig(config, "observablesUsed",                observablesUsed, Config::OPTIONAL, "", "see IONEX documentation");
    readConfig(config, "exponent",                       exponent,        Config::DEFAULT,  "-1", "factor 10^exponent is applied to all values");
    if(isCreateSchema(config)) return;

    logStatus<<"read gridded data time series file <"<<fileNameIn<<">"<<Log::endl;
    InFileGriddedDataTimeSeries infile(fileNameIn);
    GriddedDataRectangular grid;
    if(!grid.init(infile.grid()))
      throw(Exception("grid must be rectangular"));

    std::vector<Time> times = infile.times();
    if(timeSeries)
      times = timeSeries->times();

    logStatus << "write IONEX file <" << fileNameOut << ">" << Log::endl;
    OutFile outfile(fileNameOut);

    // Header: IONEX VERSION / TYPE
    outfile << "     1.0            IONOSPHERE MAPS     GNSS                IONEX VERSION / TYPE" << std::endl;

    // Header: PGM / RUN BY / DATE
    program.resize(20, ' ');
    institution.resize(20, ' ');
    outfile << program << institution << System::now()%"%D %H:%M:%S "s << "PGM / RUN BY / DATE" << std::endl;

    // Header: DESCRIPTION
    for(auto &d : description)
    {
      d.resize(60, ' ');
      outfile << d << "DESCRIPTION" << std::endl;
    }

    // Header: EPOCH OF FIRST/LAST MAP, INTERVAL, # OF MAPS IN FILE
    outfile << times.front()%"  %y    %m    %d    %H    %M    %S"s << std::string(24, ' ') << "EPOCH OF FIRST MAP" << std::endl;
    outfile << times.back()%"  %y    %m    %d    %H    %M    %S"s << std::string(24, ' ') << "EPOCH OF LAST MAP" << std::endl;
    outfile << medianSampling(times).seconds()%"% 6i"s << std::string(54, ' ') << "INTERVAL" << std::endl;
    outfile << times.size()%"% 6i"s << std::string(54, ' ') << "# OF MAPS IN FILE" << std::endl;

    // Header: MAPPING FUNCTION, ELEVATION CUTOFF, OBSERVABLES USED, BASE RADIUS, MAP DIMENSION
    mappingFunction.resize(4, ' ');
    observablesUsed.resize(60, ' ');
    outfile << "  " << mappingFunction << std::string(54, ' ') << "MAPPING FUNCTION" << std::endl;
    outfile << elevationCutoff%"%8.1f"s << std::string(52, ' ') << "ELEVATION CUTOFF" << std::endl;
    outfile << observablesUsed << "OBSERVABLES USED" << std::endl;
    outfile << (infile.grid().ellipsoid.a()/1000)%"%8.1f"s << std::string(52, ' ') << "BASE RADIUS" << std::endl;
    outfile << "     2" << std::string(54, ' ') << "MAP DIMENSION" << std::endl;

    // Header: HGT1 / HGT 2 / DHGT, LAT1 / LAT2 / DLAT, LON1 / LON2 / DLON
    const Double height = mean(Vector(grid.heights));
    UInt columns = (grid.longitudes.front()+2*PI == grid.longitudes.back() ? grid.longitudes.size()-1 : grid.longitudes.size());
    outfile << "  " << (height/1000)%"%6.1f%6.1f"s << 0%"%6.1f"s << std::string(40, ' ') << "HGT1 / HGT2 / DHGT" << std::endl;
    outfile << "  " << (grid.latitudes.front()*RAD2DEG)%"%6.1f"s << (grid.latitudes.back()*RAD2DEG)%"%6.1f"s
            << ((grid.latitudes.back()-grid.latitudes.front())/grid.latitudes.size()*RAD2DEG)%"%6.1f"s << std::string(40, ' ') << "LAT1 / LAT2 / DLAT" << std::endl;
    outfile << "  " << (grid.longitudes.front()*RAD2DEG)%"%6.1f"s << (grid.longitudes.back()*RAD2DEG)%"%6.1f"s
            << ((grid.longitudes.back()-grid.longitudes.front())/columns*RAD2DEG)%"%6.1f"s << std::string(40, ' ') << "LON1 / LON2 / DLON" << std::endl;

    // Header: EXPONENT
    outfile << exponent%"% 6i"s << std::string(54, ' ') << "EXPONENT" << std::endl;

    // Header: COMMENT
    for(auto &c : comment)
    {
      c.resize(60, ' ');
      outfile << c << "COMMENT" << std::endl;
    }

    // Header: END OF HEADER
    outfile << std::string(60, ' ') << "END OF HEADER" << std::endl;

    // Data records
    const Double factor = std::pow(10, exponent);
    for(UInt idEpoch = 0; idEpoch < times.size(); idEpoch++)
    {
      outfile << (idEpoch+1)%"% 6i"s << std::string(54, ' ') << "START OF TEC MAP" << std::endl;
      outfile << times.at(idEpoch)%"  %y    %m    %d    %H    %M    %S"s << std::string(24, ' ') << "EPOCH OF CURRENT MAP" << std::endl;
      for(UInt i = 0; i < grid.latitudes.size(); i++)
      {
        outfile << "  " << (grid.latitudes.at(i)*RAD2DEG)%"% 6.1f"s <<  (grid.longitudes.front()*RAD2DEG)%"%6.1f"s << (grid.longitudes.back()*RAD2DEG)%"%6.1f"s
                << ((grid.longitudes.back()-grid.longitudes.front())/columns*RAD2DEG)%"%6.1f"s << (height/1000)%"%6.1f"s << std::string(28, ' ') << "LAT/LON1/LON2/DLON/H" << std::endl;

        Matrix data = infile.data(times.at(idEpoch));
        VariableList varList;
        addDataVariables(data, varList);

        for(UInt j = 0; j < grid.longitudes.size(); j++)
        {
          evaluateDataVariables(data, i*grid.longitudes.size()+j, varList);
          Double value = expression->evaluate(varList);

          if(j > 0 && j%16 == 0) // start new line after 16 values
            outfile << std::endl;
          outfile << (std::isnan(value) ? " 9999" : (value/factor)%"% 5i"s);
        }
        outfile << std::endl;
      }
      outfile << (idEpoch+1)%"% 6i"s << std::string(54, ' ') << "END OF TEC MAP" << std::endl;
    }

    outfile << std::string(60, ' ') << "END OF FILE" << std::endl;
  }
  catch(std::exception &e)
  {
    GROOPS_RETHROW(e)
  }
}

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