SUMO - Simulation of Urban MObility
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GeoConvHelper.cpp
Go to the documentation of this file.
1 /****************************************************************************/
9 // static methods for processing the coordinates conversion for the current net
10 /****************************************************************************/
11 // SUMO, Simulation of Urban MObility; see http://sumo-sim.org/
12 // Copyright (C) 2001-2014 DLR (http://www.dlr.de/) and contributors
13 /****************************************************************************/
14 //
15 // This file is part of SUMO.
16 // SUMO is free software: you can redistribute it and/or modify
17 // it under the terms of the GNU General Public License as published by
18 // the Free Software Foundation, either version 3 of the License, or
19 // (at your option) any later version.
20 //
21 /****************************************************************************/
22 
23 
24 // ===========================================================================
25 // included modules
26 // ===========================================================================
27 #ifdef _MSC_VER
28 #include <windows_config.h>
29 #else
30 #include <config.h>
31 #endif
32 
33 #include <map>
34 #include <cmath>
35 #include <cassert>
36 #include <climits>
38 #include <utils/common/ToString.h>
39 #include <utils/geom/GeomHelper.h>
41 #include "GeoConvHelper.h"
42 
43 #ifdef CHECK_MEMORY_LEAKS
44 #include <foreign/nvwa/debug_new.h>
45 #endif // CHECK_MEMORY_LEAKS
46 
47 
48 // ===========================================================================
49 // static member variables
50 // ===========================================================================
55 
56 // ===========================================================================
57 // method definitions
58 // ===========================================================================
59 GeoConvHelper::GeoConvHelper(const std::string& proj, const Position& offset,
60  const Boundary& orig, const Boundary& conv, int shift, bool inverse):
61  myProjString(proj),
62 #ifdef HAVE_PROJ
63  myProjection(0),
64  myInverseProjection(0),
65  myGeoProjection(0),
66 #endif
67  myOffset(offset),
68  myGeoScale(pow(10, (double) - shift)),
69  myProjectionMethod(NONE),
70  myUseInverseProjection(inverse),
71  myOrigBoundary(orig),
72  myConvBoundary(conv) {
73  if (proj == "!") {
75  } else if (proj == "-") {
77  } else if (proj == "UTM") {
79  } else if (proj == "DHDN") {
81  } else if (proj == "DHDN_UTM") {
83 #ifdef HAVE_PROJ
84  } else {
86  myProjection = pj_init_plus(proj.c_str());
87  if (myProjection == 0) {
88  // !!! check pj_errno
89  throw ProcessError("Could not build projection!");
90  }
91 #endif
92  }
93 }
94 
95 
97 #ifdef HAVE_PROJ
98  if (myProjection != 0) {
99  pj_free(myProjection);
100  }
101  if (myInverseProjection != 0) {
102  pj_free(myInverseProjection);
103  }
104  if (myGeoProjection != 0) {
105  pj_free(myInverseProjection);
106  }
107 #endif
108 }
109 
110 
113  myProjString = orig.myProjString;
114  myOffset = orig.myOffset;
118  myGeoScale = orig.myGeoScale;
120 #ifdef HAVE_PROJ
121  if (myProjection != 0) {
122  pj_free(myProjection);
123  myProjection = 0;
124  }
125  if (myInverseProjection != 0) {
126  pj_free(myInverseProjection);
127  myInverseProjection = 0;
128  }
129  if (myGeoProjection != 0) {
130  pj_free(myGeoProjection);
131  myGeoProjection = 0;
132  }
133  if (orig.myProjection != 0) {
134  myProjection = pj_init_plus(orig.myProjString.c_str());
135  }
136  if (orig.myInverseProjection != 0) {
137  myInverseProjection = pj_init_plus(pj_get_def(orig.myInverseProjection, 0));
138  }
139  if (orig.myGeoProjection != 0) {
140  myGeoProjection = pj_init_plus(pj_get_def(orig.myGeoProjection, 0));
141  }
142 #endif
143  return *this;
144 }
145 
146 
147 bool
149  std::string proj = "!"; // the default
150  int shift = oc.getInt("proj.scale");
151  Position offset = Position(oc.getFloat("offset.x"), oc.getFloat("offset.y"));
152  bool inverse = oc.exists("proj.inverse") && oc.getBool("proj.inverse");
153 
154  if (oc.getBool("simple-projection")) {
155  proj = "-";
156  }
157 
158 #ifdef HAVE_PROJ
159  if (oc.getBool("proj.inverse") && oc.getString("proj") == "!") {
160  WRITE_ERROR("Inverse projection works only with explicit proj parameters.");
161  return false;
162  }
163  unsigned numProjections = oc.getBool("simple-projection") + oc.getBool("proj.utm") + oc.getBool("proj.dhdn") + oc.getBool("proj.dhdnutm") + (oc.getString("proj").length() > 1);
164  if (numProjections > 1) {
165  WRITE_ERROR("The projection method needs to be uniquely defined.");
166  return false;
167  }
168 
169  if (oc.getBool("proj.utm")) {
170  proj = "UTM";
171  } else if (oc.getBool("proj.dhdn")) {
172  proj = "DHDN";
173  } else if (oc.getBool("proj.dhdnutm")) {
174  proj = "DHDN_UTM";
175  } else {
176  proj = oc.getString("proj");
177  }
178 #endif
179  myProcessing = GeoConvHelper(proj, offset, Boundary(), Boundary(), shift, inverse);
181  return true;
182 }
183 
184 
185 void
186 GeoConvHelper::init(const std::string& proj,
187  const Position& offset,
188  const Boundary& orig,
189  const Boundary& conv,
190  int shift) {
191  myProcessing = GeoConvHelper(proj, offset, orig, conv, shift);
193 }
194 
195 
196 void
198  oc.addOptionSubTopic("Projection");
199 
200  oc.doRegister("simple-projection", new Option_Bool(false));
201  oc.addSynonyme("simple-projection", "proj.simple", true);
202  oc.addDescription("simple-projection", "Projection", "Uses a simple method for projection");
203 
204  oc.doRegister("proj.scale", new Option_Integer(0));
205  oc.addSynonyme("proj.scale", "proj.shift", true);
206  oc.addDescription("proj.scale", "Projection", "Number of places to shift decimal point to right in geo-coordinates");
207 
208 #ifdef HAVE_PROJ
209  oc.doRegister("proj.utm", new Option_Bool(false));
210  oc.addDescription("proj.utm", "Projection", "Determine the UTM zone (for a universal transversal mercator projection based on the WGS84 ellipsoid)");
211 
212  oc.doRegister("proj.dhdn", new Option_Bool(false));
213  oc.addDescription("proj.dhdn", "Projection", "Determine the DHDN zone (for a transversal mercator projection based on the bessel ellipsoid, \"Gauss-Krueger\")");
214 
215  oc.doRegister("proj", new Option_String("!"));
216  oc.addDescription("proj", "Projection", "Uses STR as proj.4 definition for projection");
217 
218  oc.doRegister("proj.inverse", new Option_Bool(false));
219  oc.addDescription("proj.inverse", "Projection", "Inverses projection");
220 
221  oc.doRegister("proj.dhdnutm", new Option_Bool(false));
222  oc.addDescription("proj.dhdnutm", "Projection", "Convert from Gauss-Krueger to UTM");
223 #endif // HAVE_PROJ
224 }
225 
226 
227 bool
229  return myProjectionMethod != NONE;
230 }
231 
232 
233 bool
235  return myUseInverseProjection;
236 }
237 
238 
239 void
241  cartesian.sub(getOffsetBase());
242  if (myProjectionMethod == NONE) {
243  return;
244  }
245 #ifdef HAVE_PROJ
246  projUV p;
247  p.u = cartesian.x();
248  p.v = cartesian.y();
249  p = pj_inv(p, myProjection);
251  p.u *= RAD_TO_DEG;
252  p.v *= RAD_TO_DEG;
253  cartesian.set((SUMOReal) p.u, (SUMOReal) p.v);
254 #endif
255 }
256 
257 
258 bool
259 GeoConvHelper::x2cartesian(Position& from, bool includeInBoundary) {
260  if (includeInBoundary) {
261  myOrigBoundary.add(from);
262  }
263  // init projection parameter on first use
264 #ifdef HAVE_PROJ
265  if (myProjection == 0) {
266  double x = from.x() * myGeoScale;
267  switch (myProjectionMethod) {
268  case DHDN_UTM: {
269  int zone = (int)((x - 500000.) / 1000000.);
270  if (zone < 1 || zone > 5) {
271  WRITE_WARNING("Attempt to initialize DHDN_UTM-projection on invalid longitude " + toString(x));
272  return false;
273  }
274  myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
275  " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
276  " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
277  myInverseProjection = pj_init_plus(myProjString.c_str());
278  myGeoProjection = pj_init_plus("+proj=latlong +datum=WGS84");
280  x = ((x - 500000.) / 1000000.) * 3; // continues with UTM
281  }
282  case UTM: {
283  int zone = (int)(x + 180) / 6 + 1;
284  myProjString = "+proj=utm +zone=" + toString(zone) +
285  " +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
286  myProjection = pj_init_plus(myProjString.c_str());
288  }
289  break;
290  case DHDN: {
291  int zone = (int)(x / 3);
292  if (zone < 1 || zone > 5) {
293  WRITE_WARNING("Attempt to initialize DHDN-projection on invalid longitude " + toString(x));
294  return false;
295  }
296  myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
297  " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
298  " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
299  myProjection = pj_init_plus(myProjString.c_str());
301  }
302  break;
303  default:
304  break;
305  }
306  }
307  if (myInverseProjection != 0) {
308  double x = from.x();
309  double y = from.y();
310  if (pj_transform(myInverseProjection, myGeoProjection, 1, 1, &x, &y, NULL)) {
311  WRITE_WARNING("Could not transform (" + toString(x) + "," + toString(y) + ")");
312  }
313  from.set(SUMOReal(x * RAD_TO_DEG), SUMOReal(y * RAD_TO_DEG));
314  }
315 #endif
316  // perform conversion
317  bool ok = x2cartesian_const(from);
318  if (ok) {
319  if (includeInBoundary) {
320  myConvBoundary.add(from);
321  }
322  }
323  return ok;
324 }
325 
326 
327 bool
329  double x = from.x() * myGeoScale;
330  double y = from.y() * myGeoScale;
331  if (myProjectionMethod == NONE) {
332  from.add(myOffset);
333  } else if (myUseInverseProjection) {
334  cartesian2geo(from);
335  } else {
336  if (x > 180.1 || x < -180.1 || y > 90.1 || y < -90.1) {
337  return false;
338  }
339 #ifdef HAVE_PROJ
340  if (myProjection != 0) {
341  projUV p;
342  p.u = x * DEG_TO_RAD;
343  p.v = y * DEG_TO_RAD;
344  p = pj_fwd(p, myProjection);
346  x = p.u;
347  y = p.v;
348  }
349 #endif
350  if (myProjectionMethod == SIMPLE) {
351  double ys = y;
352  x *= 111320. * cos(DEG2RAD(ys));
353  y *= 111136.;
354  from.set((SUMOReal)x, (SUMOReal)y);
356  from.add(myOffset);
357  }
358  }
361  return false;
362  }
363  if (myProjectionMethod != SIMPLE) {
364  from.set((SUMOReal)x, (SUMOReal)y);
365  from.add(myOffset);
366  }
367  return true;
368 }
369 
370 
371 void
373  myOffset.add(x, y);
374  myConvBoundary.moveby(x, y);
375 }
376 
377 
378 const Boundary&
380  return myOrigBoundary;
381 }
382 
383 
384 const Boundary&
386  return myConvBoundary;
387 }
388 
389 
390 const Position
392  return myOffset;
393 }
394 
395 
396 const Position
398  return myOffset;
399 }
400 
401 
402 const std::string&
404  return myProjString;
405 }
406 
407 
408 void
410  if (myNumLoaded == 0) {
412  } else {
414  // prefer options over loaded location
416  // let offset and boundary lead back to the original coords of the loaded data
419  // the new boundary (updated during loading)
421  }
422 }
423 
424 
425 void
427  myNumLoaded++;
428  if (myNumLoaded > 1) {
429  WRITE_WARNING("Ignoring loaded location attribute nr. " + toString(myNumLoaded) + " for tracking of original location");
430  } else {
431  myLoaded = loaded;
432  }
433 }
434 
435 
436 void
438  myNumLoaded = 0;
439 }
440 
441 
442 /****************************************************************************/
443 
void sub(SUMOReal dx, SUMOReal dy)
Substracts the given position from this one.
Definition: Position.h:139
void doRegister(const std::string &name, Option *v)
Adds an option under the given name.
Definition: OptionsCont.cpp:84
void cartesian2geo(Position &cartesian) const
Converts the given cartesian (shifted) position to its geo (lat/long) representation.
static GeoConvHelper myProcessing
coordinate transformation to use for input conversion and processing
bool x2cartesian_const(Position &from) const
Converts the given coordinate into a cartesian using the previous initialisation. ...
~GeoConvHelper()
Destructor.
void add(const Position &pos)
Adds the given position to this one.
Definition: Position.h:119
Position myOffset
The offset to apply.
bool x2cartesian(Position &from, bool includeInBoundary=true)
bool getBool(const std::string &name) const
Returns the boolean-value of the named option (only for Option_Bool)
Boundary myOrigBoundary
The boundary before conversion (x2cartesian)
static void setLoaded(const GeoConvHelper &loaded)
sets the coordinate transformation loaded from a location element
static GeoConvHelper myLoaded
coordinate transformation loaded from a location element
SUMOReal getFloat(const std::string &name) const
Returns the SUMOReal-value of the named option (only for Option_Float)
static void resetLoaded()
resets loaded location elements
bool myUseInverseProjection
Information whether inverse projection shall be used.
SUMOReal x() const
Returns the x-position.
Definition: Position.h:63
A class that stores a 2D geometrical boundary.
Definition: Boundary.h:48
#define WRITE_WARNING(msg)
Definition: MsgHandler.h:200
void addSynonyme(const std::string &name1, const std::string &name2, bool isDeprecated=false)
Adds a synonyme for an options name (any order)
double myGeoScale
The scaling to apply to geo-coordinates.
static void computeFinal()
compute the location attributes which will be used for output based on the loaded location data...
std::string getString(const std::string &name) const
Returns the string-value of the named option (only for Option_String)
#define max(a, b)
Definition: polyfonts.c:65
bool usingGeoProjection() const
Returns whether a transformation from geo to metric coordinates will be performed.
static methods for processing the coordinates conversion for the current net
Definition: GeoConvHelper.h:59
static GeoConvHelper myFinal
coordinate transformation to use for writing the location element and for tracking the original coord...
A point in 2D or 3D with translation and scaling methods.
Definition: Position.h:46
ProjectionMethod myProjectionMethod
Information whether no projection shall be done.
std::string myProjString
A proj options string describing the proj.4-projection to use.
const std::string & getProjString() const
Returns the network offset.
static bool init(OptionsCont &oc)
Initialises the processing and the final instance using the given options.
void addOptionSubTopic(const std::string &topic)
Adds an option subtopic.
const Position getOffsetBase() const
Returns the network base.
#define DEG2RAD(x)
Definition: GeomHelper.h:45
const Boundary & getConvBoundary() const
Returns the converted boundary.
std::string toString(const T &t, std::streamsize accuracy=OUTPUT_ACCURACY)
Definition: ToString.h:53
static void addProjectionOptions(OptionsCont &oc)
Adds projection options to the given container.
GeoConvHelper(OptionsCont &oc)
Constructor based on the stored options.
Boundary myConvBoundary
The boundary after conversion (x2cartesian)
#define WRITE_ERROR(msg)
Definition: MsgHandler.h:205
static int myNumLoaded
the numer of coordinate transformations loaded from location elements
void add(SUMOReal x, SUMOReal y)
Makes the boundary include the given coordinate.
Definition: Boundary.cpp:76
An integer-option.
Definition: Option.h:309
SUMOReal y() const
Returns the y-position.
Definition: Position.h:68
A storage for options typed value containers)
Definition: OptionsCont.h:108
void set(SUMOReal x, SUMOReal y)
Definition: Position.h:78
const Position getOffset() const
Returns the network offset.
#define SUMOReal
Definition: config.h:215
const Boundary & getOrigBoundary() const
Returns the original boundary.
void addDescription(const std::string &name, const std::string &subtopic, const std::string &description)
Adds a description for an option.
void moveby(SUMOReal x, SUMOReal y)
Moves the boundary by the given amount.
Definition: Boundary.cpp:249
bool usingInverseGeoProjection() const
Returns the information whether an inverse transformation will happen.
#define HAVE_PROJ
Definition: config.h:78
int getInt(const std::string &name) const
Returns the int-value of the named option (only for Option_Integer)
void moveConvertedBy(SUMOReal x, SUMOReal y)
Shifts the converted boundary by the given amounts.
bool exists(const std::string &name) const
Returns the information whether the named option is known.
GeoConvHelper & operator=(const GeoConvHelper &)
assignment operator.