Line data Source code
1 : /****************************************************************************/
2 : // Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.dev/sumo
3 : // Copyright (C) 2001-2026 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 : /****************************************************************************/
14 : /// @file GeoConvHelper.cpp
15 : /// @author Daniel Krajzewicz
16 : /// @author Jakob Erdmann
17 : /// @author Michael Behrisch
18 : /// @date 2006-08-01
19 : ///
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>
29 : #include <utils/common/MsgHandler.h>
30 : #include <utils/common/ToString.h>
31 : #include <utils/geom/GeomHelper.h>
32 : #include <utils/options/OptionsCont.h>
33 : #include <utils/iodevices/OutputDevice.h>
34 : #include "GeoConvHelper.h"
35 :
36 :
37 : // ===========================================================================
38 : // static member variables
39 : // ===========================================================================
40 :
41 : GeoConvHelper GeoConvHelper::myProcessing("!", Position(), Boundary(), Boundary());
42 : GeoConvHelper GeoConvHelper::myLoaded("!", Position(), Boundary(), Boundary());
43 : GeoConvHelper GeoConvHelper::myFinal("!", Position(), Boundary(), Boundary());
44 : int GeoConvHelper::myNumLoaded = 0;
45 : std::map<std::string, std::pair<std::string, Position> > GeoConvHelper::myLoadedPlain;
46 :
47 : // ===========================================================================
48 : // method definitions
49 : // ===========================================================================
50 :
51 206020 : GeoConvHelper::GeoConvHelper(const std::string& proj, const Position& offset,
52 206020 : const Boundary& orig, const Boundary& conv, double scale, double rot, bool inverse, bool flatten):
53 206020 : myProjString(proj),
54 : #ifdef PROJ_API_FILE
55 206020 : myProjection(nullptr),
56 206020 : myInverseProjection(nullptr),
57 206020 : myGeoProjection(nullptr),
58 206020 : myRadians(true),
59 : #endif
60 206020 : myOffset(offset),
61 206020 : myGeoScale(scale),
62 206020 : mySin(sin(DEG2RAD(-rot))), // rotate clockwise
63 206020 : myCos(cos(DEG2RAD(-rot))),
64 206020 : myProjectionMethod(NONE),
65 206020 : myUseInverseProjection(inverse),
66 206020 : myFlatten(flatten),
67 : myOrigBoundary(orig),
68 : myConvBoundary(conv) {
69 : // older PROJ libraries fail to construct inverse projection if this
70 : // string is present
71 618060 : myProjString = StringUtils::replace(myProjString, "+type=crs", "");
72 :
73 206020 : if (proj == "!") {
74 191129 : myProjectionMethod = NONE;
75 14891 : } else if (proj == "-") {
76 3 : myProjectionMethod = SIMPLE;
77 14888 : } else if (proj == "UTM") {
78 607 : myProjectionMethod = UTM;
79 14281 : } else if (proj == "DHDN") {
80 0 : myProjectionMethod = DHDN;
81 14281 : } else if (proj == "DHDN_UTM") {
82 0 : myProjectionMethod = DHDN_UTM;
83 : #ifdef PROJ_API_FILE
84 : } else {
85 14281 : myProjectionMethod = PROJ;
86 14281 : initProj(myProjString);
87 14281 : if (myProjection == nullptr) {
88 : // avoid error about missing datum shift file
89 0 : myProjString = std::regex_replace(proj, std::regex("\\+geoidgrids[^ ]*"), std::string(""));
90 0 : myProjString = std::regex_replace(myProjString, std::regex("\\+step \\+proj=vgridshift \\+grids[^ ]*"), std::string(""));
91 0 : if (myProjString != proj) {
92 0 : WRITE_WARNING(TL("Ignoring geoidgrids and vgridshift in projection"));
93 0 : initProj(myProjString);
94 : }
95 : }
96 14281 : if (myProjection == nullptr) {
97 0 : throw ProcessError(TL("Could not build projection!"));
98 : }
99 : #endif
100 : }
101 206020 : }
102 :
103 :
104 : #ifdef PROJ_API_FILE
105 : void
106 42539 : GeoConvHelper::initProj(const std::string& proj) {
107 : #ifdef PROJ_VERSION_MAJOR
108 42539 : myProjection = proj_create(PJ_DEFAULT_CTX, proj.c_str());
109 42539 : checkError(myProjection);
110 : #if PROJ_VERSION_MAJOR > 5
111 42539 : if (myProjection != nullptr) {
112 42539 : const PJ_TYPE type = proj_get_type(myProjection);
113 42539 : if (type != PJ_TYPE_TRANSFORMATION
114 : && type != PJ_TYPE_CONCATENATED_OPERATION
115 42539 : && type != PJ_TYPE_OTHER_COORDINATE_OPERATION) {
116 : // handle PROJCS WKT (i.e. from Visum) which doesn't define a transformation but only CRS
117 0 : proj_destroy(myProjection);
118 0 : PJ_CONTEXT* ctx = proj_context_create();
119 0 : proj_context_use_proj4_init_rules(ctx, 1);
120 0 : myProjection = proj_create_crs_to_crs(ctx,
121 : "+proj=longlat +datum=WGS84 +type=crs +no_defs",
122 : proj.c_str(),
123 : NULL);
124 0 : checkError(myProjection);
125 : // "modern" proj doesn't default to radians but traditional proj strings do
126 0 : myRadians = false;
127 : //auto tmp = proj_normalize_for_visualization(ctx, myProjection);
128 : //PJ_PROJ_INFO info = proj_pj_info(myProjection);
129 : //printf("ID init: %s\n", info.id);
130 : //printf("Description: %s\n", info.description);
131 : //printf("Definition: %s\n", info.definition);
132 : //printf("Has Inverse: %s\n", info.has_inverse ? "Yes" : "No");
133 : }
134 : }
135 : #endif
136 : #else
137 : myProjection = pj_init_plus(proj.c_str());
138 : #endif
139 42539 : }
140 : #endif
141 :
142 :
143 206020 : GeoConvHelper::~GeoConvHelper() {
144 : #ifdef PROJ_API_FILE
145 206020 : if (myProjection != nullptr) {
146 : #ifdef PROJ_VERSION_MAJOR
147 43114 : proj_destroy(myProjection);
148 : #else
149 : pj_free(myProjection);
150 : #endif
151 : }
152 206020 : if (myInverseProjection != nullptr) {
153 : #ifdef PROJ_VERSION_MAJOR
154 0 : proj_destroy(myInverseProjection);
155 : #else
156 : pj_free(myInverseProjection);
157 : #endif
158 : }
159 206020 : if (myGeoProjection != nullptr) {
160 : #ifdef PROJ_VERSION_MAJOR
161 0 : proj_destroy(myGeoProjection);
162 : #else
163 : pj_free(myGeoProjection);
164 : #endif
165 : }
166 : #endif
167 206020 : }
168 :
169 : bool
170 2612 : GeoConvHelper::operator==(const GeoConvHelper& o) const {
171 : return (
172 2612 : myProjString == o.myProjString &&
173 370 : myOffset == o.myOffset &&
174 740 : myProjectionMethod == o.myProjectionMethod &&
175 370 : myOrigBoundary == o.myOrigBoundary &&
176 0 : myConvBoundary == o.myConvBoundary &&
177 0 : myGeoScale == o.myGeoScale &&
178 0 : myCos == o.myCos &&
179 0 : mySin == o.mySin &&
180 2612 : myUseInverseProjection == o.myUseInverseProjection &&
181 : myFlatten == o.myFlatten
182 2612 : );
183 : }
184 :
185 : GeoConvHelper&
186 100217 : GeoConvHelper::operator=(const GeoConvHelper& orig) {
187 100217 : myProjString = orig.myProjString;
188 100217 : myOffset = orig.myOffset;
189 100217 : myProjectionMethod = orig.myProjectionMethod;
190 : myOrigBoundary = orig.myOrigBoundary;
191 : myConvBoundary = orig.myConvBoundary;
192 100217 : myGeoScale = orig.myGeoScale;
193 100217 : myCos = orig.myCos;
194 100217 : mySin = orig.mySin;
195 100217 : myUseInverseProjection = orig.myUseInverseProjection;
196 100217 : myFlatten = orig.myFlatten;
197 : #ifdef PROJ_API_FILE
198 100217 : if (myProjection != nullptr) {
199 : #ifdef PROJ_VERSION_MAJOR
200 18 : proj_destroy(myProjection);
201 : #else
202 : pj_free(myProjection);
203 : #endif
204 18 : myProjection = nullptr;
205 : }
206 100217 : if (myInverseProjection != nullptr) {
207 : #ifdef PROJ_VERSION_MAJOR
208 0 : proj_destroy(myInverseProjection);
209 : #else
210 : pj_free(myInverseProjection);
211 : #endif
212 0 : myInverseProjection = nullptr;
213 : }
214 100217 : if (myGeoProjection != nullptr) {
215 : #ifdef PROJ_VERSION_MAJOR
216 0 : proj_destroy(myGeoProjection);
217 : #else
218 : pj_free(myGeoProjection);
219 : #endif
220 0 : myGeoProjection = nullptr;
221 : }
222 100217 : if (orig.myProjection != nullptr) {
223 28258 : initProj(myProjString);
224 : }
225 100217 : if (orig.myInverseProjection != nullptr) {
226 : #ifdef PROJ_VERSION_MAJOR
227 0 : myInverseProjection = orig.myInverseProjection;
228 : #else
229 : myInverseProjection = pj_init_plus(pj_get_def(orig.myInverseProjection, 0));
230 : #endif
231 : }
232 100217 : if (orig.myGeoProjection != nullptr) {
233 : #ifdef PROJ_VERSION_MAJOR
234 0 : myGeoProjection = orig.myGeoProjection;
235 : #else
236 : myGeoProjection = pj_init_plus(pj_get_def(orig.myGeoProjection, 0));
237 : #endif
238 : }
239 : #endif
240 100217 : return *this;
241 : }
242 :
243 :
244 : bool
245 2418 : GeoConvHelper::init(OptionsCont& oc) {
246 2418 : std::string proj = "!"; // the default
247 2418 : double scale = oc.getFloat("proj.scale");
248 2418 : double rot = oc.getFloat("proj.rotate");
249 4836 : Position offset = Position(oc.getFloat("offset.x"), oc.getFloat("offset.y"), oc.getFloat("offset.z"));
250 4836 : bool inverse = oc.exists("proj.inverse") && oc.getBool("proj.inverse");
251 4836 : bool flatten = oc.exists("flatten") && oc.getBool("flatten");
252 :
253 4836 : if (oc.getBool("simple-projection")) {
254 : proj = "-";
255 : }
256 :
257 : #ifdef PROJ_API_FILE
258 2418 : if (oc.getBool("proj.inverse") && oc.getString("proj") == "!") {
259 0 : WRITE_ERROR(TL("Inverse projection works only with explicit proj parameters."));
260 0 : return false;
261 : }
262 4836 : unsigned numProjections = oc.getBool("simple-projection") + oc.getBool("proj.utm") + oc.getBool("proj.dhdn") + oc.getBool("proj.dhdnutm") + (oc.getString("proj").length() > 1);
263 2418 : if (numProjections > 1) {
264 0 : WRITE_ERROR(TL("The projection method needs to be uniquely defined."));
265 0 : return false;
266 : }
267 :
268 4836 : if (oc.getBool("proj.utm")) {
269 : proj = "UTM";
270 3624 : } else if (oc.getBool("proj.dhdn")) {
271 : proj = "DHDN";
272 3624 : } else if (oc.getBool("proj.dhdnutm")) {
273 : proj = "DHDN_UTM";
274 3624 : } else if (!oc.isDefault("proj")) {
275 4 : proj = oc.getString("proj");
276 : }
277 : #endif
278 2418 : myProcessing = GeoConvHelper(proj, offset, Boundary(), Boundary(), scale, rot, inverse, flatten);
279 2418 : myFinal = myProcessing;
280 : return true;
281 : }
282 :
283 :
284 : void
285 45936 : GeoConvHelper::init(const std::string& proj, const Position& offset, const Boundary& orig,
286 : const Boundary& conv, double scale) {
287 45936 : myProcessing = GeoConvHelper(proj, offset, orig, conv, scale);
288 45936 : myProcessing.resolveAbstractProjection();
289 45936 : myFinal = myProcessing;
290 45936 : }
291 :
292 : void
293 47255 : GeoConvHelper::resolveAbstractProjection() {
294 : #ifdef PROJ_API_FILE
295 47255 : if (myProjection == nullptr &&
296 33395 : myProjectionMethod != NONE && myProjectionMethod != SIMPLE) {
297 : const std::string origProj = myProjString;
298 : // try to initialized projection based on origBoundary
299 1 : Position tmp = myOrigBoundary.getCenter();
300 1 : x2cartesian(tmp, false);
301 1 : if (myProjection == nullptr) {
302 0 : WRITE_WARNING("Failed to intialized projection '" + origProj + "' based on origBoundary centered on '" + toString(myOrigBoundary.getCenter()) + "'");
303 0 : myProjectionMethod = NONE;
304 : }
305 : }
306 : #endif
307 47255 : }
308 :
309 : void
310 2456 : GeoConvHelper::addProjectionOptions(OptionsCont& oc) {
311 2456 : oc.addOptionSubTopic("Projection");
312 :
313 2456 : oc.doRegister("simple-projection", new Option_Bool(false));
314 4912 : oc.addSynonyme("simple-projection", "proj.simple", true);
315 4912 : oc.addDescription("simple-projection", "Projection", TL("Uses a simple method for projection"));
316 :
317 2456 : oc.doRegister("proj.scale", new Option_Float(1.0));
318 4912 : oc.addDescription("proj.scale", "Projection", TL("Scaling factor for input coordinates"));
319 :
320 2456 : oc.doRegister("proj.rotate", new Option_Float(0.0));
321 4912 : oc.addDescription("proj.rotate", "Projection", TL("Rotation (clockwise degrees) for input coordinates"));
322 :
323 : #ifdef PROJ_API_FILE
324 2456 : oc.doRegister("proj.utm", new Option_Bool(false));
325 4912 : oc.addDescription("proj.utm", "Projection", TL("Determine the UTM zone (for a universal transversal mercator projection based on the WGS84 ellipsoid)"));
326 :
327 2456 : oc.doRegister("proj.dhdn", new Option_Bool(false));
328 4912 : oc.addDescription("proj.dhdn", "Projection", "Determine the DHDN zone (for a transversal mercator projection based on the bessel ellipsoid, \"Gauss-Krueger\")");
329 :
330 4912 : oc.doRegister("proj", new Option_String("!"));
331 4912 : oc.addDescription("proj", "Projection", TL("Uses STR as proj.4 definition for projection"));
332 :
333 2456 : oc.doRegister("proj.inverse", new Option_Bool(false));
334 4912 : oc.addDescription("proj.inverse", "Projection", TL("Inverses projection"));
335 :
336 2456 : oc.doRegister("proj.dhdnutm", new Option_Bool(false));
337 4912 : oc.addDescription("proj.dhdnutm", "Projection", TL("Convert from Gauss-Krueger to UTM"));
338 : #endif // PROJ_API_FILE
339 2456 : }
340 :
341 :
342 : bool
343 620507 : GeoConvHelper::usingGeoProjection() const {
344 620507 : return myProjectionMethod != NONE;
345 : }
346 :
347 :
348 : bool
349 344 : GeoConvHelper::usingInverseGeoProjection() const {
350 344 : return myUseInverseProjection;
351 : }
352 :
353 :
354 : void
355 23140 : GeoConvHelper::cartesian2geo(Position& cartesian) const {
356 23140 : cartesian.sub(getOffsetBase());
357 23140 : if (myProjectionMethod == NONE) {
358 2875 : return;
359 : }
360 20265 : if (myProjectionMethod == SIMPLE) {
361 0 : const double y = cartesian.y() / 111136.;
362 0 : const double x = cartesian.x() / 111320. / cos(DEG2RAD(y));
363 : cartesian.set(x, y);
364 0 : return;
365 : }
366 : #ifdef PROJ_API_FILE
367 : #ifdef PROJ_VERSION_MAJOR
368 20265 : PJ_COORD c = proj_coord(cartesian.x(), cartesian.y(), cartesian.z(), 0);
369 20265 : c = proj_trans(myProjection, PJ_INV, c);
370 20265 : checkError(myProjection);
371 20265 : if (myRadians) {
372 20265 : cartesian.set(proj_todeg(c.lp.lam), proj_todeg(c.lp.phi));
373 : } else {
374 : cartesian.set(c.lp.lam, c.lp.phi);
375 : }
376 : #else
377 : projUV p;
378 : p.u = cartesian.x();
379 : p.v = cartesian.y();
380 : p = pj_inv(p, myProjection);
381 : //!!! check pj_errno
382 : p.u *= RAD_TO_DEG;
383 : p.v *= RAD_TO_DEG;
384 : cartesian.set((double) p.u, (double) p.v);
385 : #endif
386 : #endif
387 : }
388 :
389 :
390 : #ifdef PROJ_API_FILE
391 : #ifdef PROJ_VERSION_MAJOR
392 : bool
393 317097 : GeoConvHelper::checkError(projPJ projection) const {
394 317097 : const int err_no = proj_context_errno(PJ_DEFAULT_CTX);
395 : const char* err_string = "";
396 317097 : if (err_no != 0) {
397 : #if PROJ_VERSION_MAJOR > 7
398 0 : err_string = proj_context_errno_string(PJ_DEFAULT_CTX, err_no);
399 : #else
400 : err_string = proj_errno_string(err_no);
401 : #endif
402 : }
403 317097 : if (projection == nullptr) {
404 0 : if (err_no == 0) {
405 0 : WRITE_WARNING(TL("Failed to create transformation, reason unknown."));
406 : } else {
407 0 : WRITE_WARNINGF(TL("Failed to create transformation, %."), err_string);
408 : }
409 0 : return false;
410 : }
411 317097 : if (err_no != 0) {
412 0 : WRITE_WARNINGF(TL("Failed to transform, %."), err_string);
413 0 : return false;
414 : }
415 : return true;
416 : }
417 : #endif
418 : #endif
419 :
420 :
421 : bool
422 760933 : GeoConvHelper::x2cartesian(Position& from, bool includeInBoundary) {
423 760933 : if (includeInBoundary) {
424 556387 : myOrigBoundary.add(from);
425 : }
426 : // init projection parameter on first use
427 : #ifdef PROJ_API_FILE
428 760933 : if (myProjection == nullptr) {
429 508481 : double x = from.x() * myGeoScale;
430 508481 : switch (myProjectionMethod) {
431 0 : case DHDN_UTM: {
432 0 : int zone = (int)((x - 500000.) / 1000000.);
433 0 : if (zone < 1 || zone > 5) {
434 0 : WRITE_WARNING("Attempt to initialize DHDN_UTM-projection on invalid longitude " + toString(x));
435 0 : return false;
436 : }
437 0 : myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
438 0 : " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
439 0 : " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
440 : #ifdef PROJ_VERSION_MAJOR
441 0 : myInverseProjection = proj_create(PJ_DEFAULT_CTX, myProjString.c_str());
442 0 : checkError(myInverseProjection);
443 0 : myGeoProjection = proj_create(PJ_DEFAULT_CTX, "+proj=latlong +datum=WGS84");
444 0 : checkError(myGeoProjection);
445 : #else
446 : myInverseProjection = pj_init_plus(myProjString.c_str());
447 : myGeoProjection = pj_init_plus("+proj=latlong +datum=WGS84");
448 : #endif
449 0 : x = ((x - 500000.) / 1000000.) * 3; // continues with UTM
450 : }
451 : FALLTHROUGH;
452 593 : case UTM: {
453 593 : int zone = (int)(x + 180) / 6 + 1;
454 1186 : myProjString = "+proj=utm +zone=" + toString(zone) +
455 593 : " +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
456 : #ifdef PROJ_VERSION_MAJOR
457 593 : myProjection = proj_create(PJ_DEFAULT_CTX, myProjString.c_str());
458 593 : checkError(myProjection);
459 : #else
460 : myProjection = pj_init_plus(myProjString.c_str());
461 : #endif
462 : }
463 593 : break;
464 0 : case DHDN: {
465 0 : int zone = (int)(x / 3);
466 0 : if (zone < 1 || zone > 5) {
467 0 : WRITE_WARNING("Attempt to initialize DHDN-projection on invalid longitude " + toString(x));
468 0 : return false;
469 : }
470 0 : myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
471 0 : " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
472 0 : " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
473 : #ifdef PROJ_VERSION_MAJOR
474 0 : myProjection = proj_create(PJ_DEFAULT_CTX, myProjString.c_str());
475 0 : checkError(myProjection);
476 : #else
477 : myProjection = pj_init_plus(myProjString.c_str());
478 : #endif
479 : }
480 : break;
481 : default:
482 : break;
483 : }
484 : }
485 760933 : if (myInverseProjection != nullptr) {
486 : #ifdef PROJ_VERSION_MAJOR
487 0 : PJ_COORD c = proj_coord(from.x(), from.y(), from.z(), 0);
488 0 : c = proj_trans(myInverseProjection, PJ_INV, c);
489 0 : checkError(myInverseProjection);
490 0 : from.set(proj_todeg(c.lp.lam), proj_todeg(c.lp.phi));
491 : #else
492 : double x = from.x();
493 : double y = from.y();
494 : if (pj_transform(myInverseProjection, myGeoProjection, 1, 1, &x, &y, nullptr)) {
495 : WRITE_WARNINGF(TL("Could not transform (%,%)"), toString(x), toString(y));
496 : }
497 : from.set(double(x * RAD_TO_DEG), double(y * RAD_TO_DEG));
498 : #endif
499 : }
500 : #endif
501 : // perform conversion
502 760933 : bool ok = x2cartesian_const(from);
503 760933 : if (ok) {
504 760932 : if (includeInBoundary) {
505 556386 : myConvBoundary.add(from);
506 : }
507 : }
508 : return ok;
509 : }
510 :
511 :
512 : bool
513 761589 : GeoConvHelper::x2cartesian_const(Position& from) const {
514 761589 : const double x2 = from.x() * myGeoScale;
515 761589 : const double y2 = from.y() * myGeoScale;
516 761589 : double x = x2 * myCos - y2 * mySin;
517 761589 : double y = x2 * mySin + y2 * myCos;
518 761589 : if (myProjectionMethod == NONE) {
519 : // do nothing
520 253715 : } else if (myUseInverseProjection) {
521 0 : cartesian2geo(from);
522 : } else {
523 253715 : if (x > 180.1 || x < -180.1) {
524 0 : WRITE_WARNING("Invalid longitude " + toString(x));
525 0 : return false;
526 : }
527 253715 : if (y > 90.1 || y < -90.1) {
528 2 : WRITE_WARNING("Invalid latitude " + toString(y));
529 1 : return false;
530 : }
531 : #ifdef PROJ_API_FILE
532 253714 : if (myProjection != nullptr) {
533 : #ifdef PROJ_VERSION_MAJOR
534 253700 : PJ_COORD c = proj_coord(myRadians ? proj_torad(x) : x, myRadians ? proj_torad(y) : y, from.z(), 0);
535 253700 : c = proj_trans(myProjection, PJ_FWD, c);
536 253700 : checkError(myProjection);
537 253700 : x = c.xy.x;
538 253700 : y = c.xy.y;
539 : #else
540 : projUV p;
541 : p.u = x * DEG_TO_RAD;
542 : p.v = y * DEG_TO_RAD;
543 : p = pj_fwd(p, myProjection);
544 : //!!! check pj_errno
545 : x = p.u;
546 : y = p.v;
547 : #endif
548 : }
549 : #endif
550 253714 : if (myProjectionMethod == SIMPLE) {
551 : // Sinusoidal projection (https://en.wikipedia.org/wiki/Sinusoidal_projection)
552 14 : x *= 111320. * cos(DEG2RAD(y));
553 14 : y *= 111136.;
554 : //!!! recheck whether the axes are mirrored
555 : }
556 : }
557 761588 : if (x > std::numeric_limits<double>::max() ||
558 761588 : y > std::numeric_limits<double>::max()) {
559 : return false;
560 : }
561 : from.set(x, y);
562 : from.add(myOffset);
563 761588 : if (myFlatten) {
564 : from.setz(0);
565 : }
566 : return true;
567 : }
568 :
569 :
570 : void
571 1650 : GeoConvHelper::moveConvertedBy(double x, double y) {
572 : myOffset.add(x, y);
573 1650 : myConvBoundary.moveby(x, y);
574 1650 : }
575 :
576 :
577 : const Boundary&
578 5552 : GeoConvHelper::getOrigBoundary() const {
579 5552 : return myOrigBoundary;
580 : }
581 :
582 :
583 : const Boundary&
584 20741 : GeoConvHelper::getConvBoundary() const {
585 20741 : return myConvBoundary;
586 : }
587 :
588 :
589 : const Position
590 4273 : GeoConvHelper::getOffset() const {
591 4273 : return myOffset;
592 : }
593 :
594 :
595 : const Position
596 27658 : GeoConvHelper::getOffsetBase() const {
597 27658 : return myOffset;
598 : }
599 :
600 :
601 : const std::string&
602 4516 : GeoConvHelper::getProjString() const {
603 4516 : return myProjString;
604 : }
605 :
606 :
607 : void
608 2227 : GeoConvHelper::computeFinal(bool lefthand) {
609 2227 : if (myNumLoaded == 0) {
610 1161 : myFinal = myProcessing;
611 1161 : if (lefthand) {
612 : myFinal.myOffset.mul(1, -1);
613 : }
614 : } else {
615 1066 : if (lefthand) {
616 : myProcessing.myOffset.mul(1, -1);
617 : }
618 4264 : myFinal = GeoConvHelper(
619 : // prefer options over loaded location
620 2132 : myProcessing.usingGeoProjection() ? myProcessing.getProjString() : myLoaded.getProjString(),
621 : // let offset and boundary lead back to the original coords of the loaded data
622 2132 : myProcessing.getOffset() + myLoaded.getOffset(),
623 : myLoaded.getOrigBoundary(),
624 : // the new boundary (updated during loading)
625 1066 : myProcessing.getConvBoundary());
626 : }
627 1075 : if (lefthand) {
628 23 : myFinal.myConvBoundary.flipY();
629 : }
630 2227 : }
631 :
632 :
633 : void
634 1295 : GeoConvHelper::setLoaded(const GeoConvHelper& loaded) {
635 1295 : myNumLoaded++;
636 1295 : if (myNumLoaded > 1) {
637 26 : WRITE_WARNINGF(TL("Ignoring loaded location attribute nr. % for tracking of original location"), toString(myNumLoaded));
638 : } else {
639 1282 : myLoaded = loaded;
640 : }
641 1295 : }
642 :
643 :
644 : void
645 569 : GeoConvHelper::setLoadedPlain(const std::string& nodFile, const GeoConvHelper& loaded) {
646 1138 : myLoadedPlain[nodFile] = {loaded.getProjString(), loaded.getOffset()};
647 569 : }
648 :
649 :
650 : GeoConvHelper*
651 1404 : GeoConvHelper::getLoadedPlain(const std::string& plainFile, const std::string& suffix) {
652 4212 : std::string nodFile = StringUtils::replace(plainFile, suffix, ".nod.xml");
653 : auto it = myLoadedPlain.find(nodFile);
654 1404 : if (it != myLoadedPlain.end()) {
655 698 : return new GeoConvHelper(it->second.first, it->second.second, Boundary(), Boundary());
656 : } else {
657 : return nullptr;
658 : }
659 : }
660 :
661 :
662 : void
663 0 : GeoConvHelper::resetLoaded() {
664 0 : myNumLoaded = 0;
665 : myLoadedPlain.clear();
666 0 : }
667 :
668 :
669 : void
670 2302 : GeoConvHelper::writeLocation(OutputDevice& into) {
671 2302 : into.openTag(SUMO_TAG_LOCATION);
672 2302 : into.writeAttr(SUMO_ATTR_NET_OFFSET, myFinal.getOffsetBase());
673 2302 : into.writeAttr(SUMO_ATTR_CONV_BOUNDARY, myFinal.getConvBoundary());
674 2302 : if (myFinal.usingGeoProjection()) {
675 867 : into.setPrecision(gPrecisionGeo);
676 : }
677 2302 : into.writeAttr(SUMO_ATTR_ORIG_BOUNDARY, myFinal.getOrigBoundary());
678 2302 : if (myFinal.usingGeoProjection()) {
679 867 : into.setPrecision();
680 : }
681 2302 : into.writeAttr(SUMO_ATTR_ORIG_PROJ, StringUtils::escapeXML(myFinal.getProjString()));
682 2302 : into.closeTag();
683 2302 : into.lf();
684 2302 : }
685 :
686 :
687 : /****************************************************************************/
|
