1 /****************************************************************************/
2 // Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.org/sumo
3 // Copyright (C) 2001-2019 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 /****************************************************************************/
10 /// @file GeoConvHelper.cpp
11 /// @author Daniel Krajzewicz
12 /// @author Jakob Erdmann
13 /// @author Michael Behrisch
14 /// @date 2006-08-01
15 /// @version $Id$
16 ///
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>
30 #include <utils/common/MsgHandler.h>
31 #include <utils/common/ToString.h>
32 #include <utils/geom/GeomHelper.h>
33 #include <utils/options/OptionsCont.h>
34 #include <utils/iodevices/OutputDevice.h>
35 #include "GeoConvHelper.h"
36
37
38 // ===========================================================================
39 // static member variables
40 // ===========================================================================
41
42 GeoConvHelper GeoConvHelper::myProcessing("!", Position(), Boundary(), Boundary());
43 GeoConvHelper GeoConvHelper::myLoaded("!", Position(), Boundary(), Boundary());
44 GeoConvHelper GeoConvHelper::myFinal("!", Position(), Boundary(), Boundary());
45 int GeoConvHelper::myNumLoaded = 0;
46
47 // ===========================================================================
48 // method definitions
49 // ===========================================================================
50
GeoConvHelper(const std::string & proj,const Position & offset,const Boundary & orig,const Boundary & conv,double scale,double rot,bool inverse,bool flatten)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 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 == "!") {
69 myProjectionMethod = NONE;
70 } else if (proj == "-") {
71 myProjectionMethod = SIMPLE;
72 } else if (proj == "UTM") {
73 myProjectionMethod = UTM;
74 } else if (proj == "DHDN") {
75 myProjectionMethod = DHDN;
76 } else if (proj == "DHDN_UTM") {
77 myProjectionMethod = DHDN_UTM;
78 #ifdef PROJ_API_FILE
79 } else {
80 myProjectionMethod = PROJ;
81 #ifdef PROJ_VERSION_MAJOR
82 myProjection = proj_create(PJ_DEFAULT_CTX, proj.c_str());
83 #else
84 myProjection = pj_init_plus(proj.c_str());
85 #endif
86 if (myProjection == nullptr) {
87 // !!! check pj_errno
88 throw ProcessError("Could not build projection!");
89 }
90 #endif
91 }
92 }
93
94
~GeoConvHelper()95 GeoConvHelper::~GeoConvHelper() {
96 #ifdef PROJ_API_FILE
97 if (myProjection != nullptr) {
98 #ifdef PROJ_VERSION_MAJOR
99 proj_destroy(myProjection);
100 #else
101 pj_free(myProjection);
102 #endif
103 }
104 if (myInverseProjection != nullptr) {
105 #ifdef PROJ_VERSION_MAJOR
106 proj_destroy(myInverseProjection);
107 #else
108 pj_free(myInverseProjection);
109 #endif
110 }
111 if (myGeoProjection != nullptr) {
112 #ifdef PROJ_VERSION_MAJOR
113 proj_destroy(myGeoProjection);
114 #else
115 pj_free(myGeoProjection);
116 #endif
117 }
118 #endif
119 }
120
121 bool
operator ==(const GeoConvHelper & o) const122 GeoConvHelper::operator==(const GeoConvHelper& o) const {
123 return (
124 myProjString == o.myProjString &&
125 myOffset == o.myOffset &&
126 myProjectionMethod == o.myProjectionMethod &&
127 myOrigBoundary == o.myOrigBoundary &&
128 myConvBoundary == o.myConvBoundary &&
129 myGeoScale == o.myGeoScale &&
130 myCos == o.myCos &&
131 mySin == o.mySin &&
132 myUseInverseProjection == o.myUseInverseProjection &&
133 myFlatten == o.myFlatten
134 );
135 }
136
137 GeoConvHelper&
operator =(const GeoConvHelper & orig)138 GeoConvHelper::operator=(const GeoConvHelper& orig) {
139 myProjString = orig.myProjString;
140 myOffset = orig.myOffset;
141 myProjectionMethod = orig.myProjectionMethod;
142 myOrigBoundary = orig.myOrigBoundary;
143 myConvBoundary = orig.myConvBoundary;
144 myGeoScale = orig.myGeoScale;
145 myCos = orig.myCos;
146 mySin = orig.mySin;
147 myUseInverseProjection = orig.myUseInverseProjection;
148 myFlatten = orig.myFlatten;
149 #ifdef PROJ_API_FILE
150 if (myProjection != nullptr) {
151 #ifdef PROJ_VERSION_MAJOR
152 proj_destroy(myProjection);
153 #else
154 pj_free(myProjection);
155 #endif
156 myProjection = nullptr;
157 }
158 if (myInverseProjection != nullptr) {
159 #ifdef PROJ_VERSION_MAJOR
160 proj_destroy(myInverseProjection);
161 #else
162 pj_free(myInverseProjection);
163 #endif
164 myInverseProjection = nullptr;
165 }
166 if (myGeoProjection != nullptr) {
167 #ifdef PROJ_VERSION_MAJOR
168 proj_destroy(myGeoProjection);
169 #else
170 pj_free(myGeoProjection);
171 #endif
172 myGeoProjection = nullptr;
173 }
174 if (orig.myProjection != nullptr) {
175 #ifdef PROJ_VERSION_MAJOR
176 myProjection = proj_create(PJ_DEFAULT_CTX, orig.myProjString.c_str());
177 #else
178 myProjection = pj_init_plus(orig.myProjString.c_str());
179 #endif
180 }
181 if (orig.myInverseProjection != nullptr) {
182 #ifdef PROJ_VERSION_MAJOR
183 myInverseProjection = orig.myInverseProjection;
184 #else
185 myInverseProjection = pj_init_plus(pj_get_def(orig.myInverseProjection, 0));
186 #endif
187 }
188 if (orig.myGeoProjection != nullptr) {
189 #ifdef PROJ_VERSION_MAJOR
190 myGeoProjection = orig.myGeoProjection;
191 #else
192 myGeoProjection = pj_init_plus(pj_get_def(orig.myGeoProjection, 0));
193 #endif
194 }
195 #endif
196 return *this;
197 }
198
199
200 bool
init(OptionsCont & oc)201 GeoConvHelper::init(OptionsCont& oc) {
202 std::string proj = "!"; // the default
203 double scale = oc.getFloat("proj.scale");
204 double rot = oc.getFloat("proj.rotate");
205 Position offset = Position(oc.getFloat("offset.x"), oc.getFloat("offset.y"));
206 bool inverse = oc.exists("proj.inverse") && oc.getBool("proj.inverse");
207 bool flatten = oc.exists("flatten") && oc.getBool("flatten");
208
209 if (oc.getBool("simple-projection")) {
210 proj = "-";
211 }
212
213 #ifdef PROJ_API_FILE
214 if (oc.getBool("proj.inverse") && oc.getString("proj") == "!") {
215 WRITE_ERROR("Inverse projection works only with explicit proj parameters.");
216 return false;
217 }
218 unsigned numProjections = oc.getBool("simple-projection") + oc.getBool("proj.utm") + oc.getBool("proj.dhdn") + oc.getBool("proj.dhdnutm") + (oc.getString("proj").length() > 1);
219 if (numProjections > 1) {
220 WRITE_ERROR("The projection method needs to be uniquely defined.");
221 return false;
222 }
223
224 if (oc.getBool("proj.utm")) {
225 proj = "UTM";
226 } else if (oc.getBool("proj.dhdn")) {
227 proj = "DHDN";
228 } else if (oc.getBool("proj.dhdnutm")) {
229 proj = "DHDN_UTM";
230 } else if (!oc.isDefault("proj")) {
231 proj = oc.getString("proj");
232 }
233 #endif
234 myProcessing = GeoConvHelper(proj, offset, Boundary(), Boundary(), scale, rot, inverse, flatten);
235 myFinal = myProcessing;
236 return true;
237 }
238
239
240 void
init(const std::string & proj,const Position & offset,const Boundary & orig,const Boundary & conv,double scale)241 GeoConvHelper::init(const std::string& proj, const Position& offset, const Boundary& orig,
242 const Boundary& conv, double scale) {
243 myProcessing = GeoConvHelper(proj, offset, orig, conv, scale);
244 myFinal = myProcessing;
245 }
246
247
248 void
addProjectionOptions(OptionsCont & oc)249 GeoConvHelper::addProjectionOptions(OptionsCont& oc) {
250 oc.addOptionSubTopic("Projection");
251
252 oc.doRegister("simple-projection", new Option_Bool(false));
253 oc.addSynonyme("simple-projection", "proj.simple", true);
254 oc.addDescription("simple-projection", "Projection", "Uses a simple method for projection");
255
256 oc.doRegister("proj.scale", new Option_Float(1.0));
257 oc.addDescription("proj.scale", "Projection", "Scaling factor for input coordinates");
258
259 oc.doRegister("proj.rotate", new Option_Float(0.0));
260 oc.addDescription("proj.rotate", "Projection", "Rotation (clockwise degrees) for input coordinates");
261
262 #ifdef PROJ_API_FILE
263 oc.doRegister("proj.utm", new Option_Bool(false));
264 oc.addDescription("proj.utm", "Projection", "Determine the UTM zone (for a universal transversal mercator projection based on the WGS84 ellipsoid)");
265
266 oc.doRegister("proj.dhdn", new Option_Bool(false));
267 oc.addDescription("proj.dhdn", "Projection", "Determine the DHDN zone (for a transversal mercator projection based on the bessel ellipsoid, \"Gauss-Krueger\")");
268
269 oc.doRegister("proj", new Option_String("!"));
270 oc.addDescription("proj", "Projection", "Uses STR as proj.4 definition for projection");
271
272 oc.doRegister("proj.inverse", new Option_Bool(false));
273 oc.addDescription("proj.inverse", "Projection", "Inverses projection");
274
275 oc.doRegister("proj.dhdnutm", new Option_Bool(false));
276 oc.addDescription("proj.dhdnutm", "Projection", "Convert from Gauss-Krueger to UTM");
277 #endif // PROJ_API_FILE
278 }
279
280
281 bool
usingGeoProjection() const282 GeoConvHelper::usingGeoProjection() const {
283 return myProjectionMethod != NONE;
284 }
285
286
287 bool
usingInverseGeoProjection() const288 GeoConvHelper::usingInverseGeoProjection() const {
289 return myUseInverseProjection;
290 }
291
292
293 void
cartesian2geo(Position & cartesian) const294 GeoConvHelper::cartesian2geo(Position& cartesian) const {
295 cartesian.sub(getOffsetBase());
296 if (myProjectionMethod == NONE) {
297 return;
298 }
299 if (myProjectionMethod == SIMPLE) {
300 const double y = cartesian.y() / 111136.;
301 const double x = cartesian.x() / 111320. / cos(DEG2RAD(y));
302 cartesian.set(x, y);
303 return;
304 }
305 #ifdef PROJ_API_FILE
306 #ifdef PROJ_VERSION_MAJOR
307 PJ_COORD c;
308 c.xy.x = cartesian.x();
309 c.xy.y = cartesian.y();
310 c = proj_trans(myProjection, PJ_INV, c);
311 cartesian.set(proj_todeg(c.lp.lam), proj_todeg(c.lp.phi));
312 #else
313 projUV p;
314 p.u = cartesian.x();
315 p.v = cartesian.y();
316 p = pj_inv(p, myProjection);
317 //!!! check pj_errno
318 p.u *= RAD_TO_DEG;
319 p.v *= RAD_TO_DEG;
320 cartesian.set((double) p.u, (double) p.v);
321 #endif
322 #endif
323 }
324
325
326 bool
x2cartesian(Position & from,bool includeInBoundary)327 GeoConvHelper::x2cartesian(Position& from, bool includeInBoundary) {
328 if (includeInBoundary) {
329 myOrigBoundary.add(from);
330 }
331 // init projection parameter on first use
332 #ifdef PROJ_API_FILE
333 if (myProjection == nullptr) {
334 double x = from.x() * myGeoScale;
335 switch (myProjectionMethod) {
336 case DHDN_UTM: {
337 int zone = (int)((x - 500000.) / 1000000.);
338 if (zone < 1 || zone > 5) {
339 WRITE_WARNING("Attempt to initialize DHDN_UTM-projection on invalid longitude " + toString(x));
340 return false;
341 }
342 myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
343 " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
344 " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
345 #ifdef PROJ_VERSION_MAJOR
346 myInverseProjection = proj_create(PJ_DEFAULT_CTX, myProjString.c_str());
347 myGeoProjection = proj_create(PJ_DEFAULT_CTX, "+proj=latlong +datum=WGS84");
348 #else
349 myInverseProjection = pj_init_plus(myProjString.c_str());
350 myGeoProjection = pj_init_plus("+proj=latlong +datum=WGS84");
351 #endif
352 //!!! check pj_errno
353 x = ((x - 500000.) / 1000000.) * 3; // continues with UTM
354 }
355 FALLTHROUGH;
356 case UTM: {
357 int zone = (int)(x + 180) / 6 + 1;
358 myProjString = "+proj=utm +zone=" + toString(zone) +
359 " +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
360 #ifdef PROJ_VERSION_MAJOR
361 myProjection = proj_create(PJ_DEFAULT_CTX, myProjString.c_str());
362 #else
363 myProjection = pj_init_plus(myProjString.c_str());
364 #endif
365 //!!! check pj_errno
366 }
367 break;
368 case DHDN: {
369 int zone = (int)(x / 3);
370 if (zone < 1 || zone > 5) {
371 WRITE_WARNING("Attempt to initialize DHDN-projection on invalid longitude " + toString(x));
372 return false;
373 }
374 myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
375 " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
376 " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
377 #ifdef PROJ_VERSION_MAJOR
378 myProjection = proj_create(PJ_DEFAULT_CTX, myProjString.c_str());
379 #else
380 myProjection = pj_init_plus(myProjString.c_str());
381 #endif
382 //!!! check pj_errno
383 }
384 break;
385 default:
386 break;
387 }
388 }
389 if (myInverseProjection != nullptr) {
390 #ifdef PROJ_VERSION_MAJOR
391 PJ_COORD c;
392 c.xy.x = from.x();
393 c.xy.y = from.y();
394 c = proj_trans(myInverseProjection, PJ_INV, c);
395 from.set(proj_todeg(c.lp.lam), proj_todeg(c.lp.phi));
396 #else
397 double x = from.x();
398 double y = from.y();
399 if (pj_transform(myInverseProjection, myGeoProjection, 1, 1, &x, &y, nullptr)) {
400 WRITE_WARNING("Could not transform (" + toString(x) + "," + toString(y) + ")");
401 }
402 from.set(double(x * RAD_TO_DEG), double(y * RAD_TO_DEG));
403 #endif
404 }
405 #endif
406 // perform conversion
407 bool ok = x2cartesian_const(from);
408 if (ok) {
409 if (includeInBoundary) {
410 myConvBoundary.add(from);
411 }
412 }
413 return ok;
414 }
415
416
417 bool
x2cartesian_const(Position & from) const418 GeoConvHelper::x2cartesian_const(Position& from) const {
419 double x2 = from.x() * myGeoScale;
420 double y2 = from.y() * myGeoScale;
421 double x = x2 * myCos - y2 * mySin;
422 double y = x2 * mySin + y2 * myCos;
423 if (myProjectionMethod == NONE) {
424 from.add(myOffset);
425 } else if (myUseInverseProjection) {
426 cartesian2geo(from);
427 } else {
428 if (x > 180.1 || x < -180.1) {
429 WRITE_WARNING("Invalid longitude " + toString(x));
430 return false;
431 }
432 if (y > 90.1 || y < -90.1) {
433 WRITE_WARNING("Invalid latitude " + toString(y));
434 return false;
435 }
436 #ifdef PROJ_API_FILE
437 if (myProjection != nullptr) {
438 #ifdef PROJ_VERSION_MAJOR
439 PJ_COORD c;
440 c.lp.lam = proj_torad(x);
441 c.lp.phi = proj_torad(y);
442 c = proj_trans(myProjection, PJ_FWD, c);
443 //!!! check pj_errno
444 x = c.xy.x;
445 y = c.xy.y;
446 #else
447 projUV p;
448 p.u = x * DEG_TO_RAD;
449 p.v = y * DEG_TO_RAD;
450 p = pj_fwd(p, myProjection);
451 //!!! check pj_errno
452 x = p.u;
453 y = p.v;
454 #endif
455 }
456 #endif
457 if (myProjectionMethod == SIMPLE) {
458 x *= 111320. * cos(DEG2RAD(y));
459 y *= 111136.;
460 //!!! recheck whether the axes are mirrored
461 }
462 }
463 if (x > std::numeric_limits<double>::max() ||
464 y > std::numeric_limits<double>::max()) {
465 return false;
466 }
467 from.set(x, y);
468 from.add(myOffset);
469 if (myFlatten) {
470 from.setz(0);
471 }
472 return true;
473 }
474
475
476 void
moveConvertedBy(double x,double y)477 GeoConvHelper::moveConvertedBy(double x, double y) {
478 myOffset.add(x, y);
479 myConvBoundary.moveby(x, y);
480 }
481
482
483 const Boundary&
getOrigBoundary() const484 GeoConvHelper::getOrigBoundary() const {
485 return myOrigBoundary;
486 }
487
488
489 const Boundary&
getConvBoundary() const490 GeoConvHelper::getConvBoundary() const {
491 return myConvBoundary;
492 }
493
494
495 const Position
getOffset() const496 GeoConvHelper::getOffset() const {
497 return myOffset;
498 }
499
500
501 const Position
getOffsetBase() const502 GeoConvHelper::getOffsetBase() const {
503 return myOffset;
504 }
505
506
507 const std::string&
getProjString() const508 GeoConvHelper::getProjString() const {
509 return myProjString;
510 }
511
512
513 void
computeFinal(bool lefthand)514 GeoConvHelper::computeFinal(bool lefthand) {
515 if (myNumLoaded == 0) {
516 myFinal = myProcessing;
517 if (lefthand) {
518 myFinal.myOffset.mul(1, -1);
519 }
520 } else {
521 if (lefthand) {
522 myProcessing.myOffset.mul(1, -1);
523 }
524 myFinal = GeoConvHelper(
525 // prefer options over loaded location
526 myProcessing.usingGeoProjection() ? myProcessing.getProjString() : myLoaded.getProjString(),
527 // let offset and boundary lead back to the original coords of the loaded data
528 myProcessing.getOffset() + myLoaded.getOffset(),
529 myLoaded.getOrigBoundary(),
530 // the new boundary (updated during loading)
531 myProcessing.getConvBoundary());
532 }
533 if (lefthand) {
534 myFinal.myConvBoundary.flipY();
535 }
536 }
537
538
539 void
setLoaded(const GeoConvHelper & loaded)540 GeoConvHelper::setLoaded(const GeoConvHelper& loaded) {
541 myNumLoaded++;
542 if (myNumLoaded > 1) {
543 WRITE_WARNING("Ignoring loaded location attribute nr. " + toString(myNumLoaded) + " for tracking of original location");
544 } else {
545 myLoaded = loaded;
546 }
547 }
548
549
550 void
resetLoaded()551 GeoConvHelper::resetLoaded() {
552 myNumLoaded = 0;
553 }
554
555
556 void
writeLocation(OutputDevice & into)557 GeoConvHelper::writeLocation(OutputDevice& into) {
558 into.openTag(SUMO_TAG_LOCATION);
559 into.writeAttr(SUMO_ATTR_NET_OFFSET, myFinal.getOffsetBase());
560 into.writeAttr(SUMO_ATTR_CONV_BOUNDARY, myFinal.getConvBoundary());
561 if (myFinal.usingGeoProjection()) {
562 into.setPrecision(gPrecisionGeo);
563 }
564 into.writeAttr(SUMO_ATTR_ORIG_BOUNDARY, myFinal.getOrigBoundary());
565 if (myFinal.usingGeoProjection()) {
566 into.setPrecision();
567 }
568 into.writeAttr(SUMO_ATTR_ORIG_PROJ, myFinal.getProjString());
569 into.closeTag();
570 into.lf();
571 }
572
573
574 /****************************************************************************/
575