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>
34#include "GeoConvHelper.h"
35
36
37// ===========================================================================
38// static member variables
39// ===========================================================================
40
45
46// ===========================================================================
47// method definitions
48// ===========================================================================
49
50GeoConvHelper::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
100void
101GeoConvHelper::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
137bool
139 return (
141 myOffset == o.myOffset &&
145 myGeoScale == o.myGeoScale &&
146 myCos == o.myCos &&
147 mySin == o.mySin &&
150 );
151}
152
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
216bool
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
256void
257GeoConvHelper::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
264void
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
297bool
299 return myProjectionMethod != NONE;
300}
301
302
303bool
306}
307
308
309void
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
342bool
343GeoConvHelper::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 }
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
433bool
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
493void
495 myOffset.add(x, y);
497}
498
499
500const Boundary&
502 return myOrigBoundary;
503}
504
505
506const Boundary&
508 return myConvBoundary;
509}
510
511
512const Position
514 return myOffset;
515}
516
517
518const Position
520 return myOffset;
521}
522
523
524const std::string&
526 return myProjString;
527}
528
529
530void
532 if (myNumLoaded == 0) {
534 if (lefthand) {
535 myFinal.myOffset.mul(1, -1);
536 }
537 } else {
538 if (lefthand) {
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
556void
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
567void
569 myNumLoaded = 0;
570}
571
572
573void
580 }
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 & writeAttr(const SumoXMLAttr attr, const T &val)
writes a named attribute
Definition: OutputDevice.h:251
OutputDevice & openTag(const std::string &xmlElement)
Opens an XML tag.
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