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  myProjString = std::regex_replace(proj, std::regex("\\+geoidgrids[^ ]*"), std::string(""));
84  myProjString = std::regex_replace(myProjString, std::regex("\\+step \\+proj=vgridshift \\+grids[^ ]*"), std::string(""));
85  if (myProjString != proj) {
86  WRITE_WARNING(TL("Ignoring geoidgrids and vgridshift in projection"));
87  initProj(myProjString);
88  }
89  }
90  if (myProjection == nullptr) {
91  // !!! check pj_errno
92  throw ProcessError("Could not build projection!");
93  }
94 #endif
95  }
96 }
97 
98 
99 #ifdef PROJ_API_FILE
100 void
101 GeoConvHelper::initProj(const std::string& proj) {
102 #ifdef PROJ_VERSION_MAJOR
103  myProjection = proj_create(PJ_DEFAULT_CTX, proj.c_str());
104 #else
105  myProjection = pj_init_plus(proj.c_str());
106 #endif
107 }
108 #endif
109 
110 
112 #ifdef PROJ_API_FILE
113  if (myProjection != nullptr) {
114 #ifdef PROJ_VERSION_MAJOR
115  proj_destroy(myProjection);
116 #else
117  pj_free(myProjection);
118 #endif
119  }
120  if (myInverseProjection != nullptr) {
121 #ifdef PROJ_VERSION_MAJOR
122  proj_destroy(myInverseProjection);
123 #else
124  pj_free(myInverseProjection);
125 #endif
126  }
127  if (myGeoProjection != nullptr) {
128 #ifdef PROJ_VERSION_MAJOR
129  proj_destroy(myGeoProjection);
130 #else
131  pj_free(myGeoProjection);
132 #endif
133  }
134 #endif
135 }
136 
137 bool
139  return (
140  myProjString == o.myProjString &&
141  myOffset == o.myOffset &&
145  myGeoScale == o.myGeoScale &&
146  myCos == o.myCos &&
147  mySin == o.mySin &&
149  myFlatten == o.myFlatten
150  );
151 }
152 
155  myProjString = orig.myProjString;
156  myOffset = orig.myOffset;
160  myGeoScale = orig.myGeoScale;
161  myCos = orig.myCos;
162  mySin = orig.mySin;
164  myFlatten = orig.myFlatten;
165 #ifdef PROJ_API_FILE
166  if (myProjection != nullptr) {
167 #ifdef PROJ_VERSION_MAJOR
168  proj_destroy(myProjection);
169 #else
170  pj_free(myProjection);
171 #endif
172  myProjection = nullptr;
173  }
174  if (myInverseProjection != nullptr) {
175 #ifdef PROJ_VERSION_MAJOR
176  proj_destroy(myInverseProjection);
177 #else
178  pj_free(myInverseProjection);
179 #endif
180  myInverseProjection = nullptr;
181  }
182  if (myGeoProjection != nullptr) {
183 #ifdef PROJ_VERSION_MAJOR
184  proj_destroy(myGeoProjection);
185 #else
186  pj_free(myGeoProjection);
187 #endif
188  myGeoProjection = nullptr;
189  }
190  if (orig.myProjection != nullptr) {
191 #ifdef PROJ_VERSION_MAJOR
192  myProjection = proj_create(PJ_DEFAULT_CTX, orig.myProjString.c_str());
193 #else
194  myProjection = pj_init_plus(orig.myProjString.c_str());
195 #endif
196  }
197  if (orig.myInverseProjection != nullptr) {
198 #ifdef PROJ_VERSION_MAJOR
199  myInverseProjection = orig.myInverseProjection;
200 #else
201  myInverseProjection = pj_init_plus(pj_get_def(orig.myInverseProjection, 0));
202 #endif
203  }
204  if (orig.myGeoProjection != nullptr) {
205 #ifdef PROJ_VERSION_MAJOR
206  myGeoProjection = orig.myGeoProjection;
207 #else
208  myGeoProjection = pj_init_plus(pj_get_def(orig.myGeoProjection, 0));
209 #endif
210  }
211 #endif
212  return *this;
213 }
214 
215 
216 bool
218  std::string proj = "!"; // the default
219  double scale = oc.getFloat("proj.scale");
220  double rot = oc.getFloat("proj.rotate");
221  Position offset = Position(oc.getFloat("offset.x"), oc.getFloat("offset.y"), oc.getFloat("offset.z"));
222  bool inverse = oc.exists("proj.inverse") && oc.getBool("proj.inverse");
223  bool flatten = oc.exists("flatten") && oc.getBool("flatten");
224 
225  if (oc.getBool("simple-projection")) {
226  proj = "-";
227  }
228 
229 #ifdef PROJ_API_FILE
230  if (oc.getBool("proj.inverse") && oc.getString("proj") == "!") {
231  WRITE_ERROR(TL("Inverse projection works only with explicit proj parameters."));
232  return false;
233  }
234  unsigned numProjections = oc.getBool("simple-projection") + oc.getBool("proj.utm") + oc.getBool("proj.dhdn") + oc.getBool("proj.dhdnutm") + (oc.getString("proj").length() > 1);
235  if (numProjections > 1) {
236  WRITE_ERROR(TL("The projection method needs to be uniquely defined."));
237  return false;
238  }
239 
240  if (oc.getBool("proj.utm")) {
241  proj = "UTM";
242  } else if (oc.getBool("proj.dhdn")) {
243  proj = "DHDN";
244  } else if (oc.getBool("proj.dhdnutm")) {
245  proj = "DHDN_UTM";
246  } else if (!oc.isDefault("proj")) {
247  proj = oc.getString("proj");
248  }
249 #endif
250  myProcessing = GeoConvHelper(proj, offset, Boundary(), Boundary(), scale, rot, inverse, flatten);
252  return true;
253 }
254 
255 
256 void
257 GeoConvHelper::init(const std::string& proj, const Position& offset, const Boundary& orig,
258  const Boundary& conv, double scale) {
259  myProcessing = GeoConvHelper(proj, offset, orig, conv, scale);
261 }
262 
263 
264 void
266  oc.addOptionSubTopic("Projection");
267 
268  oc.doRegister("simple-projection", new Option_Bool(false));
269  oc.addSynonyme("simple-projection", "proj.simple", true);
270  oc.addDescription("simple-projection", "Projection", "Uses a simple method for projection");
271 
272  oc.doRegister("proj.scale", new Option_Float(1.0));
273  oc.addDescription("proj.scale", "Projection", "Scaling factor for input coordinates");
274 
275  oc.doRegister("proj.rotate", new Option_Float(0.0));
276  oc.addDescription("proj.rotate", "Projection", "Rotation (clockwise degrees) for input coordinates");
277 
278 #ifdef PROJ_API_FILE
279  oc.doRegister("proj.utm", new Option_Bool(false));
280  oc.addDescription("proj.utm", "Projection", "Determine the UTM zone (for a universal transversal mercator projection based on the WGS84 ellipsoid)");
281 
282  oc.doRegister("proj.dhdn", new Option_Bool(false));
283  oc.addDescription("proj.dhdn", "Projection", "Determine the DHDN zone (for a transversal mercator projection based on the bessel ellipsoid, \"Gauss-Krueger\")");
284 
285  oc.doRegister("proj", new Option_String("!"));
286  oc.addDescription("proj", "Projection", "Uses STR as proj.4 definition for projection");
287 
288  oc.doRegister("proj.inverse", new Option_Bool(false));
289  oc.addDescription("proj.inverse", "Projection", "Inverses projection");
290 
291  oc.doRegister("proj.dhdnutm", new Option_Bool(false));
292  oc.addDescription("proj.dhdnutm", "Projection", "Convert from Gauss-Krueger to UTM");
293 #endif // PROJ_API_FILE
294 }
295 
296 
297 bool
299  return myProjectionMethod != NONE;
300 }
301 
302 
303 bool
305  return myUseInverseProjection;
306 }
307 
308 
309 void
311  cartesian.sub(getOffsetBase());
312  if (myProjectionMethod == NONE) {
313  return;
314  }
315  if (myProjectionMethod == SIMPLE) {
316  const double y = cartesian.y() / 111136.;
317  const double x = cartesian.x() / 111320. / cos(DEG2RAD(y));
318  cartesian.set(x, y);
319  return;
320  }
321 #ifdef PROJ_API_FILE
322 #ifdef PROJ_VERSION_MAJOR
323  PJ_COORD c;
324  c.xy.x = cartesian.x();
325  c.xy.y = cartesian.y();
326  c = proj_trans(myProjection, PJ_INV, c);
327  cartesian.set(proj_todeg(c.lp.lam), proj_todeg(c.lp.phi));
328 #else
329  projUV p;
330  p.u = cartesian.x();
331  p.v = cartesian.y();
332  p = pj_inv(p, myProjection);
334  p.u *= RAD_TO_DEG;
335  p.v *= RAD_TO_DEG;
336  cartesian.set((double) p.u, (double) p.v);
337 #endif
338 #endif
339 }
340 
341 
342 bool
343 GeoConvHelper::x2cartesian(Position& from, bool includeInBoundary) {
344  if (includeInBoundary) {
345  myOrigBoundary.add(from);
346  }
347  // init projection parameter on first use
348 #ifdef PROJ_API_FILE
349  if (myProjection == nullptr) {
350  double x = from.x() * myGeoScale;
351  switch (myProjectionMethod) {
352  case DHDN_UTM: {
353  int zone = (int)((x - 500000.) / 1000000.);
354  if (zone < 1 || zone > 5) {
355  WRITE_WARNING("Attempt to initialize DHDN_UTM-projection on invalid longitude " + toString(x));
356  return false;
357  }
358  myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
359  " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
360  " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
361 #ifdef PROJ_VERSION_MAJOR
362  myInverseProjection = proj_create(PJ_DEFAULT_CTX, myProjString.c_str());
363  myGeoProjection = proj_create(PJ_DEFAULT_CTX, "+proj=latlong +datum=WGS84");
364 #else
365  myInverseProjection = pj_init_plus(myProjString.c_str());
366  myGeoProjection = pj_init_plus("+proj=latlong +datum=WGS84");
367 #endif
369  x = ((x - 500000.) / 1000000.) * 3; // continues with UTM
370  }
371  FALLTHROUGH;
372  case UTM: {
373  int zone = (int)(x + 180) / 6 + 1;
374  myProjString = "+proj=utm +zone=" + toString(zone) +
375  " +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
376 #ifdef PROJ_VERSION_MAJOR
377  myProjection = proj_create(PJ_DEFAULT_CTX, myProjString.c_str());
378 #else
379  myProjection = pj_init_plus(myProjString.c_str());
380 #endif
382  }
383  break;
384  case DHDN: {
385  int zone = (int)(x / 3);
386  if (zone < 1 || zone > 5) {
387  WRITE_WARNING("Attempt to initialize DHDN-projection on invalid longitude " + toString(x));
388  return false;
389  }
390  myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
391  " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
392  " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
393 #ifdef PROJ_VERSION_MAJOR
394  myProjection = proj_create(PJ_DEFAULT_CTX, myProjString.c_str());
395 #else
396  myProjection = pj_init_plus(myProjString.c_str());
397 #endif
399  }
400  break;
401  default:
402  break;
403  }
404  }
405  if (myInverseProjection != nullptr) {
406 #ifdef PROJ_VERSION_MAJOR
407  PJ_COORD c;
408  c.xy.x = from.x();
409  c.xy.y = from.y();
410  c = proj_trans(myInverseProjection, PJ_INV, c);
411  from.set(proj_todeg(c.lp.lam), proj_todeg(c.lp.phi));
412 #else
413  double x = from.x();
414  double y = from.y();
415  if (pj_transform(myInverseProjection, myGeoProjection, 1, 1, &x, &y, nullptr)) {
416  WRITE_WARNING("Could not transform (" + toString(x) + "," + toString(y) + ")");
417  }
418  from.set(double(x * RAD_TO_DEG), double(y * RAD_TO_DEG));
419 #endif
420  }
421 #endif
422  // perform conversion
423  bool ok = x2cartesian_const(from);
424  if (ok) {
425  if (includeInBoundary) {
426  myConvBoundary.add(from);
427  }
428  }
429  return ok;
430 }
431 
432 
433 bool
435  double x2 = from.x() * myGeoScale;
436  double y2 = from.y() * myGeoScale;
437  double x = x2 * myCos - y2 * mySin;
438  double y = x2 * mySin + y2 * myCos;
439  if (myProjectionMethod == NONE) {
440  // do nothing
441  } else if (myUseInverseProjection) {
442  cartesian2geo(from);
443  } else {
444  if (x > 180.1 || x < -180.1) {
445  WRITE_WARNING("Invalid longitude " + toString(x));
446  return false;
447  }
448  if (y > 90.1 || y < -90.1) {
449  WRITE_WARNING("Invalid latitude " + toString(y));
450  return false;
451  }
452 #ifdef PROJ_API_FILE
453  if (myProjection != nullptr) {
454 #ifdef PROJ_VERSION_MAJOR
455  PJ_COORD c;
456  c.lp.lam = proj_torad(x);
457  c.lp.phi = proj_torad(y);
458  c = proj_trans(myProjection, PJ_FWD, c);
460  x = c.xy.x;
461  y = c.xy.y;
462 #else
463  projUV p;
464  p.u = x * DEG_TO_RAD;
465  p.v = y * DEG_TO_RAD;
466  p = pj_fwd(p, myProjection);
468  x = p.u;
469  y = p.v;
470 #endif
471  }
472 #endif
473  if (myProjectionMethod == SIMPLE) {
474  // Sinusoidal projection (https://en.wikipedia.org/wiki/Sinusoidal_projection)
475  x *= 111320. * cos(DEG2RAD(y));
476  y *= 111136.;
478  }
479  }
480  if (x > std::numeric_limits<double>::max() ||
481  y > std::numeric_limits<double>::max()) {
482  return false;
483  }
484  from.set(x, y);
485  from.add(myOffset);
486  if (myFlatten) {
487  from.setz(0);
488  }
489  return true;
490 }
491 
492 
493 void
494 GeoConvHelper::moveConvertedBy(double x, double y) {
495  myOffset.add(x, y);
496  myConvBoundary.moveby(x, y);
497 }
498 
499 
500 const Boundary&
502  return myOrigBoundary;
503 }
504 
505 
506 const Boundary&
508  return myConvBoundary;
509 }
510 
511 
512 const Position
514  return myOffset;
515 }
516 
517 
518 const Position
520  return myOffset;
521 }
522 
523 
524 const std::string&
526  return myProjString;
527 }
528 
529 
530 void
532  if (myNumLoaded == 0) {
534  if (lefthand) {
535  myFinal.myOffset.mul(1, -1);
536  }
537  } else {
538  if (lefthand) {
539  myProcessing.myOffset.mul(1, -1);
540  }
542  // prefer options over loaded location
544  // let offset and boundary lead back to the original coords of the loaded data
547  // the new boundary (updated during loading)
549  }
550  if (lefthand) {
552  }
553 }
554 
555 
556 void
558  myNumLoaded++;
559  if (myNumLoaded > 1) {
560  WRITE_WARNING("Ignoring loaded location attribute nr. " + toString(myNumLoaded) + " for tracking of original location");
561  } else {
562  myLoaded = loaded;
563  }
564 }
565 
566 
567 void
569  myNumLoaded = 0;
570 }
571 
572 
573 void
578  if (myFinal.usingGeoProjection()) {
580  }
582  if (myFinal.usingGeoProjection()) {
583  into.setPrecision();
584  }
586  into.closeTag();
587  into.lf();
588 }
589 
590 
591 /****************************************************************************/
#define DEG2RAD(x)
Definition: GeomHelper.h:35
#define WRITE_ERROR(msg)
Definition: MsgHandler.h:274
#define WRITE_WARNING(msg)
Definition: MsgHandler.h:265
#define TL(string)
Definition: MsgHandler.h:282
@ 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:78
void moveby(double x, double y, double z=0)
Moves the boundary by the given amount.
Definition: Boundary.cpp:393
void flipY()
flips ymin and ymax
Definition: Boundary.cpp:332
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:76
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:97
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:239
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:251
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