SUMO - Simulation of Urban MObility
|
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