torstenandCopilot
342ffd3e5c
- HandlerMinelayerData: BUILD_RTI Guards für RTI-Abhängigkeiten - svcTabCommon: Export-Header-Fix via Fallback-Target - TrainingAreaBase: Stub-Header für BUILD_RTI=OFF - ddsfts: Include-Pfad und Dependency-Fixes - BUILD_RTI=OFF Compilation vollständig unterstützt Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com> Entire-Checkpoint: 01M0JNZA586SACX1PRZ2DJAE5S
110 lines
3.4 KiB
C++
110 lines
3.4 KiB
C++
// 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;
|
|
}
|
|
}
|