GeoidEval.cpp

Go to the documentation of this file.
00001 /**
00002  * \file GeoidEval.cpp
00003  * \brief Command line utility for evaluation geoid heights
00004  *
00005  * Copyright (c) Charles Karney (2009, 2010) <charles@karney.com>
00006  * and licensed under the LGPL.  For more information, see
00007  * http://geographiclib.sourceforge.net/
00008  *
00009  * Compile with -I../include and link with Geoid.o DMS.o
00010  *
00011  * See \ref geoideval for usage information.
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 }

Generated on 21 May 2010 for GeographicLib by  doxygen 1.6.1