Eclipse SUMO - Simulation of Urban MObility
GeoConvHelper.cpp
Go to the documentation of this file.
1 /****************************************************************************/
2 // Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.org/sumo
3 // Copyright (C) 2001-2019 German Aerospace Center (DLR) and others.
4 // This program and the accompanying materials
5 // are made available under the terms of the Eclipse Public License v2.0
6 // which accompanies this distribution, and is available at
7 // http://www.eclipse.org/legal/epl-v20.html
8 // SPDX-License-Identifier: EPL-2.0
9 /****************************************************************************/
17 // static methods for processing the coordinates conversion for the current net
18 /****************************************************************************/
19 
20 
21 // ===========================================================================
22 // included modules
23 // ===========================================================================
24 #include <config.h>
25 
26 #include <map>
27 #include <cmath>
28 #include <cassert>
29 #include <climits>
31 #include <utils/common/ToString.h>
32 #include <utils/geom/GeomHelper.h>
35 #include "GeoConvHelper.h"
36 
37 
38 // ===========================================================================
39 // static member variables
40 // ===========================================================================
41 
46 
47 // ===========================================================================
48 // method definitions
49 // ===========================================================================
50 
51 GeoConvHelper::GeoConvHelper(const std::string& proj, const Position& offset,
52  const Boundary& orig, const Boundary& conv, double scale, double rot, bool inverse, bool flatten):
53  myProjString(proj),
54 #ifdef PROJ_API_FILE
55  myProjection(nullptr),
56  myInverseProjection(nullptr),
57  myGeoProjection(nullptr),
58 #endif
59  myOffset(offset),
60  myGeoScale(scale),
61  mySin(sin(DEG2RAD(-rot))), // rotate clockwise
62  myCos(cos(DEG2RAD(-rot))),
63  myProjectionMethod(NONE),
64  myUseInverseProjection(inverse),
65  myFlatten(flatten),
66  myOrigBoundary(orig),
67  myConvBoundary(conv) {
68  if (proj == "!") {
70  } else if (proj == "-") {
72  } else if (proj == "UTM") {
74  } else if (proj == "DHDN") {
76  } else if (proj == "DHDN_UTM") {
78 #ifdef PROJ_API_FILE
79  } else {
81 #ifdef PROJ_VERSION_MAJOR
82  myProjection = proj_create(PJ_DEFAULT_CTX, proj.c_str());
83 #else
84  myProjection = pj_init_plus(proj.c_str());
85 #endif
86  if (myProjection == nullptr) {
87  // !!! check pj_errno
88  throw ProcessError("Could not build projection!");
89  }
90 #endif
91  }
92 }
93 
94 
96 #ifdef PROJ_API_FILE
97  if (myProjection != nullptr) {
98 #ifdef PROJ_VERSION_MAJOR
99  proj_destroy(myProjection);
100 #else
101  pj_free(myProjection);
102 #endif
103  }
104  if (myInverseProjection != nullptr) {
105 #ifdef PROJ_VERSION_MAJOR
106  proj_destroy(myInverseProjection);
107 #else
108  pj_free(myInverseProjection);
109 #endif
110  }
111  if (myGeoProjection != nullptr) {
112 #ifdef PROJ_VERSION_MAJOR
113  proj_destroy(myGeoProjection);
114 #else
115  pj_free(myGeoProjection);
116 #endif
117  }
118 #endif
119 }
120 
121 bool
123  return (
124  myProjString == o.myProjString &&
125  myOffset == o.myOffset &&
129  myGeoScale == o.myGeoScale &&
130  myCos == o.myCos &&
131  mySin == o.mySin &&
133  myFlatten == o.myFlatten
134  );
135 }
136 
139  myProjString = orig.myProjString;
140  myOffset = orig.myOffset;
144  myGeoScale = orig.myGeoScale;
145  myCos = orig.myCos;
146  mySin = orig.mySin;
148  myFlatten = orig.myFlatten;
149 #ifdef PROJ_API_FILE
150  if (myProjection != nullptr) {
151 #ifdef PROJ_VERSION_MAJOR
152  proj_destroy(myProjection);
153 #else
154  pj_free(myProjection);
155 #endif
156  myProjection = nullptr;
157  }
158  if (myInverseProjection != nullptr) {
159 #ifdef PROJ_VERSION_MAJOR
160  proj_destroy(myInverseProjection);
161 #else
162  pj_free(myInverseProjection);
163 #endif
164  myInverseProjection = nullptr;
165  }
166  if (myGeoProjection != nullptr) {
167 #ifdef PROJ_VERSION_MAJOR
168  proj_destroy(myGeoProjection);
169 #else
170  pj_free(myGeoProjection);
171 #endif
172  myGeoProjection = nullptr;
173  }
174  if (orig.myProjection != nullptr) {
175 #ifdef PROJ_VERSION_MAJOR
176  myProjection = proj_create(PJ_DEFAULT_CTX, orig.myProjString.c_str());
177 #else
178  myProjection = pj_init_plus(orig.myProjString.c_str());
179 #endif
180  }
181  if (orig.myInverseProjection != nullptr) {
182 #ifdef PROJ_VERSION_MAJOR
183  myInverseProjection = orig.myInverseProjection;
184 #else
185  myInverseProjection = pj_init_plus(pj_get_def(orig.myInverseProjection, 0));
186 #endif
187  }
188  if (orig.myGeoProjection != nullptr) {
189 #ifdef PROJ_VERSION_MAJOR
190  myGeoProjection = orig.myGeoProjection;
191 #else
192  myGeoProjection = pj_init_plus(pj_get_def(orig.myGeoProjection, 0));
193 #endif
194  }
195 #endif
196  return *this;
197 }
198 
199 
200 bool
202  std::string proj = "!"; // the default
203  double scale = oc.getFloat("proj.scale");
204  double rot = oc.getFloat("proj.rotate");
205  Position offset = Position(oc.getFloat("offset.x"), oc.getFloat("offset.y"));
206  bool inverse = oc.exists("proj.inverse") && oc.getBool("proj.inverse");
207  bool flatten = oc.exists("flatten") && oc.getBool("flatten");
208 
209  if (oc.getBool("simple-projection")) {
210  proj = "-";
211  }
212 
213 #ifdef PROJ_API_FILE
214  if (oc.getBool("proj.inverse") && oc.getString("proj") == "!") {
215  WRITE_ERROR("Inverse projection works only with explicit proj parameters.");
216  return false;
217  }
218  unsigned numProjections = oc.getBool("simple-projection") + oc.getBool("proj.utm") + oc.getBool("proj.dhdn") + oc.getBool("proj.dhdnutm") + (oc.getString("proj").length() > 1);
219  if (numProjections > 1) {
220  WRITE_ERROR("The projection method needs to be uniquely defined.");
221  return false;
222  }
223 
224  if (oc.getBool("proj.utm")) {
225  proj = "UTM";
226  } else if (oc.getBool("proj.dhdn")) {
227  proj = "DHDN";
228  } else if (oc.getBool("proj.dhdnutm")) {
229  proj = "DHDN_UTM";
230  } else if (!oc.isDefault("proj")) {
231  proj = oc.getString("proj");
232  }
233 #endif
234  myProcessing = GeoConvHelper(proj, offset, Boundary(), Boundary(), scale, rot, inverse, flatten);
236  return true;
237 }
238 
239 
240 void
241 GeoConvHelper::init(const std::string& proj, const Position& offset, const Boundary& orig,
242  const Boundary& conv, double scale) {
243  myProcessing = GeoConvHelper(proj, offset, orig, conv, scale);
245 }
246 
247 
248 void
250  oc.addOptionSubTopic("Projection");
251 
252  oc.doRegister("simple-projection", new Option_Bool(false));
253  oc.addSynonyme("simple-projection", "proj.simple", true);
254  oc.addDescription("simple-projection", "Projection", "Uses a simple method for projection");
255 
256  oc.doRegister("proj.scale", new Option_Float(1.0));
257  oc.addDescription("proj.scale", "Projection", "Scaling factor for input coordinates");
258 
259  oc.doRegister("proj.rotate", new Option_Float(0.0));
260  oc.addDescription("proj.rotate", "Projection", "Rotation (clockwise degrees) for input coordinates");
261 
262 #ifdef PROJ_API_FILE
263  oc.doRegister("proj.utm", new Option_Bool(false));
264  oc.addDescription("proj.utm", "Projection", "Determine the UTM zone (for a universal transversal mercator projection based on the WGS84 ellipsoid)");
265 
266  oc.doRegister("proj.dhdn", new Option_Bool(false));
267  oc.addDescription("proj.dhdn", "Projection", "Determine the DHDN zone (for a transversal mercator projection based on the bessel ellipsoid, \"Gauss-Krueger\")");
268 
269  oc.doRegister("proj", new Option_String("!"));
270  oc.addDescription("proj", "Projection", "Uses STR as proj.4 definition for projection");
271 
272  oc.doRegister("proj.inverse", new Option_Bool(false));
273  oc.addDescription("proj.inverse", "Projection", "Inverses projection");
274 
275  oc.doRegister("proj.dhdnutm", new Option_Bool(false));
276  oc.addDescription("proj.dhdnutm", "Projection", "Convert from Gauss-Krueger to UTM");
277 #endif // PROJ_API_FILE
278 }
279 
280 
281 bool
283  return myProjectionMethod != NONE;
284 }
285 
286 
287 bool
289  return myUseInverseProjection;
290 }
291 
292 
293 void
295  cartesian.sub(getOffsetBase());
296  if (myProjectionMethod == NONE) {
297  return;
298  }
299  if (myProjectionMethod == SIMPLE) {
300  const double y = cartesian.y() / 111136.;
301  const double x = cartesian.x() / 111320. / cos(DEG2RAD(y));
302  cartesian.set(x, y);
303  return;
304  }
305 #ifdef PROJ_API_FILE
306 #ifdef PROJ_VERSION_MAJOR
307  PJ_COORD c;
308  c.xy.x = cartesian.x();
309  c.xy.y = cartesian.y();
310  c = proj_trans(myProjection, PJ_INV, c);
311  cartesian.set(proj_todeg(c.lp.lam), proj_todeg(c.lp.phi));
312 #else
313  projUV p;
314  p.u = cartesian.x();
315  p.v = cartesian.y();
316  p = pj_inv(p, myProjection);
318  p.u *= RAD_TO_DEG;
319  p.v *= RAD_TO_DEG;
320  cartesian.set((double) p.u, (double) p.v);
321 #endif
322 #endif
323 }
324 
325 
326 bool
327 GeoConvHelper::x2cartesian(Position& from, bool includeInBoundary) {
328  if (includeInBoundary) {
329  myOrigBoundary.add(from);
330  }
331  // init projection parameter on first use
332 #ifdef PROJ_API_FILE
333  if (myProjection == nullptr) {
334  double x = from.x() * myGeoScale;
335  switch (myProjectionMethod) {
336  case DHDN_UTM: {
337  int zone = (int)((x - 500000.) / 1000000.);
338  if (zone < 1 || zone > 5) {
339  WRITE_WARNING("Attempt to initialize DHDN_UTM-projection on invalid longitude " + toString(x));
340  return false;
341  }
342  myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
343  " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
344  " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
345 #ifdef PROJ_VERSION_MAJOR
346  myInverseProjection = proj_create(PJ_DEFAULT_CTX, myProjString.c_str());
347  myGeoProjection = proj_create(PJ_DEFAULT_CTX, "+proj=latlong +datum=WGS84");
348 #else
349  myInverseProjection = pj_init_plus(myProjString.c_str());
350  myGeoProjection = pj_init_plus("+proj=latlong +datum=WGS84");
351 #endif
352  x = ((x - 500000.) / 1000000.) * 3; // continues with UTM
354  }
355  FALLTHROUGH;
356  case UTM: {
357  int zone = (int)(x + 180) / 6 + 1;
358  myProjString = "+proj=utm +zone=" + toString(zone) +
359  " +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
360 #ifdef PROJ_VERSION_MAJOR
361  myProjection = proj_create(PJ_DEFAULT_CTX, myProjString.c_str());
362 #else
363  myProjection = pj_init_plus(myProjString.c_str());
364 #endif
365  }
367  break;
368  case DHDN: {
369  int zone = (int)(x / 3);
370  if (zone < 1 || zone > 5) {
371  WRITE_WARNING("Attempt to initialize DHDN-projection on invalid longitude " + toString(x));
372  return false;
373  }
374  myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
375  " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
376  " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
377 #ifdef PROJ_VERSION_MAJOR
378  myProjection = proj_create(PJ_DEFAULT_CTX, myProjString.c_str());
379 #else
380  myProjection = pj_init_plus(myProjString.c_str());
381 #endif
382  }
384  break;
385  default:
386  break;
387  }
388  }
389  if (myInverseProjection != nullptr) {
390 #ifdef PROJ_VERSION_MAJOR
391  PJ_COORD c;
392  c.xy.x = from.x();
393  c.xy.y = from.y();
394  c = proj_trans(myInverseProjection, PJ_INV, c);
395  from.set(proj_todeg(c.lp.lam), proj_todeg(c.lp.phi));
396 #else
397  double x = from.x();
398  double y = from.y();
399  if (pj_transform(myInverseProjection, myGeoProjection, 1, 1, &x, &y, nullptr)) {
400  WRITE_WARNING("Could not transform (" + toString(x) + "," + toString(y) + ")");
401  }
402  from.set(double(x * RAD_TO_DEG), double(y * RAD_TO_DEG));
403 #endif
404  }
405 #endif
406  // perform conversion
407  bool ok = x2cartesian_const(from);
408  if (ok) {
409  if (includeInBoundary) {
410  myConvBoundary.add(from);
411  }
412  }
413  return ok;
414 }
415 
416 
417 bool
419  double x2 = from.x() * myGeoScale;
420  double y2 = from.y() * myGeoScale;
421  double x = x2 * myCos - y2 * mySin;
422  double y = x2 * mySin + y2 * myCos;
423  if (myProjectionMethod == NONE) {
424  from.add(myOffset);
425  } else if (myUseInverseProjection) {
426  cartesian2geo(from);
427  } else {
428  if (x > 180.1 || x < -180.1) {
429  WRITE_WARNING("Invalid longitude " + toString(x));
430  return false;
431  }
432  if (y > 90.1 || y < -90.1) {
433  WRITE_WARNING("Invalid latitude " + toString(y));
434  return false;
435  }
436 #ifdef PROJ_API_FILE
437  if (myProjection != nullptr) {
438 #ifdef PROJ_VERSION_MAJOR
439  PJ_COORD c;
440  c.lp.lam = proj_torad(x);
441  c.lp.phi = proj_torad(y);
442  c = proj_trans(myProjection, PJ_FWD, c);
444  x = c.xy.x;
445  y = c.xy.y;
446 #else
447  projUV p;
448  p.u = x * DEG_TO_RAD;
449  p.v = y * DEG_TO_RAD;
450  p = pj_fwd(p, myProjection);
452  x = p.u;
453  y = p.v;
454 #endif
455  }
456 #endif
457  if (myProjectionMethod == SIMPLE) {
458  x *= 111320. * cos(DEG2RAD(y));
459  y *= 111136.;
461  }
462  }
463  if (x > std::numeric_limits<double>::max() ||
464  y > std::numeric_limits<double>::max()) {
465  return false;
466  }
467  from.set(x, y);
468  from.add(myOffset);
469  if (myFlatten) {
470  from.setz(0);
471  }
472  return true;
473 }
474 
475 
476 void
477 GeoConvHelper::moveConvertedBy(double x, double y) {
478  myOffset.add(x, y);
479  myConvBoundary.moveby(x, y);
480 }
481 
482 
483 const Boundary&
485  return myOrigBoundary;
486 }
487 
488 
489 const Boundary&
491  return myConvBoundary;
492 }
493 
494 
495 const Position
497  return myOffset;
498 }
499 
500 
501 const Position
503  return myOffset;
504 }
505 
506 
507 const std::string&
509  return myProjString;
510 }
511 
512 
513 void
515  if (myNumLoaded == 0) {
517  if (lefthand) {
518  myFinal.myOffset.mul(1, -1);
519  }
520  } else {
521  if (lefthand) {
522  myProcessing.myOffset.mul(1, -1);
523  }
525  // prefer options over loaded location
527  // let offset and boundary lead back to the original coords of the loaded data
530  // the new boundary (updated during loading)
532  }
533  if (lefthand) {
535  }
536 }
537 
538 
539 void
541  myNumLoaded++;
542  if (myNumLoaded > 1) {
543  WRITE_WARNING("Ignoring loaded location attribute nr. " + toString(myNumLoaded) + " for tracking of original location");
544  } else {
545  myLoaded = loaded;
546  }
547 }
548 
549 
550 void
552  myNumLoaded = 0;
553 }
554 
555 
556 void
561  if (myFinal.usingGeoProjection()) {
563  }
565  if (myFinal.usingGeoProjection()) {
566  into.setPrecision();
567  }
569  into.closeTag();
570  into.lf();
571 }
572 
573 
574 /****************************************************************************/
void doRegister(const std::string &name, Option *v)
Adds an option under the given name.
Definition: OptionsCont.cpp:75
OutputDevice & writeAttr(const SumoXMLAttr attr, const T &val)
writes a named attribute
Definition: OutputDevice.h:256
static void writeLocation(OutputDevice &into)
writes the location element
static GeoConvHelper myProcessing
coordinate transformation to use for input conversion and processing
~GeoConvHelper()
Destructor.
const Boundary & getConvBoundary() const
Returns the converted boundary.
void add(const Position &pos)
Adds the given position to this one.
Definition: Position.h:127
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...
bool x2cartesian(Position &from, bool includeInBoundary=true)
Converts the given coordinate into a cartesian and optionally update myConvBoundary.
double y() const
Returns the y-position.
Definition: Position.h:62
bool usingGeoProjection() const
Returns whether a transformation from geo to metric coordinates will be performed.
void moveby(double x, double y, double z=0)
Moves the boundary by the given amount.
Definition: Boundary.cpp:369
double mySin
The rotation to apply to geo-coordinates.
Boundary myOrigBoundary
The boundary before conversion (x2cartesian)
static void setLoaded(const GeoConvHelper &loaded)
sets the coordinate transformation loaded from a location element
double x() const
Returns the x-position.
Definition: Position.h:57
void setPrecision(int precision=gPrecision)
Sets the precison or resets it to default.
const Boundary & getOrigBoundary() const
Returns the original boundary.
bool myFlatten
whether to discard z-data
static GeoConvHelper myLoaded
coordinate transformation loaded from a location element
bool getBool(const std::string &name) const
Returns the boolean-value of the named option (only for Option_Bool)
const std::string & getProjString() const
Returns the original projection definition.
static void resetLoaded()
resets loaded location elements
bool isDefault(const std::string &name) const
Returns the information whether the named option has still the default value.
void set(double x, double y)
set positions x and y
Definition: Position.h:87
bool myUseInverseProjection
Information whether inverse projection shall be used.
void moveConvertedBy(double x, double y)
Shifts the converted boundary by the given amounts.
A class that stores a 2D geometrical boundary.
Definition: Boundary.h:42
#define WRITE_WARNING(msg)
Definition: MsgHandler.h:239
void addSynonyme(const std::string &name1, const std::string &name2, bool isDeprecated=false)
Adds a synonyme for an options name (any order)
Definition: OptionsCont.cpp:96
void cartesian2geo(Position &cartesian) const
Converts the given cartesian (shifted) position to its geo (lat/long) representation.
double myGeoScale
The scaling to apply to geo-coordinates.
static methods for processing the coordinates conversion for the current net
Definition: GeoConvHelper.h:56
std::string toString(const T &t, std::streamsize accuracy=gPrecision)
Definition: ToString.h:48
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:39
int gPrecisionGeo
Definition: StdDefs.cpp:28
ProjectionMethod myProjectionMethod
Information whether no projection shall be done.
std::string myProjString
A proj options string describing the proj.4-projection to use.
std::string getString(const std::string &name) const
Returns the string-value of the named option (only for Option_String)
bool exists(const std::string &name) const
Returns the information whether the named option is known.
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.
#define DEG2RAD(x)
Definition: GeomHelper.h:38
static void addProjectionOptions(OptionsCont &oc)
Adds projection options to the given container.
GeoConvHelper(OptionsCont &oc)
Constructor based on the stored options.
double getFloat(const std::string &name) const
Returns the double-value of the named option (only for Option_Float)
Boundary myConvBoundary
The boundary after conversion (x2cartesian)
#define WRITE_ERROR(msg)
Definition: MsgHandler.h:245
static int myNumLoaded
the numer of coordinate transformations loaded from location elements
bool operator==(const GeoConvHelper &o) const
void flipY()
flips ymin and ymax
Definition: Boundary.cpp:323
A storage for options typed value containers)
Definition: OptionsCont.h:90
bool usingInverseGeoProjection() const
Returns the information whether an inverse transformation will happen.
const Position getOffsetBase() const
Returns the network base.
Static storage of an output device and its base (abstract) implementation.
Definition: OutputDevice.h:64
bool x2cartesian_const(Position &from) const
Converts the given coordinate into a cartesian using the previous initialisation. ...
bool closeTag(const std::string &comment="")
Closes the most recently opened tag and optionally adds a comment.
void mul(double val)
Multiplies both positions with the given value.
Definition: Position.h:107
const Position getOffset() const
Returns the network offset.
void addDescription(const std::string &name, const std::string &subtopic, const std::string &description)
Adds a description for an option.
void add(double x, double y, double z=0)
Makes the boundary include the given coordinate.
Definition: Boundary.cpp:79
OutputDevice & openTag(const std::string &xmlElement)
Opens an XML tag.
void lf()
writes a line feed if applicable
Definition: OutputDevice.h:234
GeoConvHelper & operator=(const GeoConvHelper &)
make assignment operator private
void setz(double z)
set position z
Definition: Position.h:82
void sub(double dx, double dy)
Substracts the given position from this one.
Definition: Position.h:147
#define FALLTHROUGH
Definition: StdDefs.h:37