GeoidToGTX.cpp
Go to the documentation of this file.
1 // Write out a gtx file of geoid heights above the ellipsoid. For egm2008 at
2 // 1' resolution this takes about 40 mins on a 8-processor Intel 2.66 GHz
3 // machine using OpenMP.
4 //
5 // For the format of gtx files, see
6 // http://vdatum.noaa.gov/dev/gtx_info.html#dev_gtx_binary
7 //
8 // data is binary big-endian:
9 // south latitude edge (degrees double)
10 // west longitude edge (degrees double)
11 // delta latitude (degrees double)
12 // delta longitude (degrees double)
13 // nlat = number of latitude rows (integer)
14 // nlong = number of longitude columns (integer)
15 // nlat * nlong geoid heights (meters float)
16 
17 #include <vector>
18 #include <iostream>
19 #include <fstream>
20 #include <string>
21 #include <algorithm>
22 
23 #if defined(_OPENMP)
24 #define HAVE_OPENMP 1
25 #else
26 #define HAVE_OPENMP 0
27 #endif
28 
29 #if HAVE_OPENMP
30 # include <omp.h>
31 #endif
32 
36 
37 using namespace std;
38 using namespace GeographicLib;
39 
40 int main(int argc, const char* const argv[]) {
41  // Hardwired for 3 args:
42  // 1 = the gravity model (e.g., egm2008)
43  // 2 = intervals per degree
44  // 3 = output GTX file
45  if (argc != 4) {
46  cerr << "Usage: " << argv[0]
47  << " gravity-model intervals-per-degree output.gtx\n";
48  return 1;
49  }
50  try {
51  // Will need to set the precision for each thread, so save return value
52  int ndigits = Utility::set_digits();
53  string model(argv[1]);
54  // Number of intervals per degree
55  int ndeg = Utility::val<int>(string(argv[2]));
56  string filename(argv[3]);
57  GravityModel g(model);
58  int
59  nlat = 180 * ndeg + 1,
60  nlon = 360 * ndeg;
62  delta = 1 / Math::real(ndeg), // Grid spacing
63  latorg = -90,
64  lonorg = -180;
65  // Write results as floats in binary mode
66  ofstream file(filename.c_str(), ios::binary);
67 
68  // Write header
69  {
70  Math::real transform[] = {latorg, lonorg, delta, delta};
71  unsigned sizes[] = {unsigned(nlat), unsigned(nlon)};
72  Utility::writearray<double, Math::real, true>(file, transform, 4);
73  Utility::writearray<unsigned, unsigned, true>(file, sizes, 2);
74  }
75 
76  // Compute and store results for nbatch latitudes at a time
77  const int nbatch = 64;
78  vector< vector<float> > N(nbatch, vector<float>(nlon));
79 
80  for (int ilat0 = 0; ilat0 < nlat; ilat0 += nbatch) { // Loop over batches
81  int nlat0 = min(nlat, ilat0 + nbatch);
82 
83 #if HAVE_OPENMP
84 # pragma omp parallel for
85 #endif
86  for (int ilat = ilat0; ilat < nlat0; ++ilat) { // Loop over latitudes
87  Utility::set_digits(ndigits); // Set the precision
89  lat = latorg + (ilat / ndeg) + delta * (ilat - ndeg * (ilat / ndeg)),
90  h = 0;
91  GravityCircle c(g.Circle(lat, h, GravityModel::GEOID_HEIGHT));
92  for (int ilon = 0; ilon < nlon; ++ilon) { // Loop over longitudes
93  Math::real lon = lonorg
94  + (ilon / ndeg) + delta * (ilon - ndeg * (ilon / ndeg));
95  N[ilat - ilat0][ilon] = float(c.GeoidHeight(lon));
96  } // longitude loop
97  } // latitude loop -- end of parallel section
98 
99  for (int ilat = ilat0; ilat < nlat0; ++ilat) // write out data
100  Utility::writearray<float, float, true>(file, N[ilat - ilat0]);
101  } // batch loop
102  }
103  catch (const exception& e) {
104  cerr << "Caught exception: " << e.what() << "\n";
105  return 1;
106  }
107  catch (...) {
108  cerr << "Caught unknown exception\n";
109  return 1;
110  }
111 }
float real
Definition: datatypes.h:10
static const double lat
noiseModel::Diagonal::shared_ptr model
Header for GeographicLib::Utility class.
#define min(a, b)
Definition: datatypes.h:19
std::vector< Array2i > sizes
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
Definition: BFloat16.h:88
Header for GeographicLib::GravityModel class.
#define N
Definition: gksort.c:12
void g(const string &key, int i)
Definition: testBTree.cpp:41
int main(int argc, const char *const argv[])
Definition: GeoidToGTX.cpp:40
Namespace for GeographicLib.
Array< double, 1, 3 > e(1./3., 0.5, 2.)
EIGEN_DONT_INLINE void transform(const Transformation &t, Data &data)
Definition: geometry.cpp:25
Model of the earth&#39;s gravity field.
const double h
GravityCircle Circle(real lat, real h, unsigned caps=ALL) const
static const double lon
Header for GeographicLib::GravityCircle class.
Gravity on a circle of latitude.


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:18