GeoidEval.cpp
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014 #include "GeographicLib/Geoid.hpp"
00015 #include "GeographicLib/DMS.hpp"
00016 #include "GeographicLib/GeoCoords.hpp"
00017 #include <iostream>
00018 #include <sstream>
00019 #include <iomanip>
00020
00021 int usage(int retval) {
00022 using namespace GeographicLib;
00023 std::string
00024 geoidpath = Geoid::GeoidPath(),
00025 defaultpath = Geoid::DefaultPath();
00026 if (geoidpath.empty())
00027 geoidpath = "UNDEFINED";
00028 ( retval ? std::cerr : std::cout ) <<
00029 "Usage:\n\
00030 GeoidEval [-n name] [-d dir] [-l] [-a] [-c south west north east] [-v] [-h]\n\
00031 $Id: GeoidEval.cpp 6827 2010-05-20 19:56:18Z karney $\n\
00032 \n\
00033 Read in positions on standard input and print out the corresponding\n\
00034 geoid heights on standard output. In addition print the northly and\n\
00035 easterly gradients of the geoid height (i.e., the rate at which the\n\
00036 geoid height changes per unit distance along the WGS84 ellipsoid in\n\
00037 the specified directions).\n\
00038 \n\
00039 Positions are given as latitude and longitude, UTM/UPS, or MGRS, in\n\
00040 any of the formats accepted by GeoConvert.\n\
00041 \n\
00042 By default the EGM96 geoid is used with a 5\' grid. This may be\n\
00043 overriden with the -n option. The name specified should be one of\n\
00044 \n\
00045 bilinear error cubic error\n\
00046 name geoid grid max rms max rms\n\
00047 egm84-30 EGM84 30\' 1.546m 70mm 0.274m 14mm\n\
00048 egm84-15 EGM84 15\' 0.413m 18mm 0.020m 1mm\n\
00049 egm96-15 EGM96 15\' 1.152m 40mm 0.169m 7mm\n\
00050 egm96-5 EGM96 5\' 0.140m 5mm 0.003m 1mm\n\
00051 egm2008-5 EGM2008 5\' 0.478m 12mm 0.294m 5mm\n\
00052 egm2008-2_5 EGM2008 2.5\' 0.135m 3mm 0.031m 1mm\n\
00053 egm2008-1 EGM2008 1\' 0.025m 1mm 0.003m 1mm\n\
00054 \n\
00055 (Some of the geoids may not be available.) The errors listed here\n\
00056 are estimates of the quantization and interpolation errors in the\n\
00057 reported heights compared to the specified geoid.\n\
00058 \n\
00059 Cubic interpolation is used to compute the geoid height unless\n\
00060 -l is specified in which case bilinear interpolation is used.\n\
00061 Cubic interpolation is more accurate; however it results in\n\
00062 small discontinuities in the returned height on cell boundaries.\n\
00063 The gradients are computed by differentiating the interpolated\n\
00064 results.\n\
00065 \n\
00066 GeoidEval will load the geoid data from the directory specified by\n\
00067 the -d option. If this is not provided, it will look up the value of\n\
00068 GEOID_PATH (currently " << geoidpath << ") in the environment.\n\
00069 If this is not defined, it will use the compile-time value of\n"
00070 << defaultpath << ".\n\
00071 \n\
00072 By default, the data file is randomly read to compute the geoid\n\
00073 heights at the input positions. Usually this is sufficient for\n\
00074 interactive use. If many heights are to be computed, GeoidEval\n\
00075 allows a block of data to be read into memory and heights within the\n\
00076 corresponding rectangle can then be computed without any disk acces.\n\
00077 If -a is specified all the geoid data is read; in the case of\n\
00078 egm2008-1, this requires about 0.5 GB of RAM. The -c option allows\n\
00079 a rectangle of data to be cached. The evaluation of heights\n\
00080 outside the cached rectangle causes the necessary data to be read\n\
00081 from disk.\n\
00082 \n\
00083 Regardless of whether any cache is requested (with the -a or -c\n\
00084 options), the data for the last grid cell in cached. This allows\n\
00085 the geoid height along a continuous path to be returned with little\n\
00086 disk overhead.\n\
00087 \n\
00088 The -v option causes the data about the current geoid to be printed\n\
00089 to standard error.\n\
00090 \n\
00091 -h prints this help.\n";
00092 return retval;
00093 }
00094
00095 int main(int argc, char* argv[]) {
00096 using namespace GeographicLib;
00097 typedef Math::real real;
00098 bool cacheall = false, cachearea = false, verbose = false, cubic = true;
00099 real caches, cachew, cachen, cachee;
00100 std::string dir;
00101 std::string geoid = "egm96-5";
00102 for (int m = 1; m < argc; ++m) {
00103 std::string arg(argv[m]);
00104 if (arg == "-a") {
00105 cacheall = true;
00106 cachearea = false;
00107 }
00108 else if (arg == "-c") {
00109 if (m + 4 >= argc) return usage(1);
00110 cacheall = false;
00111 cachearea = true;
00112 try {
00113 DMS::DecodeLatLon(std::string(argv[m + 1]), std::string(argv[m + 2]),
00114 caches, cachew);
00115 DMS::DecodeLatLon(std::string(argv[m + 3]), std::string(argv[m + 4]),
00116 cachen, cachee);
00117 }
00118 catch (const std::exception& e) {
00119 std::cerr << "Error decoding argument of -c: " << e.what() << "\n";
00120 return 1;
00121 }
00122 m += 4;
00123 } else if (arg == "-n") {
00124 if (++m == argc) return usage(1);
00125 geoid = argv[m];
00126 } else if (arg == "-d") {
00127 if (++m == argc) return usage(1);
00128 dir = argv[m];
00129 } else if (arg == "-l") {
00130 cubic = false;
00131 } else if (arg == "-v")
00132 verbose = true;
00133 else
00134 return usage(arg != "-h");
00135 }
00136
00137 int retval = 0;
00138 try {
00139 Geoid g(geoid, dir, cubic);
00140 try {
00141 if (cacheall)
00142 g.CacheAll();
00143 else if (cachearea)
00144 g.CacheArea(caches, cachew, cachen, cachee);
00145 }
00146 catch (const std::exception& e) {
00147 std::cerr << "ERROR: " << e.what() << "\nProceeding without a cache\n";
00148 }
00149 if (verbose) {
00150 std::cerr << "Geoid file: " << g.GeoidFile() << "\n"
00151 << "Description: " << g.Description() << "\n"
00152 << "Interpolation: " << g.Interpolation() << "\n"
00153 << "Date & Time: " << g.DateTime() << "\n"
00154 << "Offset (m): " << g.Offset() << "\n"
00155 << "Scale (m): " << g.Scale() << "\n"
00156 << "Max error (m): " << g.MaxError() << "\n"
00157 << "RMS error (m): " << g.RMSError() << "\n";
00158 if (g.Cache())
00159 std::cerr << "Caching:"
00160 << "\n SW Corner: " << g.CacheSouth() << " " << g.CacheWest()
00161 << "\n NE Corner: " << g.CacheNorth() << " " << g.CacheEast()
00162 << "\n";
00163 }
00164
00165 std::cout << std::fixed;
00166 GeoCoords p;
00167 std::string s;
00168 while (std::getline(std::cin, s)) {
00169 try {
00170 p.Reset(s);
00171 real gradn, grade;
00172 real h = g(p.Latitude(), p.Longitude(), gradn, grade);
00173 std::cout << std::setprecision(4) << h << " " << std::setprecision(2)
00174 << gradn * 1e6 << "e-6 " << grade * 1e6 << "e-6\n";
00175 }
00176 catch (const std::exception& e) {
00177 std::cout << "ERROR: " << e.what() << "\n";
00178 retval = 1;
00179 }
00180 }
00181 }
00182 catch (const std::exception& e) {
00183 std::cerr << "Error reading " << geoid << ": " << e.what() << "\n";
00184 retval = 1;
00185 }
00186 return retval;
00187 }