00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 #include "GeographicLib/Geodesic.hpp"
00017 #include "GeographicLib/AzimuthalEquidistant.hpp"
00018 #include "GeographicLib/CassiniSoldner.hpp"
00019 #include "GeographicLib/DMS.hpp"
00020 #include <iostream>
00021 #include <iomanip>
00022 #include <sstream>
00023
00024 int usage(int retval) {
00025 ( retval ? std::cerr : std::cout ) <<
00026 "Usage: EquidistantTest [-c lat0 lon0] [-z lat0 lon0] [-r] [-h]\n\
00027 $Id: EquidistantTest.cpp 6827 2010-05-20 19:56:18Z karney $\n\
00028 \n\
00029 Convert geodetic coordinates to either azimuthal equidistant or\n\
00030 Cassini-Soldner coordinates. The center of the projection (lat0, lon0)\n\
00031 is specified by either the -c option (for Cassini-Soldner) or -z option\n\
00032 (for azimuthal equidistant). At least one of these options must be\n\
00033 given (the last one given is used). The WGS84 model of the earth is\n\
00034 used.\n\
00035 \n\
00036 Geodetic coordinates are provided on standard input as a set of lines\n\
00037 containing (blank separated) latitude and longitude (degrees or DMS).\n\
00038 For each set of geodetic coordinates, the corresponding projected\n\
00039 coordinates x, y (meters) are printed on standard output together with\n\
00040 the azimuth azi (degrees) and reciprocal scale rk. For Cassini-Soldner,\n\
00041 azi is the bearing of the easting direction and the scale in the\n\
00042 northing direction is 1/rk. For azimuthal equidistant, azi is the\n\
00043 bearing of the radial direction and the scale in the azimuthal direction\n\
00044 is 1/rk.\n\
00045 \n\
00046 If -r is given the reverse transformation is performed. x and y are\n\
00047 given on standard input and each line of the standard output gives\n\
00048 latitude, longitude, azi, and rk.\n\
00049 \n\
00050 -h prints this help.\n";
00051 return retval;
00052 }
00053
00054 int main(int argc, char* argv[]) {
00055 using namespace GeographicLib;
00056 typedef Math::real real;
00057 bool azimuthal = false, cassini = false, reverse = false;
00058 real lat0 = 0, lon0 = 0;
00059 for (int m = 1; m < argc; ++m) {
00060 std::string arg(argv[m]);
00061 if (arg == "-r")
00062 reverse = true;
00063 else if (arg == "-c" || arg == "-z") {
00064 cassini = arg == "-c";
00065 azimuthal = arg != "-c";
00066 if (m + 2 >= argc) return usage(1);
00067 try {
00068 DMS::DecodeLatLon(std::string(argv[m + 1]), std::string(argv[m + 2]),
00069 lat0, lon0);
00070 }
00071 catch (const std::exception& e) {
00072 std::cerr << "Error decoding arguments of " << arg << ": "
00073 << e.what() << "\n";
00074 return 1;
00075 }
00076 m += 2;
00077 } else
00078 return usage(arg != "-h");
00079 }
00080
00081 if (!(azimuthal || cassini)) {
00082 std::cerr << "Must specify \"-c lat lon\" or \"-z lat lon\"\n";
00083 return 1;
00084 }
00085
00086 const CassiniSoldner cs = cassini ?
00087 CassiniSoldner(lat0, lon0, Geodesic::WGS84) :
00088 CassiniSoldner(Geodesic::WGS84);
00089 const AzimuthalEquidistant az(Geodesic::WGS84);
00090
00091 std::string s;
00092 int retval = 0;
00093 std::cout << std::setprecision(16);
00094 while (std::getline(std::cin, s)) {
00095 try {
00096 std::istringstream str(s);
00097 real lat, lon, x, y, a, m;
00098 std::string stra, strb;
00099 if (!(str >> stra >> strb))
00100 throw GeographicErr("Incomplete input: " + s);
00101 if (reverse) {
00102 x = DMS::Decode(stra);
00103 y = DMS::Decode(strb);
00104 } else
00105 DMS::DecodeLatLon(stra, strb, lat, lon);
00106 std::string strc;
00107 if (str >> strc)
00108 throw GeographicErr("Extraneous input: " + strc);
00109 if (reverse) {
00110 if (cassini)
00111 cs.Reverse(x, y, lat, lon, a, m);
00112 else
00113 az.Reverse(lat0, lon0, x, y, lat, lon, a, m);
00114 std::cout << lat << " " << lon << " " << a << " " << m << "\n";
00115 } else {
00116 if ( !(-90 <= lat && lat <= 90) )
00117 throw GeographicErr("Latitude not in range [-90, 90]");
00118 if ( !(-180 <= lon && lat <= 360) )
00119 throw GeographicErr("Longitude not in range [-180, 360]");
00120 if (cassini)
00121 cs.Forward(lat, lon, x, y, a, m);
00122 else
00123 az.Forward(lat0, lon0, lat, lon, x, y, a, m);
00124 std::cout << x << " " << y << " " << a << " " << m << "\n";
00125 }
00126 }
00127 catch (const std::exception& e) {
00128 std::cout << "ERROR: " << e.what() << "\n";
00129 retval = 1;
00130 }
00131 }
00132
00133 return retval;
00134 }