SUMO - Simulation of Urban MObility
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules 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.dlr.de/
12 // Copyright (C) 2001-2016 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>
42 #include "GeoConvHelper.h"
43 
44 #ifdef CHECK_MEMORY_LEAKS
45 #include <foreign/nvwa/debug_new.h>
46 #endif // CHECK_MEMORY_LEAKS
47 
48 
49 // ===========================================================================
50 // static member variables
51 // ===========================================================================
56 
57 // ===========================================================================
58 // method definitions
59 // ===========================================================================
60 GeoConvHelper::GeoConvHelper(const std::string& proj, const Position& offset,
61  const Boundary& orig, const Boundary& conv, int shift, bool inverse):
62  myProjString(proj),
63 #ifdef HAVE_PROJ
64  myProjection(0),
65  myInverseProjection(0),
66  myGeoProjection(0),
67 #endif
68  myOffset(offset),
69  myGeoScale(pow(10, (double) - shift)),
70  myProjectionMethod(NONE),
71  myUseInverseProjection(inverse),
72  myOrigBoundary(orig),
73  myConvBoundary(conv) {
74  if (proj == "!") {
76  } else if (proj == "-") {
78  } else if (proj == "UTM") {
80  } else if (proj == "DHDN") {
82  } else if (proj == "DHDN_UTM") {
84 #ifdef HAVE_PROJ
85  } else {
87  myProjection = pj_init_plus(proj.c_str());
88  if (myProjection == 0) {
89  // !!! check pj_errno
90  throw ProcessError("Could not build projection!");
91  }
92 #endif
93  }
94 }
95 
96 
98 #ifdef HAVE_PROJ
99  if (myProjection != 0) {
100  pj_free(myProjection);
101  }
102  if (myInverseProjection != 0) {
103  pj_free(myInverseProjection);
104  }
105  if (myGeoProjection != 0) {
106  pj_free(myInverseProjection);
107  }
108 #endif
109 }
110 
111 
114  myProjString = orig.myProjString;
115  myOffset = orig.myOffset;
119  myGeoScale = orig.myGeoScale;
121 #ifdef HAVE_PROJ
122  if (myProjection != 0) {
123  pj_free(myProjection);
124  myProjection = 0;
125  }
126  if (myInverseProjection != 0) {
127  pj_free(myInverseProjection);
128  myInverseProjection = 0;
129  }
130  if (myGeoProjection != 0) {
131  pj_free(myGeoProjection);
132  myGeoProjection = 0;
133  }
134  if (orig.myProjection != 0) {
135  myProjection = pj_init_plus(orig.myProjString.c_str());
136  }
137  if (orig.myInverseProjection != 0) {
138  myInverseProjection = pj_init_plus(pj_get_def(orig.myInverseProjection, 0));
139  }
140  if (orig.myGeoProjection != 0) {
141  myGeoProjection = pj_init_plus(pj_get_def(orig.myGeoProjection, 0));
142  }
143 #endif
144  return *this;
145 }
146 
147 
148 bool
150  std::string proj = "!"; // the default
151  int shift = oc.getInt("proj.scale");
152  Position offset = Position(oc.getFloat("offset.x"), oc.getFloat("offset.y"));
153  bool inverse = oc.exists("proj.inverse") && oc.getBool("proj.inverse");
154 
155  if (oc.getBool("simple-projection")) {
156  proj = "-";
157  }
158 
159 #ifdef HAVE_PROJ
160  if (oc.getBool("proj.inverse") && oc.getString("proj") == "!") {
161  WRITE_ERROR("Inverse projection works only with explicit proj parameters.");
162  return false;
163  }
164  unsigned numProjections = oc.getBool("simple-projection") + oc.getBool("proj.utm") + oc.getBool("proj.dhdn") + oc.getBool("proj.dhdnutm") + (oc.getString("proj").length() > 1);
165  if (numProjections > 1) {
166  WRITE_ERROR("The projection method needs to be uniquely defined.");
167  return false;
168  }
169 
170  if (oc.getBool("proj.utm")) {
171  proj = "UTM";
172  } else if (oc.getBool("proj.dhdn")) {
173  proj = "DHDN";
174  } else if (oc.getBool("proj.dhdnutm")) {
175  proj = "DHDN_UTM";
176  } else if (!oc.isDefault("proj")) {
177  proj = oc.getString("proj");
178  }
179 #endif
180  myProcessing = GeoConvHelper(proj, offset, Boundary(), Boundary(), shift, inverse);
182  return true;
183 }
184 
185 
186 void
187 GeoConvHelper::init(const std::string& proj,
188  const Position& offset,
189  const Boundary& orig,
190  const Boundary& conv,
191  int shift) {
192  myProcessing = GeoConvHelper(proj, offset, orig, conv, shift);
194 }
195 
196 
197 void
199  oc.addOptionSubTopic("Projection");
200 
201  oc.doRegister("simple-projection", new Option_Bool(false));
202  oc.addSynonyme("simple-projection", "proj.simple", true);
203  oc.addDescription("simple-projection", "Projection", "Uses a simple method for projection");
204 
205  oc.doRegister("proj.scale", new Option_Integer(0));
206  oc.addSynonyme("proj.scale", "proj.shift", true);
207  oc.addDescription("proj.scale", "Projection", "Number of places to shift decimal point to right in geo-coordinates");
208 
209 #ifdef HAVE_PROJ
210  oc.doRegister("proj.utm", new Option_Bool(false));
211  oc.addDescription("proj.utm", "Projection", "Determine the UTM zone (for a universal transversal mercator projection based on the WGS84 ellipsoid)");
212 
213  oc.doRegister("proj.dhdn", new Option_Bool(false));
214  oc.addDescription("proj.dhdn", "Projection", "Determine the DHDN zone (for a transversal mercator projection based on the bessel ellipsoid, \"Gauss-Krueger\")");
215 
216  oc.doRegister("proj", new Option_String("!"));
217  oc.addDescription("proj", "Projection", "Uses STR as proj.4 definition for projection");
218 
219  oc.doRegister("proj.inverse", new Option_Bool(false));
220  oc.addDescription("proj.inverse", "Projection", "Inverses projection");
221 
222  oc.doRegister("proj.dhdnutm", new Option_Bool(false));
223  oc.addDescription("proj.dhdnutm", "Projection", "Convert from Gauss-Krueger to UTM");
224 #endif // HAVE_PROJ
225 }
226 
227 
228 bool
230  return myProjectionMethod != NONE;
231 }
232 
233 
234 bool
236  return myUseInverseProjection;
237 }
238 
239 
240 void
242  cartesian.sub(getOffsetBase());
243  if (myProjectionMethod == NONE) {
244  return;
245  }
246 #ifdef HAVE_PROJ
247  projUV p;
248  p.u = cartesian.x();
249  p.v = cartesian.y();
250  p = pj_inv(p, myProjection);
252  p.u *= RAD_TO_DEG;
253  p.v *= RAD_TO_DEG;
254  cartesian.set((SUMOReal) p.u, (SUMOReal) p.v);
255 #endif
256 }
257 
258 
259 bool
260 GeoConvHelper::x2cartesian(Position& from, bool includeInBoundary) {
261  if (includeInBoundary) {
262  myOrigBoundary.add(from);
263  }
264  // init projection parameter on first use
265 #ifdef HAVE_PROJ
266  if (myProjection == 0) {
267  double x = from.x() * myGeoScale;
268  switch (myProjectionMethod) {
269  case DHDN_UTM: {
270  int zone = (int)((x - 500000.) / 1000000.);
271  if (zone < 1 || zone > 5) {
272  WRITE_WARNING("Attempt to initialize DHDN_UTM-projection on invalid longitude " + toString(x));
273  return false;
274  }
275  myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
276  " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
277  " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
278  myInverseProjection = pj_init_plus(myProjString.c_str());
279  myGeoProjection = pj_init_plus("+proj=latlong +datum=WGS84");
281  x = ((x - 500000.) / 1000000.) * 3; // continues with UTM
282  }
283  case UTM: {
284  int zone = (int)(x + 180) / 6 + 1;
285  myProjString = "+proj=utm +zone=" + toString(zone) +
286  " +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
287  myProjection = pj_init_plus(myProjString.c_str());
289  }
290  break;
291  case DHDN: {
292  int zone = (int)(x / 3);
293  if (zone < 1 || zone > 5) {
294  WRITE_WARNING("Attempt to initialize DHDN-projection on invalid longitude " + toString(x));
295  return false;
296  }
297  myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
298  " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
299  " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
300  myProjection = pj_init_plus(myProjString.c_str());
302  }
303  break;
304  default:
305  break;
306  }
307  }
308  if (myInverseProjection != 0) {
309  double x = from.x();
310  double y = from.y();
311  if (pj_transform(myInverseProjection, myGeoProjection, 1, 1, &x, &y, NULL)) {
312  WRITE_WARNING("Could not transform (" + toString(x) + "," + toString(y) + ")");
313  }
314  from.set(SUMOReal(x * RAD_TO_DEG), SUMOReal(y * RAD_TO_DEG));
315  }
316 #endif
317  // perform conversion
318  bool ok = x2cartesian_const(from);
319  if (ok) {
320  if (includeInBoundary) {
321  myConvBoundary.add(from);
322  }
323  }
324  return ok;
325 }
326 
327 
328 bool
330  double x = from.x() * myGeoScale;
331  double y = from.y() * myGeoScale;
332  if (myProjectionMethod == NONE) {
333  from.add(myOffset);
334  } else if (myUseInverseProjection) {
335  cartesian2geo(from);
336  } else {
337  if (x > 180.1 || x < -180.1) {
338  WRITE_WARNING("Invalid longitude " + toString(x));
339  return false;
340  }
341  if (y > 90.1 || y < -90.1) {
342  WRITE_WARNING("Invalid latitude " + toString(y));
343  return false;
344  }
345 #ifdef HAVE_PROJ
346  if (myProjection != 0) {
347  projUV p;
348  p.u = x * DEG_TO_RAD;
349  p.v = y * DEG_TO_RAD;
350  p = pj_fwd(p, myProjection);
352  x = p.u;
353  y = p.v;
354  }
355 #endif
356  if (myProjectionMethod == SIMPLE) {
357  double ys = y;
358  x *= 111320. * cos(DEG2RAD(ys));
359  y *= 111136.;
360  from.set((SUMOReal)x, (SUMOReal)y);
362  from.add(myOffset);
363  }
364  }
367  return false;
368  }
369  if (myProjectionMethod != SIMPLE) {
370  from.set((SUMOReal)x, (SUMOReal)y);
371  from.add(myOffset);
372  }
373  return true;
374 }
375 
376 
377 void
379  myOffset.add(x, y);
380  myConvBoundary.moveby(x, y);
381 }
382 
383 
384 const Boundary&
386  return myOrigBoundary;
387 }
388 
389 
390 const Boundary&
392  return myConvBoundary;
393 }
394 
395 
396 const Position
398  return myOffset;
399 }
400 
401 
402 const Position
404  return myOffset;
405 }
406 
407 
408 const std::string&
410  return myProjString;
411 }
412 
413 
414 void
416  if (myNumLoaded == 0) {
418  } else {
420  // prefer options over loaded location
422  // let offset and boundary lead back to the original coords of the loaded data
425  // the new boundary (updated during loading)
427  }
428  if (lefthand) {
429  myFinal.myOffset.mul(1, -1);
431  }
432 }
433 
434 
435 void
437  myNumLoaded++;
438  if (myNumLoaded > 1) {
439  WRITE_WARNING("Ignoring loaded location attribute nr. " + toString(myNumLoaded) + " for tracking of original location");
440  } else {
441  myLoaded = loaded;
442  }
443 }
444 
445 
446 void
448  myNumLoaded = 0;
449 }
450 
451 
452 void
457  if (myFinal.usingGeoProjection()) {
459  }
461  if (myFinal.usingGeoProjection()) {
462  into.setPrecision();
463  }
465  into.closeTag();
466  into.lf();
467 }
468 
469 
470 
471 /****************************************************************************/
472 
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:86
OutputDevice & writeAttr(const SumoXMLAttr attr, const T &val)
writes a named attribute
Definition: OutputDevice.h:257
static void writeLocation(OutputDevice &into)
writes the location element
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.
static void computeFinal(bool lefthand=false)
compute the location attributes which will be used for output based on the loaded location data...
#define GEO_OUTPUT_ACCURACY
Definition: config.h:16
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)
void setPrecision(int precision=OUTPUT_ACCURACY)
Sets the precison or resets it to default.
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.
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:60
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.
bool isDefault(const std::string &name) const
Returns the information whether the named option has still the default value.
#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:55
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:206
static int myNumLoaded
the numer of coordinate transformations loaded from location elements
An integer-option.
Definition: Option.h:313
void add(SUMOReal x, SUMOReal y, SUMOReal z=0)
Makes the boundary include the given coordinate.
Definition: Boundary.cpp:90
void flipY()
flips ymin and ymax
Definition: Boundary.cpp:254
SUMOReal y() const
Returns the y-position.
Definition: Position.h:68
A storage for options typed value containers)
Definition: OptionsCont.h:99
void set(SUMOReal x, SUMOReal y)
Definition: Position.h:78
void mul(SUMOReal val)
Multiplies both positions with the given value.
Definition: Position.h:99
const Position getOffset() const
Returns the network offset.
Static storage of an output device and its base (abstract) implementation.
Definition: OutputDevice.h:71
bool closeTag()
Closes the most recently opened tag.
#define SUMOReal
Definition: config.h:214
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.
bool usingInverseGeoProjection() const
Returns the information whether an inverse transformation will happen.
#define HAVE_PROJ
Definition: config.h:83
int getInt(const std::string &name) const
Returns the int-value of the named option (only for Option_Integer)
void moveby(SUMOReal x, SUMOReal y, SUMOReal z=0)
Moves the boundary by the given amount.
Definition: Boundary.cpp:281
OutputDevice & openTag(const std::string &xmlElement)
Opens an XML tag.
void lf()
writes a line feed if applicable
Definition: OutputDevice.h:235
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.