SUMO - Simulation of Urban MObility
GeoConvHelper.cpp
Go to the documentation of this file.
00001 /****************************************************************************/
00009 // static methods for processing the coordinates conversion for the current net
00010 /****************************************************************************/
00011 // SUMO, Simulation of Urban MObility; see http://sumo.sourceforge.net/
00012 // Copyright (C) 2001-2012 DLR (http://www.dlr.de/) and contributors
00013 /****************************************************************************/
00014 //
00015 //   This file is part of SUMO.
00016 //   SUMO is free software: you can redistribute it and/or modify
00017 //   it under the terms of the GNU General Public License as published by
00018 //   the Free Software Foundation, either version 3 of the License, or
00019 //   (at your option) any later version.
00020 //
00021 /****************************************************************************/
00022 
00023 
00024 // ===========================================================================
00025 // included modules
00026 // ===========================================================================
00027 #ifdef _MSC_VER
00028 #include <windows_config.h>
00029 #else
00030 #include <config.h>
00031 #endif
00032 
00033 #include <map>
00034 #include <cmath>
00035 #include <cassert>
00036 #include <utils/common/MsgHandler.h>
00037 #include <utils/common/ToString.h>
00038 #include <utils/geom/GeomHelper.h>
00039 #include <utils/options/OptionsCont.h>
00040 #include "GeoConvHelper.h"
00041 
00042 #ifdef CHECK_MEMORY_LEAKS
00043 #include <foreign/nvwa/debug_new.h>
00044 #endif // CHECK_MEMORY_LEAKS
00045 
00046 
00047 // ===========================================================================
00048 // static member variables
00049 // ===========================================================================
00050 GeoConvHelper GeoConvHelper::myProcessing("!", Position(), Boundary(), Boundary());
00051 GeoConvHelper GeoConvHelper::myLoaded("!", Position(), Boundary(), Boundary());
00052 GeoConvHelper GeoConvHelper::myFinal("!", Position(), Boundary(), Boundary());
00053 int GeoConvHelper::myNumLoaded = 0;
00054 
00055 // ===========================================================================
00056 // method definitions
00057 // ===========================================================================
00058 GeoConvHelper::GeoConvHelper(const std::string& proj, const Position& offset,
00059                              const Boundary& orig, const Boundary& conv, int shift, bool inverse, bool baseFound):
00060     myProjString(proj),
00061 #ifdef HAVE_PROJ
00062     myProjection(0),
00063 #endif
00064     myOffset(offset),
00065     myProjectionMethod(NONE),
00066     myOrigBoundary(orig),
00067     myConvBoundary(conv),
00068     myGeoScale(pow(10, (double) - shift)),
00069     myUseInverseProjection(inverse),
00070     myBaseFound(baseFound),
00071     myBaseX(0),
00072     myBaseY(0) {
00073     if (proj == "!") {
00074         myProjectionMethod = NONE;
00075     } else if (proj == "-") {
00076         myProjectionMethod = SIMPLE;
00077     } else if (proj == "UTM") {
00078         myProjectionMethod = UTM;
00079     } else if (proj == "DHDN") {
00080         myProjectionMethod = DHDN;
00081 #ifdef HAVE_PROJ
00082     } else {
00083         myProjectionMethod = PROJ;
00084         myProjection = pj_init_plus(proj.c_str());
00085         if (myProjection == 0) {
00086             // !!! check pj_errno
00087             throw ProcessError("Could not build projection!");
00088         }
00089 #endif
00090     }
00091 }
00092 
00093 
00094 GeoConvHelper::~GeoConvHelper() {
00095 #ifdef HAVE_PROJ
00096     if (myProjection != 0) {
00097         pj_free(myProjection);
00098     }
00099 #endif
00100 }
00101 
00102 
00103 GeoConvHelper&
00104 GeoConvHelper::operator=(const GeoConvHelper& orig) {
00105     myProjString = orig.myProjString;
00106     myOffset = orig.myOffset;
00107     myProjectionMethod = orig.myProjectionMethod;
00108     myOrigBoundary = orig.myOrigBoundary;
00109     myConvBoundary = orig.myConvBoundary;
00110     myGeoScale = orig.myGeoScale;
00111     myUseInverseProjection = orig.myUseInverseProjection;
00112     myBaseFound = orig.myBaseFound;
00113     myBaseX = orig.myBaseX;
00114     myBaseY = orig.myBaseY;
00115 #ifdef HAVE_PROJ
00116     if (myProjection != 0) {
00117         pj_free(myProjection);
00118         myProjection = 0;
00119     }
00120     if (orig.myProjection != 0) {
00121         myProjection = pj_init_plus(orig.myProjString.c_str());
00122     }
00123 #endif
00124     return *this;
00125 }
00126 
00127 
00128 bool
00129 GeoConvHelper::init(OptionsCont& oc) {
00130     std::string proj = "!"; // the default
00131     int shift = oc.getInt("proj.scale");
00132     Position offset = Position(oc.getFloat("offset.x"), oc.getFloat("offset.y"));
00133     bool inverse = oc.exists("proj.inverse") && oc.getBool("proj.inverse");
00134     bool baseFound = !oc.exists("offset.disable-normalization") ||
00135                      oc.getBool("offset.disable-normalization") ||
00136                      !oc.isDefault("offset.x") ||
00137                      !oc.isDefault("offset.y");
00138 
00139     if (oc.getBool("simple-projection")) {
00140         proj = "-";
00141     }
00142 
00143 #ifdef HAVE_PROJ
00144     if (oc.getBool("proj.inverse") && oc.getString("proj") == "!") {
00145         WRITE_ERROR("Inverse projection works only with explicit proj parameters.");
00146         return false;
00147     }
00148     unsigned numProjections = oc.getBool("simple-projection") + oc.getBool("proj.utm") + oc.getBool("proj.dhdn") + (oc.getString("proj").length() > 1);
00149     if (numProjections > 1) {
00150         WRITE_ERROR("The projection method needs to be uniquely defined.");
00151         return false;
00152     }
00153 
00154     if (oc.getBool("proj.utm")) {
00155         proj = "UTM";
00156     } else if (oc.getBool("proj.dhdn")) {
00157         proj = "DHDN";
00158     } else {
00159         proj = oc.getString("proj");
00160     }
00161 #endif
00162     myProcessing = GeoConvHelper(proj, offset, Boundary(), Boundary(), shift, inverse, baseFound);
00163     myFinal = myProcessing;
00164     return true;
00165 }
00166 
00167 
00168 void
00169 GeoConvHelper::init(const std::string& proj,
00170                     const Position& offset,
00171                     const Boundary& orig,
00172                     const Boundary& conv) {
00173     myProcessing = GeoConvHelper(proj, offset, orig, conv);
00174     myFinal = myProcessing;
00175 }
00176 
00177 
00178 void
00179 GeoConvHelper::addProjectionOptions(OptionsCont& oc) {
00180     oc.addOptionSubTopic("Projection");
00181 
00182     oc.doRegister("simple-projection", new Option_Bool(false));
00183     oc.addSynonyme("simple-projection", "proj.simple", true);
00184     oc.addDescription("simple-projection", "Projection", "Uses a simple method for projection");
00185 
00186     oc.doRegister("proj.scale", new Option_Integer(0));
00187     oc.addSynonyme("proj.scale", "proj.shift", true);
00188     oc.addDescription("proj.scale", "Projection", "Number of places to shift decimal point to right in geo-coordinates");
00189 
00190 #ifdef HAVE_PROJ
00191     oc.doRegister("proj.utm", new Option_Bool(false));
00192     oc.addDescription("proj.utm", "Projection", "Determine the UTM zone (for a universal transversal mercator projection based on the WGS84 ellipsoid)");
00193 
00194     oc.doRegister("proj.dhdn", new Option_Bool(false));
00195     oc.addDescription("proj.dhdn", "Projection", "Determine the DHDN zone (for a transversal mercator projection based on the bessel ellipsoid)");
00196 
00197     oc.doRegister("proj", new Option_String("!"));
00198     oc.addDescription("proj", "Projection", "Uses STR as proj.4 definition for projection");
00199 
00200     oc.doRegister("proj.inverse", new Option_Bool(false));
00201     oc.addDescription("proj.inverse", "Projection", "Inverses projection");
00202 #endif // HAVE_PROJ
00203 }
00204 
00205 
00206 bool
00207 GeoConvHelper::usingGeoProjection() const {
00208     return myProjectionMethod != NONE;
00209 }
00210 
00211 
00212 bool
00213 GeoConvHelper::usingInverseGeoProjection() const {
00214     return myUseInverseProjection;
00215 }
00216 
00217 
00218 void
00219 GeoConvHelper::cartesian2geo(Position& cartesian) const {
00220     cartesian.sub(getOffsetBase());
00221     if (myProjectionMethod == NONE) {
00222         return;
00223     }
00224 #ifdef HAVE_PROJ
00225     projUV p;
00226     p.u = cartesian.x();
00227     p.v = cartesian.y();
00228     p = pj_inv(p, myProjection);
00230     p.u *= RAD_TO_DEG;
00231     p.v *= RAD_TO_DEG;
00232     cartesian.set((SUMOReal) p.u, (SUMOReal) p.v);
00233 #endif
00234 }
00235 
00236 
00237 bool
00238 GeoConvHelper::x2cartesian(Position& from, bool includeInBoundary) {
00239     myOrigBoundary.add(from);
00240     // init projection parameter on first use
00241 #ifdef HAVE_PROJ
00242     if (myProjection == 0) {
00243         const double x = from.x() * myGeoScale;
00244         switch(myProjectionMethod) {
00245             case UTM: {
00246                 int zone = (int)(x + 180) / 6 + 1;
00247                 myProjString = "+proj=utm +zone=" + toString(zone) +
00248                     " +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
00249                 myProjection = pj_init_plus(myProjString.c_str());
00251             }
00252             break;
00253             case DHDN: {
00254                 int zone = (int)(x / 3);
00255                 if (zone < 1 || zone > 5) {
00256                     WRITE_WARNING("Attempt to initialize DHDN-projection on invalid longitude " + toString(x));
00257                     return false;
00258                 }
00259                 myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
00260                     " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
00261                     " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
00262                 myProjection = pj_init_plus(myProjString.c_str());
00264             } 
00265             break;
00266             default:
00267             break;
00268         }
00269     }
00270 #endif
00271     // perform conversion
00272     bool ok = x2cartesian_const(from);
00273     if (ok) {
00274         if (!myBaseFound) {
00275             // avoid very large coordinates to reduce numerical errors
00276             if (myProjectionMethod == SIMPLE || from.x() > 100000 || from.y() > 100000) {
00277                 myBaseX = from.x();
00278                 myBaseY = from.y();
00279                 from.set(0,0);
00280             }
00281             myBaseFound = true;
00282         }
00283 
00284         if (includeInBoundary) {
00285             myConvBoundary.add(from);
00286         }
00287     }
00288     return ok;
00289 }
00290 
00291 
00292 bool
00293 GeoConvHelper::x2cartesian_const(Position& from) const {
00294     double x = from.x() * myGeoScale;
00295     double y = from.y() * myGeoScale;
00296     if (myProjectionMethod == NONE) {
00297         from.add(myOffset);
00298     } else if (myUseInverseProjection) {
00299         cartesian2geo(from);
00300     } else {
00301         if (x > 180.1 || x < -180.1 || y > 90.1 || y < -90.1) {
00302             return false;
00303         }
00304 #ifdef HAVE_PROJ
00305         if (myProjection != 0) {
00306             projUV p;
00307             p.u = x * DEG_TO_RAD;
00308             p.v = y * DEG_TO_RAD;
00309             p = pj_fwd(p, myProjection);
00311             x = p.u;
00312             y = p.v;
00313         }
00314 #endif
00315         if (myProjectionMethod == SIMPLE) {
00316             double ys = y;
00317             x -= myBaseX;
00318             y -= myBaseY;
00319             x *= 111320. * cos(ys * PI / 180.0);
00320             y *= 111136.;
00321             from.set((SUMOReal)x, (SUMOReal)y);
00323             from.add(myOffset);
00324         }
00325     }
00326     if (myProjectionMethod != SIMPLE) {
00327         x -= myBaseX;
00328         y -= myBaseY;
00329         from.set((SUMOReal)x, (SUMOReal)y);
00330         from.add(myOffset);
00331     }
00332     return true;
00333 }
00334 
00335 
00336 void
00337 GeoConvHelper::moveConvertedBy(SUMOReal x, SUMOReal y) {
00338     myOffset.add(x, y);
00339     myConvBoundary.moveby(x, y);
00340 }
00341 
00342 
00343 const Boundary&
00344 GeoConvHelper::getOrigBoundary() const {
00345     return myOrigBoundary;
00346 }
00347 
00348 
00349 const Boundary&
00350 GeoConvHelper::getConvBoundary() const {
00351     return myConvBoundary;
00352 }
00353 
00354 
00355 const Position
00356 GeoConvHelper::getOffset() const {
00357     return myOffset;
00358 }
00359 
00360 
00361 const Position
00362 GeoConvHelper::getOffsetBase() const {
00363     return Position(myOffset.x() - myBaseX, myOffset.y() - myBaseY);
00364 }
00365 
00366 
00367 const std::string&
00368 GeoConvHelper::getProjString() const {
00369     return myProjString;
00370 }
00371 
00372 
00373 void 
00374 GeoConvHelper::computeFinal() {
00375     if (myNumLoaded == 0) {
00376         myFinal = myProcessing;
00377     } else  {
00378         myFinal = GeoConvHelper(
00379                 // prefer options over loaded location
00380                 myProcessing.usingGeoProjection() ? myProcessing.getProjString() : myLoaded.getProjString(),
00381                 // let offset and boundary lead back to the original coords of the loaded data
00382                 myProcessing.getOffset() + myLoaded.getOffset(),
00383                 myLoaded.getOrigBoundary(),
00384                 // the new boundary (updated during loading)
00385                 myProcessing.getConvBoundary());
00386     }
00387 }
00388 
00389 
00390 void 
00391 GeoConvHelper::setLoaded(const GeoConvHelper& loaded) {
00392     myNumLoaded++;
00393     if (myNumLoaded > 1) {
00394         WRITE_WARNING("Ignoring loaded location attribute nr. " + toString(myNumLoaded) + " for tracking of original location");
00395     } else {
00396         myLoaded = loaded;
00397     }
00398 }
00399 
00400 
00401 void 
00402 GeoConvHelper::resetLoaded() {
00403     myNumLoaded = 0;
00404 }
00405 
00406 
00407 /****************************************************************************/
00408 
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Defines