Eclipse SUMO - Simulation of Urban MObility
Loading...
Searching...
No Matches
GeoConvHelper.cpp
Go to the documentation of this file.
1/****************************************************************************/
2// Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.dev/sumo
3// Copyright (C) 2001-2023 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
45std::map<std::string, std::pair<std::string, Position> > GeoConvHelper::myLoadedPlain;
46
47// ===========================================================================
48// method definitions
49// ===========================================================================
50
51GeoConvHelper::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 PROJ_API_FILE
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 PROJ_API_FILE
79 } else {
81 initProj(proj);
82 if (myProjection == nullptr) {
83 // avoid error about missing datum shift file
84 myProjString = std::regex_replace(proj, std::regex("\\+geoidgrids[^ ]*"), std::string(""));
85 myProjString = std::regex_replace(myProjString, std::regex("\\+step \\+proj=vgridshift \\+grids[^ ]*"), std::string(""));
86 if (myProjString != proj) {
87 WRITE_WARNING(TL("Ignoring geoidgrids and vgridshift in projection"));
88 initProj(myProjString);
89 }
90 }
91 if (myProjection == nullptr) {
92 // !!! check pj_errno
93 throw ProcessError(TL("Could not build projection!"));
94 }
95#endif
96 }
97}
98
99
100#ifdef PROJ_API_FILE
101void
102GeoConvHelper::initProj(const std::string& proj) {
103#ifdef PROJ_VERSION_MAJOR
104 myProjection = proj_create(PJ_DEFAULT_CTX, proj.c_str());
105#else
106 myProjection = pj_init_plus(proj.c_str());
107#endif
108}
109#endif
110
111
113#ifdef PROJ_API_FILE
114 if (myProjection != nullptr) {
115#ifdef PROJ_VERSION_MAJOR
116 proj_destroy(myProjection);
117#else
118 pj_free(myProjection);
119#endif
120 }
121 if (myInverseProjection != nullptr) {
122#ifdef PROJ_VERSION_MAJOR
123 proj_destroy(myInverseProjection);
124#else
125 pj_free(myInverseProjection);
126#endif
127 }
128 if (myGeoProjection != nullptr) {
129#ifdef PROJ_VERSION_MAJOR
130 proj_destroy(myGeoProjection);
131#else
132 pj_free(myGeoProjection);
133#endif
134 }
135#endif
136}
137
138bool
153
157 myOffset = orig.myOffset;
161 myGeoScale = orig.myGeoScale;
162 myCos = orig.myCos;
163 mySin = orig.mySin;
165 myFlatten = orig.myFlatten;
166#ifdef PROJ_API_FILE
167 if (myProjection != nullptr) {
168#ifdef PROJ_VERSION_MAJOR
169 proj_destroy(myProjection);
170#else
171 pj_free(myProjection);
172#endif
173 myProjection = nullptr;
174 }
175 if (myInverseProjection != nullptr) {
176#ifdef PROJ_VERSION_MAJOR
177 proj_destroy(myInverseProjection);
178#else
179 pj_free(myInverseProjection);
180#endif
181 myInverseProjection = nullptr;
182 }
183 if (myGeoProjection != nullptr) {
184#ifdef PROJ_VERSION_MAJOR
185 proj_destroy(myGeoProjection);
186#else
187 pj_free(myGeoProjection);
188#endif
189 myGeoProjection = nullptr;
190 }
191 if (orig.myProjection != nullptr) {
192#ifdef PROJ_VERSION_MAJOR
193 myProjection = proj_create(PJ_DEFAULT_CTX, orig.myProjString.c_str());
194#else
195 myProjection = pj_init_plus(orig.myProjString.c_str());
196#endif
197 }
198 if (orig.myInverseProjection != nullptr) {
199#ifdef PROJ_VERSION_MAJOR
200 myInverseProjection = orig.myInverseProjection;
201#else
202 myInverseProjection = pj_init_plus(pj_get_def(orig.myInverseProjection, 0));
203#endif
204 }
205 if (orig.myGeoProjection != nullptr) {
206#ifdef PROJ_VERSION_MAJOR
207 myGeoProjection = orig.myGeoProjection;
208#else
209 myGeoProjection = pj_init_plus(pj_get_def(orig.myGeoProjection, 0));
210#endif
211 }
212#endif
213 return *this;
214}
215
216
217bool
219 std::string proj = "!"; // the default
220 double scale = oc.getFloat("proj.scale");
221 double rot = oc.getFloat("proj.rotate");
222 Position offset = Position(oc.getFloat("offset.x"), oc.getFloat("offset.y"), oc.getFloat("offset.z"));
223 bool inverse = oc.exists("proj.inverse") && oc.getBool("proj.inverse");
224 bool flatten = oc.exists("flatten") && oc.getBool("flatten");
225
226 if (oc.getBool("simple-projection")) {
227 proj = "-";
228 }
229
230#ifdef PROJ_API_FILE
231 if (oc.getBool("proj.inverse") && oc.getString("proj") == "!") {
232 WRITE_ERROR(TL("Inverse projection works only with explicit proj parameters."));
233 return false;
234 }
235 unsigned numProjections = oc.getBool("simple-projection") + oc.getBool("proj.utm") + oc.getBool("proj.dhdn") + oc.getBool("proj.dhdnutm") + (oc.getString("proj").length() > 1);
236 if (numProjections > 1) {
237 WRITE_ERROR(TL("The projection method needs to be uniquely defined."));
238 return false;
239 }
240
241 if (oc.getBool("proj.utm")) {
242 proj = "UTM";
243 } else if (oc.getBool("proj.dhdn")) {
244 proj = "DHDN";
245 } else if (oc.getBool("proj.dhdnutm")) {
246 proj = "DHDN_UTM";
247 } else if (!oc.isDefault("proj")) {
248 proj = oc.getString("proj");
249 }
250#endif
251 myProcessing = GeoConvHelper(proj, offset, Boundary(), Boundary(), scale, rot, inverse, flatten);
253 return true;
254}
255
256
257void
258GeoConvHelper::init(const std::string& proj, const Position& offset, const Boundary& orig,
259 const Boundary& conv, double scale) {
260 myProcessing = GeoConvHelper(proj, offset, orig, conv, scale);
263}
264
265void
267#ifdef PROJ_API_FILE
268 if (myProjection == nullptr &&
270 const std::string origProj = myProjString;
271 // try to initialized projection based on origBoundary
273 x2cartesian(tmp, false);
274 if (myProjection == nullptr) {
275 WRITE_WARNING("Failed to intialized projection '" + origProj + "' based on origBoundary centered on '" + toString(myOrigBoundary.getCenter()) + "'");
277 }
278 }
279#endif
280}
281
282void
284 oc.addOptionSubTopic("Projection");
285
286 oc.doRegister("simple-projection", new Option_Bool(false));
287 oc.addSynonyme("simple-projection", "proj.simple", true);
288 oc.addDescription("simple-projection", "Projection", TL("Uses a simple method for projection"));
289
290 oc.doRegister("proj.scale", new Option_Float(1.0));
291 oc.addDescription("proj.scale", "Projection", TL("Scaling factor for input coordinates"));
292
293 oc.doRegister("proj.rotate", new Option_Float(0.0));
294 oc.addDescription("proj.rotate", "Projection", TL("Rotation (clockwise degrees) for input coordinates"));
295
296#ifdef PROJ_API_FILE
297 oc.doRegister("proj.utm", new Option_Bool(false));
298 oc.addDescription("proj.utm", "Projection", TL("Determine the UTM zone (for a universal transversal mercator projection based on the WGS84 ellipsoid)"));
299
300 oc.doRegister("proj.dhdn", new Option_Bool(false));
301 oc.addDescription("proj.dhdn", "Projection", "Determine the DHDN zone (for a transversal mercator projection based on the bessel ellipsoid, \"Gauss-Krueger\")");
302
303 oc.doRegister("proj", new Option_String("!"));
304 oc.addDescription("proj", "Projection", TL("Uses STR as proj.4 definition for projection"));
305
306 oc.doRegister("proj.inverse", new Option_Bool(false));
307 oc.addDescription("proj.inverse", "Projection", TL("Inverses projection"));
308
309 oc.doRegister("proj.dhdnutm", new Option_Bool(false));
310 oc.addDescription("proj.dhdnutm", "Projection", TL("Convert from Gauss-Krueger to UTM"));
311#endif // PROJ_API_FILE
312}
313
314
315bool
319
320
321bool
325
326
327void
329 cartesian.sub(getOffsetBase());
330 if (myProjectionMethod == NONE) {
331 return;
332 }
333 if (myProjectionMethod == SIMPLE) {
334 const double y = cartesian.y() / 111136.;
335 const double x = cartesian.x() / 111320. / cos(DEG2RAD(y));
336 cartesian.set(x, y);
337 return;
338 }
339#ifdef PROJ_API_FILE
340#ifdef PROJ_VERSION_MAJOR
341 PJ_COORD c;
342 c.xy.x = cartesian.x();
343 c.xy.y = cartesian.y();
344 c = proj_trans(myProjection, PJ_INV, c);
345 cartesian.set(proj_todeg(c.lp.lam), proj_todeg(c.lp.phi));
346#else
347 projUV p;
348 p.u = cartesian.x();
349 p.v = cartesian.y();
350 p = pj_inv(p, myProjection);
352 p.u *= RAD_TO_DEG;
353 p.v *= RAD_TO_DEG;
354 cartesian.set((double) p.u, (double) p.v);
355#endif
356#endif
357}
358
359
360bool
361GeoConvHelper::x2cartesian(Position& from, bool includeInBoundary) {
362 if (includeInBoundary) {
363 myOrigBoundary.add(from);
364 }
365 // init projection parameter on first use
366#ifdef PROJ_API_FILE
367 if (myProjection == nullptr) {
368 double x = from.x() * myGeoScale;
369 switch (myProjectionMethod) {
370 case DHDN_UTM: {
371 int zone = (int)((x - 500000.) / 1000000.);
372 if (zone < 1 || zone > 5) {
373 WRITE_WARNING("Attempt to initialize DHDN_UTM-projection on invalid longitude " + toString(x));
374 return false;
375 }
376 myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
377 " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
378 " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
379#ifdef PROJ_VERSION_MAJOR
380 myInverseProjection = proj_create(PJ_DEFAULT_CTX, myProjString.c_str());
381 myGeoProjection = proj_create(PJ_DEFAULT_CTX, "+proj=latlong +datum=WGS84");
382#else
383 myInverseProjection = pj_init_plus(myProjString.c_str());
384 myGeoProjection = pj_init_plus("+proj=latlong +datum=WGS84");
385#endif
387 x = ((x - 500000.) / 1000000.) * 3; // continues with UTM
388 }
390 case UTM: {
391 int zone = (int)(x + 180) / 6 + 1;
392 myProjString = "+proj=utm +zone=" + toString(zone) +
393 " +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
394#ifdef PROJ_VERSION_MAJOR
395 myProjection = proj_create(PJ_DEFAULT_CTX, myProjString.c_str());
396#else
397 myProjection = pj_init_plus(myProjString.c_str());
398#endif
400 }
401 break;
402 case DHDN: {
403 int zone = (int)(x / 3);
404 if (zone < 1 || zone > 5) {
405 WRITE_WARNING("Attempt to initialize DHDN-projection on invalid longitude " + toString(x));
406 return false;
407 }
408 myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
409 " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
410 " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
411#ifdef PROJ_VERSION_MAJOR
412 myProjection = proj_create(PJ_DEFAULT_CTX, myProjString.c_str());
413#else
414 myProjection = pj_init_plus(myProjString.c_str());
415#endif
417 }
418 break;
419 default:
420 break;
421 }
422 }
423 if (myInverseProjection != nullptr) {
424#ifdef PROJ_VERSION_MAJOR
425 PJ_COORD c;
426 c.xy.x = from.x();
427 c.xy.y = from.y();
428 c = proj_trans(myInverseProjection, PJ_INV, c);
429 from.set(proj_todeg(c.lp.lam), proj_todeg(c.lp.phi));
430#else
431 double x = from.x();
432 double y = from.y();
433 if (pj_transform(myInverseProjection, myGeoProjection, 1, 1, &x, &y, nullptr)) {
434 WRITE_WARNINGF(TL("Could not transform (%,%)"), toString(x), toString(y));
435 }
436 from.set(double(x * RAD_TO_DEG), double(y * RAD_TO_DEG));
437#endif
438 }
439#endif
440 // perform conversion
441 bool ok = x2cartesian_const(from);
442 if (ok) {
443 if (includeInBoundary) {
444 myConvBoundary.add(from);
445 }
446 }
447 return ok;
448}
449
450
451bool
453 double x2 = from.x() * myGeoScale;
454 double y2 = from.y() * myGeoScale;
455 double x = x2 * myCos - y2 * mySin;
456 double y = x2 * mySin + y2 * myCos;
457 if (myProjectionMethod == NONE) {
458 // do nothing
459 } else if (myUseInverseProjection) {
460 cartesian2geo(from);
461 } else {
462 if (x > 180.1 || x < -180.1) {
463 WRITE_WARNING("Invalid longitude " + toString(x));
464 return false;
465 }
466 if (y > 90.1 || y < -90.1) {
467 WRITE_WARNING("Invalid latitude " + toString(y));
468 return false;
469 }
470#ifdef PROJ_API_FILE
471 if (myProjection != nullptr) {
472#ifdef PROJ_VERSION_MAJOR
473 PJ_COORD c;
474 c.lp.lam = proj_torad(x);
475 c.lp.phi = proj_torad(y);
476 c = proj_trans(myProjection, PJ_FWD, c);
478 x = c.xy.x;
479 y = c.xy.y;
480#else
481 projUV p;
482 p.u = x * DEG_TO_RAD;
483 p.v = y * DEG_TO_RAD;
484 p = pj_fwd(p, myProjection);
486 x = p.u;
487 y = p.v;
488#endif
489 }
490#endif
491 if (myProjectionMethod == SIMPLE) {
492 // Sinusoidal projection (https://en.wikipedia.org/wiki/Sinusoidal_projection)
493 x *= 111320. * cos(DEG2RAD(y));
494 y *= 111136.;
496 }
497 }
498 if (x > std::numeric_limits<double>::max() ||
499 y > std::numeric_limits<double>::max()) {
500 return false;
501 }
502 from.set(x, y);
503 from.add(myOffset);
504 if (myFlatten) {
505 from.setz(0);
506 }
507 return true;
508}
509
510
511void
513 myOffset.add(x, y);
515}
516
517
518const Boundary&
522
523
524const Boundary&
528
529
530const Position
532 return myOffset;
533}
534
535
536const Position
538 return myOffset;
539}
540
541
542const std::string&
546
547
548void
550 if (myNumLoaded == 0) {
552 if (lefthand) {
553 myFinal.myOffset.mul(1, -1);
554 }
555 } else {
556 if (lefthand) {
558 }
560 // prefer options over loaded location
562 // let offset and boundary lead back to the original coords of the loaded data
565 // the new boundary (updated during loading)
567 }
568 if (lefthand) {
570 }
571}
572
573
574void
576 myNumLoaded++;
577 if (myNumLoaded > 1) {
578 WRITE_WARNINGF(TL("Ignoring loaded location attribute nr. % for tracking of original location"), toString(myNumLoaded));
579 } else {
580 myLoaded = loaded;
581 }
582}
583
584
585void
586GeoConvHelper::setLoadedPlain(const std::string& nodFile, const GeoConvHelper& loaded) {
587 myLoadedPlain[nodFile] = {loaded.getProjString(), loaded.getOffset()};
588}
589
590
592GeoConvHelper::getLoadedPlain(const std::string& edgFile) {
593 std::string nodFile = StringUtils::replace(edgFile, ".edg.xml", ".nod.xml");
594 auto it = myLoadedPlain.find(nodFile);
595 if (it != myLoadedPlain.end()) {
596 return new GeoConvHelper(it->second.first, it->second.second, Boundary(), Boundary());
597 } else {
598 return nullptr;
599 }
600}
601
602
603void
608
609
610void
626
627
628/****************************************************************************/
#define DEG2RAD(x)
Definition GeomHelper.h:35
#define WRITE_WARNINGF(...)
Definition MsgHandler.h:271
#define WRITE_ERROR(msg)
Definition MsgHandler.h:279
#define WRITE_WARNING(msg)
Definition MsgHandler.h:270
#define TL(string)
Definition MsgHandler.h:287
@ SUMO_TAG_LOCATION
@ SUMO_ATTR_CONV_BOUNDARY
@ SUMO_ATTR_NET_OFFSET
@ SUMO_ATTR_ORIG_BOUNDARY
@ SUMO_ATTR_ORIG_PROJ
int gPrecisionGeo
Definition StdDefs.cpp:27
#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
Position getCenter() const
Returns the center of the boundary.
Definition Boundary.cpp:112
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
static void resetLoaded()
resets loaded location elements
static void setLoadedPlain(const std::string &nodFile, const GeoConvHelper &loaded)
registers the coordinate transformation as having been loaded from the given file
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.
static GeoConvHelper * getLoadedPlain(const std::string &edgFile)
bool usingInverseGeoProjection() const
Returns the information whether an inverse transformation will happen.
bool operator==(const GeoConvHelper &o) const
void resolveAbstractProjection()
init projString such as 'UTM' in loaded projection
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 std::map< std::string, std::pair< std::string, Position > > myLoadedPlain
the projections loaded from .nod.xml (to be re-used when loading edg.xml)
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.
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)
bool isDefault(const std::string &name) const
Returns the information whether the named option has still the default value.
void doRegister(const std::string &name, Option *o)
Adds an option under the given name.
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.
void lf()
writes a line feed if applicable
OutputDevice & writeAttr(const SumoXMLAttr attr, const T &val)
writes a named attribute
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)
Subtracts 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
static std::string replace(std::string str, const std::string &what, const std::string &by)
Replaces all occurrences of the second string by the third string within the first string.