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