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