coursepointer 0.3.4

Converts waypoints into Garmin FIT course points
Documentation 
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

// Write out a gtx file of geoid heights above the ellipsoid.  For egm2008 at
// 1' resolution this takes about 10 mins on a 8-processor Intel 3.0 GHz
// machine using OpenMP.
//
// For the format of gtx files, see
// https://vdatum.noaa.gov/docs/gtx_info.html#dev_gtx_binary
//
// data is binary big-endian:
//   south latitude edge (degrees double)
//   west longitude edge (degrees double)
//   delta latitude (degrees double)
//   delta longitude (degrees double)
//   nlat = number of latitude rows (integer)
//   nlong = number of longitude columns (integer)
//   nlat * nlong geoid heights (meters float)

#include <vector>
#include <iostream>
#include <fstream>
#include <string>
#include <algorithm>

#if defined(_OPENMP)
#define HAVE_OPENMP 1
#else
#define HAVE_OPENMP 0
#endif

#if HAVE_OPENMP
#  include <omp.h>
#endif

#include <GeographicLib/GravityModel.hpp>
#include <GeographicLib/GravityCircle.hpp>
#include <GeographicLib/Utility.hpp>

using namespace std;
using namespace GeographicLib;

int main(int argc, const char* const argv[]) {
  // Hardwired for 3 args:
  // 1 = the gravity model (e.g., egm2008)
  // 2 = intervals per degree
  // 3 = output GTX file
  if (argc != 4) {
    cerr << "Usage: " << argv[0]
         << " gravity-model intervals-per-degree output.gtx\n";
    return 1;
  }
  try {
    // Will need to set the precision for each thread, so save return value
    int ndigits = Utility::set_digits();
    string model(argv[1]);
    // Number of intervals per degree
    int ndeg = Utility::val<int>(string(argv[2]));
    string filename(argv[3]);
    GravityModel g(model);
    int
      nlat = Math::hd * ndeg + 1,
      nlon = Math::td * ndeg;
    Math::real
      delta = 1 / Math::real(ndeg), // Grid spacing
      latorg = -Math::qd,
      lonorg = -Math::hd;
    // Write results as floats in binary mode
    ofstream file(filename.c_str(), ios::binary);

    // Write header
    {
      Math::real transform[] = {latorg, lonorg,
                                1/Math::real(ndeg), 1/Math::real(ndeg)};
      unsigned sizes[] = {unsigned(nlat), unsigned(nlon)};
      Utility::writearray<double, Math::real, true>(file, transform, 4);
      Utility::writearray<unsigned, unsigned, true>(file, sizes, 2);
    }

    // Compute and store results for nbatch latitudes at a time
    const int nbatch = 64;
    vector< vector<float> > N(nbatch, vector<float>(nlon));

    for (int ilat0 = 0; ilat0 < nlat; ilat0 += nbatch) { // Loop over batches
      int nlat0 = min(nlat, ilat0 + nbatch);

#if HAVE_OPENMP
#  pragma omp parallel for
#endif
      for (int ilat = ilat0; ilat < nlat0; ++ilat) { // Loop over latitudes
        Utility::set_digits(ndigits);                // Set the precision
        Math::real lat = (latorg * ndeg + ilat) / ndeg, h = 0;
        GravityCircle c(g.Circle(lat, h, GravityModel::GEOID_HEIGHT));
        for (int ilon = 0; ilon < nlon; ++ilon) { // Loop over longitudes
          Math::real lon = (lonorg * ndeg + ilon) / ndeg;
          N[ilat - ilat0][ilon] = float(c.GeoidHeight(lon));
        } // longitude loop
      }   // latitude loop -- end of parallel section

      for (int ilat = ilat0; ilat < nlat0; ++ilat) // write out data
        Utility::writearray<float, float, true>(file, N[ilat - ilat0]);
    } // batch loop
  }
  catch (const exception& e) {
    cerr << "Caught exception: " << e.what() << "\n";
    return 1;
  }
  catch (...) {
    cerr << "Caught unknown exception\n";
    return 1;
  }
}