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-2020 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>
29 #include <utils/common/ToString.h>
30 #include <utils/geom/GeomHelper.h>
33 #include "GeoConvHelper.h"
34 
35 
36 // ===========================================================================
37 // static member variables
38 // ===========================================================================
39 
44 
45 // ===========================================================================
46 // method definitions
47 // ===========================================================================
48 
49 GeoConvHelper::GeoConvHelper(const std::string& proj, const Position& offset,
50  const Boundary& orig, const Boundary& conv, double scale, double rot, bool inverse, bool flatten):
51  myProjString(proj),
52 #ifdef PROJ_API_FILE
53  myProjection(nullptr),
54  myInverseProjection(nullptr),
55  myGeoProjection(nullptr),
56 #endif
57  myOffset(offset),
58  myGeoScale(scale),
59  mySin(sin(DEG2RAD(-rot))), // rotate clockwise
60  myCos(cos(DEG2RAD(-rot))),
61  myProjectionMethod(NONE),
62  myUseInverseProjection(inverse),
63  myFlatten(flatten),
64  myOrigBoundary(orig),
65  myConvBoundary(conv) {
66  if (proj == "!") {
68  } else if (proj == "-") {
70  } else if (proj == "UTM") {
72  } else if (proj == "DHDN") {
74  } else if (proj == "DHDN_UTM") {
76 #ifdef PROJ_API_FILE
77  } else {
79 #ifdef PROJ_VERSION_MAJOR
80  myProjection = proj_create(PJ_DEFAULT_CTX, proj.c_str());
81 #else
82  myProjection = pj_init_plus(proj.c_str());
83 #endif
84  if (myProjection == nullptr) {
85  // !!! check pj_errno
86  throw ProcessError("Could not build projection!");
87  }
88 #endif
89  }
90 }
91 
92 
94 #ifdef PROJ_API_FILE
95  if (myProjection != nullptr) {
96 #ifdef PROJ_VERSION_MAJOR
97  proj_destroy(myProjection);
98 #else
99  pj_free(myProjection);
100 #endif
101  }
102  if (myInverseProjection != nullptr) {
103 #ifdef PROJ_VERSION_MAJOR
104  proj_destroy(myInverseProjection);
105 #else
106  pj_free(myInverseProjection);
107 #endif
108  }
109  if (myGeoProjection != nullptr) {
110 #ifdef PROJ_VERSION_MAJOR
111  proj_destroy(myGeoProjection);
112 #else
113  pj_free(myGeoProjection);
114 #endif
115  }
116 #endif
117 }
118 
119 bool
121  return (
122  myProjString == o.myProjString &&
123  myOffset == o.myOffset &&
127  myGeoScale == o.myGeoScale &&
128  myCos == o.myCos &&
129  mySin == o.mySin &&
131  myFlatten == o.myFlatten
132  );
133 }
134 
137  myProjString = orig.myProjString;
138  myOffset = orig.myOffset;
142  myGeoScale = orig.myGeoScale;
143  myCos = orig.myCos;
144  mySin = orig.mySin;
146  myFlatten = orig.myFlatten;
147 #ifdef PROJ_API_FILE
148  if (myProjection != nullptr) {
149 #ifdef PROJ_VERSION_MAJOR
150  proj_destroy(myProjection);
151 #else
152  pj_free(myProjection);
153 #endif
154  myProjection = nullptr;
155  }
156  if (myInverseProjection != nullptr) {
157 #ifdef PROJ_VERSION_MAJOR
158  proj_destroy(myInverseProjection);
159 #else
160  pj_free(myInverseProjection);
161 #endif
162  myInverseProjection = nullptr;
163  }
164  if (myGeoProjection != nullptr) {
165 #ifdef PROJ_VERSION_MAJOR
166  proj_destroy(myGeoProjection);
167 #else
168  pj_free(myGeoProjection);
169 #endif
170  myGeoProjection = nullptr;
171  }
172  if (orig.myProjection != nullptr) {
173 #ifdef PROJ_VERSION_MAJOR
174  myProjection = proj_create(PJ_DEFAULT_CTX, orig.myProjString.c_str());
175 #else
176  myProjection = pj_init_plus(orig.myProjString.c_str());
177 #endif
178  }
179  if (orig.myInverseProjection != nullptr) {
180 #ifdef PROJ_VERSION_MAJOR
181  myInverseProjection = orig.myInverseProjection;
182 #else
183  myInverseProjection = pj_init_plus(pj_get_def(orig.myInverseProjection, 0));
184 #endif
185  }
186  if (orig.myGeoProjection != nullptr) {
187 #ifdef PROJ_VERSION_MAJOR
188  myGeoProjection = orig.myGeoProjection;
189 #else
190  myGeoProjection = pj_init_plus(pj_get_def(orig.myGeoProjection, 0));
191 #endif
192  }
193 #endif
194  return *this;
195 }
196 
197 
198 bool
200  std::string proj = "!"; // the default
201  double scale = oc.getFloat("proj.scale");
202  double rot = oc.getFloat("proj.rotate");
203  Position offset = Position(oc.getFloat("offset.x"), oc.getFloat("offset.y"), oc.getFloat("offset.z"));
204  bool inverse = oc.exists("proj.inverse") && oc.getBool("proj.inverse");
205  bool flatten = oc.exists("flatten") && oc.getBool("flatten");
206 
207  if (oc.getBool("simple-projection")) {
208  proj = "-";
209  }
210 
211 #ifdef PROJ_API_FILE
212  if (oc.getBool("proj.inverse") && oc.getString("proj") == "!") {
213  WRITE_ERROR("Inverse projection works only with explicit proj parameters.");
214  return false;
215  }
216  unsigned numProjections = oc.getBool("simple-projection") + oc.getBool("proj.utm") + oc.getBool("proj.dhdn") + oc.getBool("proj.dhdnutm") + (oc.getString("proj").length() > 1);
217  if (numProjections > 1) {
218  WRITE_ERROR("The projection method needs to be uniquely defined.");
219  return false;
220  }
221 
222  if (oc.getBool("proj.utm")) {
223  proj = "UTM";
224  } else if (oc.getBool("proj.dhdn")) {
225  proj = "DHDN";
226  } else if (oc.getBool("proj.dhdnutm")) {
227  proj = "DHDN_UTM";
228  } else if (!oc.isDefault("proj")) {
229  proj = oc.getString("proj");
230  }
231 #endif
232  myProcessing = GeoConvHelper(proj, offset, Boundary(), Boundary(), scale, rot, inverse, flatten);
234  return true;
235 }
236 
237 
238 void
239 GeoConvHelper::init(const std::string& proj, const Position& offset, const Boundary& orig,
240  const Boundary& conv, double scale) {
241  myProcessing = GeoConvHelper(proj, offset, orig, conv, scale);
243 }
244 
245 
246 void
248  oc.addOptionSubTopic("Projection");
249 
250  oc.doRegister("simple-projection", new Option_Bool(false));
251  oc.addSynonyme("simple-projection", "proj.simple", true);
252  oc.addDescription("simple-projection", "Projection", "Uses a simple method for projection");
253 
254  oc.doRegister("proj.scale", new Option_Float(1.0));
255  oc.addDescription("proj.scale", "Projection", "Scaling factor for input coordinates");
256 
257  oc.doRegister("proj.rotate", new Option_Float(0.0));
258  oc.addDescription("proj.rotate", "Projection", "Rotation (clockwise degrees) for input coordinates");
259 
260 #ifdef PROJ_API_FILE
261  oc.doRegister("proj.utm", new Option_Bool(false));
262  oc.addDescription("proj.utm", "Projection", "Determine the UTM zone (for a universal transversal mercator projection based on the WGS84 ellipsoid)");
263 
264  oc.doRegister("proj.dhdn", new Option_Bool(false));
265  oc.addDescription("proj.dhdn", "Projection", "Determine the DHDN zone (for a transversal mercator projection based on the bessel ellipsoid, \"Gauss-Krueger\")");
266 
267  oc.doRegister("proj", new Option_String("!"));
268  oc.addDescription("proj", "Projection", "Uses STR as proj.4 definition for projection");
269 
270  oc.doRegister("proj.inverse", new Option_Bool(false));
271  oc.addDescription("proj.inverse", "Projection", "Inverses projection");
272 
273  oc.doRegister("proj.dhdnutm", new Option_Bool(false));
274  oc.addDescription("proj.dhdnutm", "Projection", "Convert from Gauss-Krueger to UTM");
275 #endif // PROJ_API_FILE
276 }
277 
278 
279 bool
281  return myProjectionMethod != NONE;
282 }
283 
284 
285 bool
287  return myUseInverseProjection;
288 }
289 
290 
291 void
293  cartesian.sub(getOffsetBase());
294  if (myProjectionMethod == NONE) {
295  return;
296  }
297  if (myProjectionMethod == SIMPLE) {
298  const double y = cartesian.y() / 111136.;
299  const double x = cartesian.x() / 111320. / cos(DEG2RAD(y));
300  cartesian.set(x, y);
301  return;
302  }
303 #ifdef PROJ_API_FILE
304 #ifdef PROJ_VERSION_MAJOR
305  PJ_COORD c;
306  c.xy.x = cartesian.x();
307  c.xy.y = cartesian.y();
308  c = proj_trans(myProjection, PJ_INV, c);
309  cartesian.set(proj_todeg(c.lp.lam), proj_todeg(c.lp.phi));
310 #else
311  projUV p;
312  p.u = cartesian.x();
313  p.v = cartesian.y();
314  p = pj_inv(p, myProjection);
316  p.u *= RAD_TO_DEG;
317  p.v *= RAD_TO_DEG;
318  cartesian.set((double) p.u, (double) p.v);
319 #endif
320 #endif
321 }
322 
323 
324 bool
325 GeoConvHelper::x2cartesian(Position& from, bool includeInBoundary) {
326  if (includeInBoundary) {
327  myOrigBoundary.add(from);
328  }
329  // init projection parameter on first use
330 #ifdef PROJ_API_FILE
331  if (myProjection == nullptr) {
332  double x = from.x() * myGeoScale;
333  switch (myProjectionMethod) {
334  case DHDN_UTM: {
335  int zone = (int)((x - 500000.) / 1000000.);
336  if (zone < 1 || zone > 5) {
337  WRITE_WARNING("Attempt to initialize DHDN_UTM-projection on invalid longitude " + toString(x));
338  return false;
339  }
340  myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
341  " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
342  " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
343 #ifdef PROJ_VERSION_MAJOR
344  myInverseProjection = proj_create(PJ_DEFAULT_CTX, myProjString.c_str());
345  myGeoProjection = proj_create(PJ_DEFAULT_CTX, "+proj=latlong +datum=WGS84");
346 #else
347  myInverseProjection = pj_init_plus(myProjString.c_str());
348  myGeoProjection = pj_init_plus("+proj=latlong +datum=WGS84");
349 #endif
351  x = ((x - 500000.) / 1000000.) * 3; // continues with UTM
352  }
353  FALLTHROUGH;
354  case UTM: {
355  int zone = (int)(x + 180) / 6 + 1;
356  myProjString = "+proj=utm +zone=" + toString(zone) +
357  " +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
358 #ifdef PROJ_VERSION_MAJOR
359  myProjection = proj_create(PJ_DEFAULT_CTX, myProjString.c_str());
360 #else
361  myProjection = pj_init_plus(myProjString.c_str());
362 #endif
364  }
365  break;
366  case DHDN: {
367  int zone = (int)(x / 3);
368  if (zone < 1 || zone > 5) {
369  WRITE_WARNING("Attempt to initialize DHDN-projection on invalid longitude " + toString(x));
370  return false;
371  }
372  myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
373  " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
374  " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
375 #ifdef PROJ_VERSION_MAJOR
376  myProjection = proj_create(PJ_DEFAULT_CTX, myProjString.c_str());
377 #else
378  myProjection = pj_init_plus(myProjString.c_str());
379 #endif
381  }
382  break;
383  default:
384  break;
385  }
386  }
387  if (myInverseProjection != nullptr) {
388 #ifdef PROJ_VERSION_MAJOR
389  PJ_COORD c;
390  c.xy.x = from.x();
391  c.xy.y = from.y();
392  c = proj_trans(myInverseProjection, PJ_INV, c);
393  from.set(proj_todeg(c.lp.lam), proj_todeg(c.lp.phi));
394 #else
395  double x = from.x();
396  double y = from.y();
397  if (pj_transform(myInverseProjection, myGeoProjection, 1, 1, &x, &y, nullptr)) {
398  WRITE_WARNING("Could not transform (" + toString(x) + "," + toString(y) + ")");
399  }
400  from.set(double(x * RAD_TO_DEG), double(y * RAD_TO_DEG));
401 #endif
402  }
403 #endif
404  // perform conversion
405  bool ok = x2cartesian_const(from);
406  if (ok) {
407  if (includeInBoundary) {
408  myConvBoundary.add(from);
409  }
410  }
411  return ok;
412 }
413 
414 
415 bool
417  double x2 = from.x() * myGeoScale;
418  double y2 = from.y() * myGeoScale;
419  double x = x2 * myCos - y2 * mySin;
420  double y = x2 * mySin + y2 * myCos;
421  if (myProjectionMethod == NONE) {
422  // do nothing
423  } else if (myUseInverseProjection) {
424  cartesian2geo(from);
425  } else {
426  if (x > 180.1 || x < -180.1) {
427  WRITE_WARNING("Invalid longitude " + toString(x));
428  return false;
429  }
430  if (y > 90.1 || y < -90.1) {
431  WRITE_WARNING("Invalid latitude " + toString(y));
432  return false;
433  }
434 #ifdef PROJ_API_FILE
435  if (myProjection != nullptr) {
436 #ifdef PROJ_VERSION_MAJOR
437  PJ_COORD c;
438  c.lp.lam = proj_torad(x);
439  c.lp.phi = proj_torad(y);
440  c = proj_trans(myProjection, PJ_FWD, c);
442  x = c.xy.x;
443  y = c.xy.y;
444 #else
445  projUV p;
446  p.u = x * DEG_TO_RAD;
447  p.v = y * DEG_TO_RAD;
448  p = pj_fwd(p, myProjection);
450  x = p.u;
451  y = p.v;
452 #endif
453  }
454 #endif
455  if (myProjectionMethod == SIMPLE) {
456  x *= 111320. * cos(DEG2RAD(y));
457  y *= 111136.;
459  }
460  }
461  if (x > std::numeric_limits<double>::max() ||
462  y > std::numeric_limits<double>::max()) {
463  return false;
464  }
465  from.set(x, y);
466  from.add(myOffset);
467  if (myFlatten) {
468  from.setz(0);
469  }
470  return true;
471 }
472 
473 
474 void
475 GeoConvHelper::moveConvertedBy(double x, double y) {
476  myOffset.add(x, y);
477  myConvBoundary.moveby(x, y);
478 }
479 
480 
481 const Boundary&
483  return myOrigBoundary;
484 }
485 
486 
487 const Boundary&
489  return myConvBoundary;
490 }
491 
492 
493 const Position
495  return myOffset;
496 }
497 
498 
499 const Position
501  return myOffset;
502 }
503 
504 
505 const std::string&
507  return myProjString;
508 }
509 
510 
511 void
513  if (myNumLoaded == 0) {
515  if (lefthand) {
516  myFinal.myOffset.mul(1, -1);
517  }
518  } else {
519  if (lefthand) {
520  myProcessing.myOffset.mul(1, -1);
521  }
523  // prefer options over loaded location
525  // let offset and boundary lead back to the original coords of the loaded data
528  // the new boundary (updated during loading)
530  }
531  if (lefthand) {
533  }
534 }
535 
536 
537 void
539  myNumLoaded++;
540  if (myNumLoaded > 1) {
541  WRITE_WARNING("Ignoring loaded location attribute nr. " + toString(myNumLoaded) + " for tracking of original location");
542  } else {
543  myLoaded = loaded;
544  }
545 }
546 
547 
548 void
550  myNumLoaded = 0;
551 }
552 
553 
554 void
559  if (myFinal.usingGeoProjection()) {
561  }
563  if (myFinal.usingGeoProjection()) {
564  into.setPrecision();
565  }
567  into.closeTag();
568  into.lf();
569 }
570 
571 
572 /****************************************************************************/
#define DEG2RAD(x)
Definition: GeomHelper.h:35
#define WRITE_ERROR(msg)
Definition: MsgHandler.h:284
#define WRITE_WARNING(msg)
Definition: MsgHandler.h:276
@ 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:34
std::string toString(const T &t, std::streamsize accuracy=gPrecision)
Definition: ToString.h:44
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:367
void flipY()
flips ymin and ymax
Definition: Boundary.cpp:321
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:60
void lf()
writes a line feed if applicable
Definition: OutputDevice.h:227
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:239
bool closeTag(const std::string &comment="")
Closes the most recently opened tag and optionally adds a comment.
void setPrecision(int precision=gPrecision)
Sets the precison or resets it to default.
A point in 2D or 3D with translation and scaling methods.
Definition: Position.h:36
void set(double x, double y)
set positions x and y
Definition: Position.h:84
void sub(double dx, double dy)
Substracts the given position from this one.
Definition: Position.h:144
double x() const
Returns the x-position.
Definition: Position.h:54
void add(const Position &pos)
Adds the given position to this one.
Definition: Position.h:124
void setz(double z)
set position z
Definition: Position.h:79
void mul(double val)
Multiplies both positions with the given value.
Definition: Position.h:104
double y() const
Returns the y-position.
Definition: Position.h:59