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-2018 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 HAVE_PROJ
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 HAVE_PROJ
79  } else {
81  myProjection = pj_init_plus(proj.c_str());
82  if (myProjection == nullptr) {
83  // !!! check pj_errno
84  throw ProcessError("Could not build projection!");
85  }
86 #endif
87  }
88 }
89 
90 
92 #ifdef HAVE_PROJ
93  if (myProjection != nullptr) {
94  pj_free(myProjection);
95  }
96  if (myInverseProjection != nullptr) {
97  pj_free(myInverseProjection);
98  }
99  if (myGeoProjection != nullptr) {
100  pj_free(myInverseProjection);
101  }
102 #endif
103 }
104 
105 bool
107  return (
108  myProjString == o.myProjString &&
109  myOffset == o.myOffset &&
113  myGeoScale == o.myGeoScale &&
114  myCos == o.myCos &&
115  mySin == o.mySin &&
117  myFlatten == o.myFlatten
118  );
119 }
120 
123  myProjString = orig.myProjString;
124  myOffset = orig.myOffset;
128  myGeoScale = orig.myGeoScale;
129  myCos = orig.myCos;
130  mySin = orig.mySin;
132  myFlatten = orig.myFlatten;
133 #ifdef HAVE_PROJ
134  if (myProjection != nullptr) {
135  pj_free(myProjection);
136  myProjection = nullptr;
137  }
138  if (myInverseProjection != nullptr) {
139  pj_free(myInverseProjection);
140  myInverseProjection = nullptr;
141  }
142  if (myGeoProjection != nullptr) {
143  pj_free(myGeoProjection);
144  myGeoProjection = nullptr;
145  }
146  if (orig.myProjection != nullptr) {
147  myProjection = pj_init_plus(orig.myProjString.c_str());
148  }
149  if (orig.myInverseProjection != nullptr) {
150  myInverseProjection = pj_init_plus(pj_get_def(orig.myInverseProjection, 0));
151  }
152  if (orig.myGeoProjection != nullptr) {
153  myGeoProjection = pj_init_plus(pj_get_def(orig.myGeoProjection, 0));
154  }
155 #endif
156  return *this;
157 }
158 
159 
160 bool
162  std::string proj = "!"; // the default
163  double scale = oc.getFloat("proj.scale");
164  double rot = oc.getFloat("proj.rotate");
165  Position offset = Position(oc.getFloat("offset.x"), oc.getFloat("offset.y"));
166  bool inverse = oc.exists("proj.inverse") && oc.getBool("proj.inverse");
167  bool flatten = oc.exists("flatten") && oc.getBool("flatten");
168 
169  if (oc.getBool("simple-projection")) {
170  proj = "-";
171  }
172 
173 #ifdef HAVE_PROJ
174  if (oc.getBool("proj.inverse") && oc.getString("proj") == "!") {
175  WRITE_ERROR("Inverse projection works only with explicit proj parameters.");
176  return false;
177  }
178  unsigned numProjections = oc.getBool("simple-projection") + oc.getBool("proj.utm") + oc.getBool("proj.dhdn") + oc.getBool("proj.dhdnutm") + (oc.getString("proj").length() > 1);
179  if (numProjections > 1) {
180  WRITE_ERROR("The projection method needs to be uniquely defined.");
181  return false;
182  }
183 
184  if (oc.getBool("proj.utm")) {
185  proj = "UTM";
186  } else if (oc.getBool("proj.dhdn")) {
187  proj = "DHDN";
188  } else if (oc.getBool("proj.dhdnutm")) {
189  proj = "DHDN_UTM";
190  } else if (!oc.isDefault("proj")) {
191  proj = oc.getString("proj");
192  }
193 #endif
194  myProcessing = GeoConvHelper(proj, offset, Boundary(), Boundary(), scale, rot, inverse, flatten);
196  return true;
197 }
198 
199 
200 void
201 GeoConvHelper::init(const std::string& proj, const Position& offset, const Boundary& orig,
202  const Boundary& conv, double scale) {
203  myProcessing = GeoConvHelper(proj, offset, orig, conv, scale);
205 }
206 
207 
208 void
210  oc.addOptionSubTopic("Projection");
211 
212  oc.doRegister("simple-projection", new Option_Bool(false));
213  oc.addSynonyme("simple-projection", "proj.simple", true);
214  oc.addDescription("simple-projection", "Projection", "Uses a simple method for projection");
215 
216  oc.doRegister("proj.scale", new Option_Float(1.0));
217  oc.addDescription("proj.scale", "Projection", "Scaling factor for input coordinates");
218 
219  oc.doRegister("proj.rotate", new Option_Float(0.0));
220  oc.addDescription("proj.rotate", "Projection", "Rotation (clockwise degrees) for input coordinates");
221 
222 #ifdef HAVE_PROJ
223  oc.doRegister("proj.utm", new Option_Bool(false));
224  oc.addDescription("proj.utm", "Projection", "Determine the UTM zone (for a universal transversal mercator projection based on the WGS84 ellipsoid)");
225 
226  oc.doRegister("proj.dhdn", new Option_Bool(false));
227  oc.addDescription("proj.dhdn", "Projection", "Determine the DHDN zone (for a transversal mercator projection based on the bessel ellipsoid, \"Gauss-Krueger\")");
228 
229  oc.doRegister("proj", new Option_String("!"));
230  oc.addDescription("proj", "Projection", "Uses STR as proj.4 definition for projection");
231 
232  oc.doRegister("proj.inverse", new Option_Bool(false));
233  oc.addDescription("proj.inverse", "Projection", "Inverses projection");
234 
235  oc.doRegister("proj.dhdnutm", new Option_Bool(false));
236  oc.addDescription("proj.dhdnutm", "Projection", "Convert from Gauss-Krueger to UTM");
237 #endif // HAVE_PROJ
238 }
239 
240 
241 bool
243  return myProjectionMethod != NONE;
244 }
245 
246 
247 bool
249  return myUseInverseProjection;
250 }
251 
252 
253 void
255  cartesian.sub(getOffsetBase());
256  if (myProjectionMethod == NONE) {
257  return;
258  }
259  if (myProjectionMethod == SIMPLE) {
260  const double y = cartesian.y() / 111136.;
261  const double x = cartesian.x() / 111320. / cos(DEG2RAD(y));
262  cartesian.set(x, y);
263  return;
264  }
265 #ifdef HAVE_PROJ
266  projUV p;
267  p.u = cartesian.x();
268  p.v = cartesian.y();
269  p = pj_inv(p, myProjection);
271  p.u *= RAD_TO_DEG;
272  p.v *= RAD_TO_DEG;
273  cartesian.set((double) p.u, (double) p.v);
274 #endif
275 }
276 
277 
278 bool
279 GeoConvHelper::x2cartesian(Position& from, bool includeInBoundary) {
280  if (includeInBoundary) {
281  myOrigBoundary.add(from);
282  }
283  // init projection parameter on first use
284 #ifdef HAVE_PROJ
285  if (myProjection == nullptr) {
286  double x = from.x() * myGeoScale;
287  switch (myProjectionMethod) {
288  case DHDN_UTM: {
289  int zone = (int)((x - 500000.) / 1000000.);
290  if (zone < 1 || zone > 5) {
291  WRITE_WARNING("Attempt to initialize DHDN_UTM-projection on invalid longitude " + toString(x));
292  return false;
293  }
294  myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
295  " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
296  " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
297  myInverseProjection = pj_init_plus(myProjString.c_str());
298  myGeoProjection = pj_init_plus("+proj=latlong +datum=WGS84");
300  x = ((x - 500000.) / 1000000.) * 3; // continues with UTM
301  }
302  FALLTHROUGH;
303  case UTM: {
304  int zone = (int)(x + 180) / 6 + 1;
305  myProjString = "+proj=utm +zone=" + toString(zone) +
306  " +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
307  myProjection = pj_init_plus(myProjString.c_str());
309  }
310  break;
311  case DHDN: {
312  int zone = (int)(x / 3);
313  if (zone < 1 || zone > 5) {
314  WRITE_WARNING("Attempt to initialize DHDN-projection on invalid longitude " + toString(x));
315  return false;
316  }
317  myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
318  " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
319  " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
320  myProjection = pj_init_plus(myProjString.c_str());
322  }
323  break;
324  default:
325  break;
326  }
327  }
328  if (myInverseProjection != nullptr) {
329  double x = from.x();
330  double y = from.y();
331  if (pj_transform(myInverseProjection, myGeoProjection, 1, 1, &x, &y, nullptr)) {
332  WRITE_WARNING("Could not transform (" + toString(x) + "," + toString(y) + ")");
333  }
334  from.set(double(x * RAD_TO_DEG), double(y * RAD_TO_DEG));
335  }
336 #endif
337  // perform conversion
338  bool ok = x2cartesian_const(from);
339  if (ok) {
340  if (includeInBoundary) {
341  myConvBoundary.add(from);
342  }
343  }
344  return ok;
345 }
346 
347 
348 bool
350  double x2 = from.x() * myGeoScale;
351  double y2 = from.y() * myGeoScale;
352  double x = x2 * myCos - y2 * mySin;
353  double y = x2 * mySin + y2 * myCos;
354  if (myProjectionMethod == NONE) {
355  from.add(myOffset);
356  } else if (myUseInverseProjection) {
357  cartesian2geo(from);
358  } else {
359  if (x > 180.1 || x < -180.1) {
360  WRITE_WARNING("Invalid longitude " + toString(x));
361  return false;
362  }
363  if (y > 90.1 || y < -90.1) {
364  WRITE_WARNING("Invalid latitude " + toString(y));
365  return false;
366  }
367 #ifdef HAVE_PROJ
368  if (myProjection != nullptr) {
369  projUV p;
370  p.u = x * DEG_TO_RAD;
371  p.v = y * DEG_TO_RAD;
372  p = pj_fwd(p, myProjection);
374  x = p.u;
375  y = p.v;
376  }
377 #endif
378  if (myProjectionMethod == SIMPLE) {
379  x *= 111320. * cos(DEG2RAD(y));
380  y *= 111136.;
382  }
383  }
384  if (x > std::numeric_limits<double>::max() ||
385  y > std::numeric_limits<double>::max()) {
386  return false;
387  }
388  from.set(x, y);
389  from.add(myOffset);
390  if (myFlatten) {
391  from.setz(0);
392  }
393  return true;
394 }
395 
396 
397 void
398 GeoConvHelper::moveConvertedBy(double x, double y) {
399  myOffset.add(x, y);
400  myConvBoundary.moveby(x, y);
401 }
402 
403 
404 const Boundary&
406  return myOrigBoundary;
407 }
408 
409 
410 const Boundary&
412  return myConvBoundary;
413 }
414 
415 
416 const Position
418  return myOffset;
419 }
420 
421 
422 const Position
424  return myOffset;
425 }
426 
427 
428 const std::string&
430  return myProjString;
431 }
432 
433 
434 void
436  if (myNumLoaded == 0) {
438  if (lefthand) {
439  myFinal.myOffset.mul(1, -1);
440  }
441  } else {
442  if (lefthand) {
443  myProcessing.myOffset.mul(1, -1);
444  }
446  // prefer options over loaded location
448  // let offset and boundary lead back to the original coords of the loaded data
451  // the new boundary (updated during loading)
453  }
454  if (lefthand) {
456  }
457 }
458 
459 
460 void
462  myNumLoaded++;
463  if (myNumLoaded > 1) {
464  WRITE_WARNING("Ignoring loaded location attribute nr. " + toString(myNumLoaded) + " for tracking of original location");
465  } else {
466  myLoaded = loaded;
467  }
468 }
469 
470 
471 void
473  myNumLoaded = 0;
474 }
475 
476 
477 void
482  if (myFinal.usingGeoProjection()) {
484  }
486  if (myFinal.usingGeoProjection()) {
487  into.setPrecision();
488  }
490  into.closeTag();
491  into.lf();
492 }
493 
494 
495 
496 /****************************************************************************/
497 
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:241
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:53
std::string toString(const T &t, std::streamsize accuracy=gPrecision)
Definition: ToString.h:49
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:247
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:92
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:38