GeographicLib  1.21
GeoidEval.cpp
Go to the documentation of this file.
00001 /**
00002  * \file GeoidEval.cpp
00003  * \brief Command line utility for evaluating geoid heights
00004  *
00005  * Copyright (c) Charles Karney (2009-2012) <charles@karney.com> and licensed
00006  * under the MIT/X11 License.  For more information, see
00007  * http://geographiclib.sourceforge.net/
00008  *
00009  * Compile and link with
00010  *   g++ -g -O3 -I../include -I../man -o GeoidEval \
00011  *       GeoidEval.cpp \
00012  *       ../src/DMS.cpp \
00013  *       ../src/GeoCoords.cpp \
00014  *       ../src/Geoid.cpp \
00015  *       ../src/MGRS.cpp \
00016  *       ../src/PolarStereographic.cpp \
00017  *       ../src/TransverseMercator.cpp \
00018  *       ../src/UTMUPS.cpp
00019  *
00020  * See the <a href="GeoidEval.1.html">man page</a> for usage
00021  * information.
00022  **********************************************************************/
00023 
00024 #include <iostream>
00025 #include <string>
00026 #include <sstream>
00027 #include <fstream>
00028 #include <GeographicLib/Geoid.hpp>
00029 #include <GeographicLib/DMS.hpp>
00030 #include <GeographicLib/Utility.hpp>
00031 #include <GeographicLib/GeoCoords.hpp>
00032 
00033 #include "GeoidEval.usage"
00034 
00035 int main(int argc, char* argv[]) {
00036   try {
00037     using namespace GeographicLib;
00038     typedef Math::real real;
00039     bool cacheall = false, cachearea = false, verbose = false,
00040       cubic = true, gradp = false;
00041     real caches, cachew, cachen, cachee;
00042     std::string dir;
00043     std::string geoid = Geoid::DefaultGeoidName();
00044     Geoid::convertflag heightmult = Geoid::NONE;
00045     std::string istring, ifile, ofile, cdelim;
00046     char lsep = ';';
00047     bool northp = false;
00048     int zonenum = UTMUPS::INVALID;
00049 
00050     for (int m = 1; m < argc; ++m) {
00051       std::string arg(argv[m]);
00052       if (arg == "-a") {
00053         cacheall = true;
00054         cachearea = false;
00055       }
00056       else if (arg == "-c") {
00057         if (m + 4 >= argc) return usage(1, true);
00058         cacheall = false;
00059         cachearea = true;
00060         try {
00061           DMS::DecodeLatLon(std::string(argv[m + 1]), std::string(argv[m + 2]),
00062                             caches, cachew);
00063           DMS::DecodeLatLon(std::string(argv[m + 3]), std::string(argv[m + 4]),
00064                             cachen, cachee);
00065         }
00066         catch (const std::exception& e) {
00067           std::cerr << "Error decoding argument of -c: " << e.what() << "\n";
00068           return 1;
00069         }
00070         m += 4;
00071       } else if (arg == "--msltohae")
00072         heightmult = Geoid::GEOIDTOELLIPSOID;
00073       else if (arg == "--haetomsl")
00074         heightmult = Geoid::ELLIPSOIDTOGEOID;
00075       else if (arg == "-z") {
00076         if (++m == argc) return usage(1, true);
00077         std::string zone = argv[m];
00078         try {
00079           UTMUPS::DecodeZone(zone, zonenum, northp);
00080         }
00081         catch (const std::exception& e) {
00082           std::cerr << "Error decoding zone: " << e.what() << "\n";
00083           return 1;
00084         }
00085         if (!(zonenum >= UTMUPS::MINZONE && zonenum <= UTMUPS::MAXZONE)) {
00086           std::cerr << "Illegal zone " << zone << "\n";
00087           return 1;
00088         }
00089       } else if (arg == "-n") {
00090         if (++m == argc) return usage(1, true);
00091         geoid = argv[m];
00092       } else if (arg == "-d") {
00093         if (++m == argc) return usage(1, true);
00094         dir = argv[m];
00095       } else if (arg == "-l")
00096         cubic = false;
00097       else if (arg == "-g")
00098         gradp = true;
00099       else if (arg == "-v")
00100         verbose = true;
00101       else if (arg == "--input-string") {
00102         if (++m == argc) return usage(1, true);
00103         istring = argv[m];
00104       } else if (arg == "--input-file") {
00105         if (++m == argc) return usage(1, true);
00106         ifile = argv[m];
00107       } else if (arg == "--output-file") {
00108         if (++m == argc) return usage(1, true);
00109         ofile = argv[m];
00110       } else if (arg == "--line-separator") {
00111         if (++m == argc) return usage(1, true);
00112         if (std::string(argv[m]).size() != 1) {
00113           std::cerr << "Line separator must be a single character\n";
00114           return 1;
00115         }
00116         lsep = argv[m][0];
00117       } else if (arg == "--comment-delimiter") {
00118         if (++m == argc) return usage(1, true);
00119         cdelim = argv[m];
00120       } else if (arg == "--version") {
00121         std::cout
00122           << argv[0]
00123           << ": $Id: 6db1ff0b8309a39d0d9b0250dd73be964c5efb7c $\n"
00124           << "GeographicLib version " << GEOGRAPHICLIB_VERSION_STRING << "\n";
00125         return 0;
00126       } else {
00127         int retval = usage(!(arg == "-h" || arg == "--help"), arg != "--help");
00128         if (arg == "-h")
00129           std::cout<< "\nDefault geoid path = \""   << Geoid::DefaultGeoidPath()
00130                    << "\"\nDefault geoid name = \"" << Geoid::DefaultGeoidName()
00131                    << "\"\n";
00132         return retval;
00133       }
00134     }
00135 
00136     if (!ifile.empty() && !istring.empty()) {
00137       std::cerr << "Cannot specify --input-string and --input-file together\n";
00138       return 1;
00139     }
00140     if (ifile == "-") ifile.clear();
00141     std::ifstream infile;
00142     std::istringstream instring;
00143     if (!ifile.empty()) {
00144       infile.open(ifile.c_str());
00145       if (!infile.is_open()) {
00146         std::cerr << "Cannot open " << ifile << " for reading\n";
00147         return 1;
00148       }
00149     } else if (!istring.empty()) {
00150       std::string::size_type m = 0;
00151       while (true) {
00152         m = istring.find(lsep, m);
00153         if (m == std::string::npos)
00154           break;
00155         istring[m] = '\n';
00156       }
00157       instring.str(istring);
00158     }
00159     std::istream* input = !ifile.empty() ? &infile :
00160       (!istring.empty() ? &instring : &std::cin);
00161 
00162     std::ofstream outfile;
00163     if (ofile == "-") ofile.clear();
00164     if (!ofile.empty()) {
00165       outfile.open(ofile.c_str());
00166       if (!outfile.is_open()) {
00167         std::cerr << "Cannot open " << ofile << " for writing\n";
00168         return 1;
00169       }
00170     }
00171     std::ostream* output = !ofile.empty() ? &outfile : &std::cout;
00172 
00173     int retval = 0;
00174     try {
00175       const Geoid g(geoid, dir, cubic);
00176       try {
00177         if (cacheall)
00178           g.CacheAll();
00179         else if (cachearea)
00180           g.CacheArea(caches, cachew, cachen, cachee);
00181       }
00182       catch (const std::exception& e) {
00183         std::cerr << "ERROR: " << e.what() << "\nProceeding without a cache\n";
00184       }
00185       if (verbose) {
00186         std::cerr << "Geoid file: "    << g.GeoidFile()     << "\n"
00187                   << "Description: "   << g.Description()   << "\n"
00188                   << "Interpolation: " << g.Interpolation() << "\n"
00189                   << "Date & Time: "   << g.DateTime()      << "\n"
00190                   << "Offset (m): "    << g.Offset()        << "\n"
00191                   << "Scale (m): "     << g.Scale()         << "\n"
00192                   << "Max error (m): " << g.MaxError()      << "\n"
00193                   << "RMS error (m): " << g.RMSError()      << "\n";
00194         if (g.Cache())
00195           std::cerr<< "Caching:"
00196                    << "\n SW Corner: " << g.CacheSouth() << " " << g.CacheWest()
00197                    << "\n NE Corner: " << g.CacheNorth() << " " << g.CacheEast()
00198                    << "\n";
00199       }
00200 
00201       GeoCoords p;
00202       std::string s, suff;
00203       const char* spaces = " \t\n\v\f\r,"; // Include comma as space
00204       while (std::getline(*input, s)) {
00205         try {
00206           std::string eol("\n");
00207           if (!cdelim.empty()) {
00208             std::string::size_type m = s.find(cdelim);
00209             if (m != std::string::npos) {
00210               eol = " " + s.substr(m) + "\n";
00211               std::string::size_type m1 =
00212                 m > 0 ? s.find_last_not_of(spaces, m - 1) : std::string::npos;
00213               s = s.substr(0, m1 != std::string::npos ? m1 + 1 : m);
00214             }
00215           }
00216           real height = 0;
00217           if (zonenum != UTMUPS::INVALID) {
00218             // Expect "easting northing" if heightmult == 0, or
00219             // "easting northing height" if heightmult != 0.
00220             std::string::size_type pa = 0, pb = 0;
00221             real easting = 0, northing = 0;
00222             for (int i = 0; i < (heightmult ? 3 : 2); ++i) {
00223               if (pb == std::string::npos)
00224                 throw GeographicErr("Incomplete input: " + s);
00225               // Start of i'th token
00226               pa = s.find_first_not_of(spaces, pb);
00227               if (pa == std::string::npos)
00228                 throw GeographicErr("Incomplete input: " + s);
00229               // End of i'th token
00230               pb = s.find_first_of(spaces, pa);
00231               (i == 2 ? height : (i == 0 ? easting : northing)) =
00232                 Utility::num<real>(s.substr(pa, (pb == std::string::npos ?
00233                                                  pb : pb - pa)));
00234             }
00235             p.Reset(zonenum, northp, easting, northing);
00236             if (heightmult) {
00237               suff = pb == std::string::npos ? "" : s.substr(pb);
00238               s = s.substr(0, pa);
00239             }
00240           } else {
00241             if (heightmult) {
00242               // Treat last token as height
00243               // pb = last char of last token
00244               // pa = last char preceding white space
00245               // px = last char of 2nd last token
00246               std::string::size_type pb = s.find_last_not_of(spaces);
00247               std::string::size_type pa = s.find_last_of(spaces, pb);
00248               if (pa == std::string::npos || pb == std::string::npos)
00249                 throw GeographicErr("Incomplete input: " + s);
00250               height = Utility::num<real>(s.substr(pa + 1, pb - pa));
00251               s = s.substr(0, pa + 1);
00252             }
00253             p.Reset(s);
00254           }
00255           if (heightmult) {
00256             real h = g(p.Latitude(), p.Longitude());
00257             *output << s
00258                     << Utility::str<real>(height + real(heightmult) * h, 4)
00259                     << suff << eol;
00260           } else {
00261             if (gradp) {
00262             real gradn, grade;
00263             real h = g(p.Latitude(), p.Longitude(), gradn, grade);
00264             *output << Utility::str<real>(h, 4) << " "
00265                     << Utility::str<real>(gradn * 1e6, 2)
00266                     << (Math::isnan(gradn) ? " " : "e-6 ")
00267                     << Utility::str<real>(grade * 1e6, 2)
00268                     << (Math::isnan(grade) ? "" : "e-6")
00269                     << eol;
00270             } else {
00271             real h = g(p.Latitude(), p.Longitude());
00272             *output << Utility::str<real>(h, 4) << eol;
00273             }
00274           }
00275         }
00276         catch (const std::exception& e) {
00277           *output << "ERROR: " << e.what() << "\n";
00278           retval = 1;
00279         }
00280       }
00281     }
00282     catch (const std::exception& e) {
00283       std::cerr << "Error reading " << geoid << ": " << e.what() << "\n";
00284       retval = 1;
00285     }
00286     return retval;
00287   }
00288   catch (const std::exception& e) {
00289     std::cerr << "Caught exception: " << e.what() << "\n";
00290     return 1;
00291   }
00292   catch (...) {
00293     std::cerr << "Caught unknown exception\n";
00294     return 1;
00295   }
00296 }