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, 2011) <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 the <a href="GeoidEval.1.html">man page</a> for usage
00012  * information.
00013  **********************************************************************/
00014 
00015 #include "GeographicLib/Geoid.hpp"
00016 #include "GeographicLib/DMS.hpp"
00017 #include "GeographicLib/GeoCoords.hpp"
00018 #include <iostream>
00019 
00020 #include "GeoidEval.usage"
00021 
00022 int main(int argc, char* argv[]) {
00023   using namespace GeographicLib;
00024   typedef Math::real real;
00025   bool cacheall = false, cachearea = false, verbose = false, cubic = true;
00026   real caches, cachew, cachen, cachee;
00027   std::string dir;
00028   std::string geoid = Geoid::DefaultGeoidName();
00029   std::string zone;
00030   Geoid::convertflag heightmult = Geoid::NONE;
00031   for (int m = 1; m < argc; ++m) {
00032     std::string arg(argv[m]);
00033     if (arg == "-a") {
00034       cacheall = true;
00035       cachearea = false;
00036     }
00037     else if (arg == "-c") {
00038       if (m + 4 >= argc) return usage(1, true);
00039       cacheall = false;
00040       cachearea = true;
00041       try {
00042         DMS::DecodeLatLon(std::string(argv[m + 1]), std::string(argv[m + 2]),
00043                           caches, cachew);
00044         DMS::DecodeLatLon(std::string(argv[m + 3]), std::string(argv[m + 4]),
00045                           cachen, cachee);
00046       }
00047       catch (const std::exception& e) {
00048         std::cerr << "Error decoding argument of -c: " << e.what() << "\n";
00049         return 1;
00050       }
00051       m += 4;
00052     } else if (arg == "--msltohae" || arg == "-msltohae")
00053       heightmult = Geoid::GEOIDTOELLIPSOID;
00054     else if (arg == "--haetomsl" || arg == "-haetomsl")
00055       heightmult = Geoid::ELLIPSOIDTOGEOID;
00056     else if (arg == "-z") {
00057       if (++m == argc) return usage(1, true);
00058       zone = argv[m];
00059       bool northp;
00060       int zonenum;
00061       try {
00062         UTMUPS::DecodeZone(zone, zonenum, northp);
00063         zone += ' ';
00064       }
00065       catch (const std::exception& e) {
00066         std::cerr << "Error decoding zone: " << e.what() << "\n";
00067         return 1;
00068       }
00069     } else if (arg == "-n") {
00070       if (++m == argc) return usage(1, true);
00071       geoid = argv[m];
00072     } else if (arg == "-d") {
00073       if (++m == argc) return usage(1, true);
00074       dir = argv[m];
00075     } else if (arg == "-l") {
00076       cubic = false;
00077     } else if (arg == "-v")
00078       verbose = true;
00079     else if (arg == "--version") {
00080       std::cout
00081         << PROGRAM_NAME
00082         << ": $Id: GeoidEval.cpp 6978 2011-02-21 22:42:11Z karney $\n"
00083         << "GeographicLib version " << GEOGRAPHICLIB_VERSION << "\n";
00084       return 0;
00085     } else
00086       return usage(!(arg == "-h" || arg == "--help"), arg != "--help");
00087   }
00088 
00089   int retval = 0;
00090   try {
00091     const Geoid g(geoid, dir, cubic);
00092     try {
00093       if (cacheall)
00094         g.CacheAll();
00095       else if (cachearea)
00096         g.CacheArea(caches, cachew, cachen, cachee);
00097     }
00098     catch (const std::exception& e) {
00099       std::cerr << "ERROR: " << e.what() << "\nProceeding without a cache\n";
00100     }
00101     if (verbose) {
00102       std::cerr << "Geoid file: "    << g.GeoidFile()     << "\n"
00103                 << "Description: "   << g.Description()   << "\n"
00104                 << "Interpolation: " << g.Interpolation() << "\n"
00105                 << "Date & Time: "   << g.DateTime()      << "\n"
00106                 << "Offset (m): "    << g.Offset()        << "\n"
00107                 << "Scale (m): "     << g.Scale()         << "\n"
00108                 << "Max error (m): " << g.MaxError()      << "\n"
00109                 << "RMS error (m): " << g.RMSError()      << "\n";
00110       if (g.Cache())
00111         std::cerr << "Caching:"
00112                   << "\n SW Corner: " << g.CacheSouth() << " " << g.CacheWest()
00113                   << "\n NE Corner: " << g.CacheNorth() << " " << g.CacheEast()
00114                   << "\n";
00115     }
00116 
00117     GeoCoords p;
00118     std::string s;
00119     const char* spaces = " \t\n\v\f\r,"; // Include comma as space
00120     while (std::getline(std::cin, s)) {
00121       try {
00122         real height = 0;
00123         if (heightmult) {
00124           std::string::size_type pb = s.find_last_not_of(spaces);
00125           std::string::size_type pa = s.find_last_of(spaces, pb);
00126           std::string::size_type px = s.find_last_not_of(spaces, pa);
00127           if (pa == std::string::npos || pb == std::string::npos ||
00128               px == std::string::npos)
00129             throw GeographicErr("Incomplete input: " + s);
00130           height = DMS::Decode(s.substr(pa + 1, pb - pa));
00131           s = s.substr(0, px + 1);
00132         }
00133         p.Reset(zone + s);
00134         if (heightmult) {
00135           real h = g(p.Latitude(), p.Longitude());
00136           std::cout << s << " "
00137                     << DMS::Encode(height + real(heightmult) * h,
00138                                    3, DMS::NUMBER) << "\n";
00139         } else {
00140           real gradn, grade;
00141           real h = g(p.Latitude(), p.Longitude(), gradn, grade);
00142           std::cout << DMS::Encode(h, 4, DMS::NUMBER) << " "
00143                     << DMS::Encode(gradn * 1e6, 2, DMS::NUMBER)
00144                     << (Math::isnan(gradn) ? " " : "e-6 ")
00145                     << DMS::Encode(grade * 1e6, 2, DMS::NUMBER)
00146                     << (Math::isnan(grade) ? "\n" : "e-6\n");
00147         }
00148       }
00149       catch (const std::exception& e) {
00150         std::cout << "ERROR: " << e.what() << "\n";
00151         retval = 1;
00152       }
00153     }
00154   }
00155   catch (const std::exception& e) {
00156     std::cerr << "Error reading " << geoid << ": " << e.what() << "\n";
00157     retval = 1;
00158   }
00159   return retval;
00160 }