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 210360 : GeoConvHelper::GeoConvHelper(const std::string& proj, const Position& offset,
52 210360 : const Boundary& orig, const Boundary& conv, double scale, double rot, bool inverse, bool flatten):
53 210360 : myProjString(proj),
54 : #ifdef PROJ_API_FILE
55 210360 : myProjection(nullptr),
56 210360 : myInverseProjection(nullptr),
57 210360 : myGeoProjection(nullptr),
58 210360 : myRadians(true),
59 : #endif
60 210360 : myOffset(offset),
61 210360 : myGeoScale(scale),
62 210360 : mySin(sin(DEG2RAD(-rot))), // rotate clockwise
63 210360 : myCos(cos(DEG2RAD(-rot))),
64 210360 : myProjectionMethod(NONE),
65 210360 : myUseInverseProjection(inverse),
66 210360 : myFlatten(flatten),
67 : myOrigBoundary(orig),
68 : myConvBoundary(conv) {
69 : // older PROJ libraries fail to construct inverse projection if this
70 : // string is present
71 631080 : myProjString = StringUtils::replace(myProjString, "+type=crs", "");
72 :
73 210360 : if (proj == "!") {
74 195264 : myProjectionMethod = NONE;
75 15096 : } else if (proj == "-") {
76 3 : myProjectionMethod = SIMPLE;
77 15093 : } else if (proj == "UTM") {
78 609 : myProjectionMethod = UTM;
79 14484 : } else if (proj == "DHDN") {
80 0 : myProjectionMethod = DHDN;
81 14484 : } else if (proj == "DHDN_UTM") {
82 0 : myProjectionMethod = DHDN_UTM;
83 : #ifdef PROJ_API_FILE
84 : } else {
85 14484 : myProjectionMethod = PROJ;
86 14484 : initProj(myProjString);
87 14484 : 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 14484 : if (myProjection == nullptr) {
97 0 : throw ProcessError(TL("Could not build projection!"));
98 : }
99 : #endif
100 : }
101 210360 : }
102 :
103 :
104 : #ifdef PROJ_API_FILE
105 : void
106 43100 : GeoConvHelper::initProj(const std::string& proj) {
107 : #ifdef PROJ_VERSION_MAJOR
108 43100 : myProjection = proj_create(PJ_DEFAULT_CTX, proj.c_str());
109 43100 : checkError(myProjection);
110 : #if PROJ_VERSION_MAJOR > 5
111 43100 : if (myProjection != nullptr) {
112 43100 : const PJ_TYPE type = proj_get_type(myProjection);
113 43100 : if (type != PJ_TYPE_TRANSFORMATION
114 : && type != PJ_TYPE_CONCATENATED_OPERATION
115 43100 : && 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 43100 : }
140 : #endif
141 :
142 :
143 210360 : GeoConvHelper::~GeoConvHelper() {
144 : #ifdef PROJ_API_FILE
145 210360 : if (myProjection != nullptr) {
146 : #ifdef PROJ_VERSION_MAJOR
147 43677 : proj_destroy(myProjection);
148 : #else
149 : pj_free(myProjection);
150 : #endif
151 : }
152 210360 : if (myInverseProjection != nullptr) {
153 : #ifdef PROJ_VERSION_MAJOR
154 0 : proj_destroy(myInverseProjection);
155 : #else
156 : pj_free(myInverseProjection);
157 : #endif
158 : }
159 210360 : 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 210360 : }
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 102351 : GeoConvHelper::operator=(const GeoConvHelper& orig) {
187 102351 : myProjString = orig.myProjString;
188 102351 : myOffset = orig.myOffset;
189 102351 : myProjectionMethod = orig.myProjectionMethod;
190 : myOrigBoundary = orig.myOrigBoundary;
191 : myConvBoundary = orig.myConvBoundary;
192 102351 : myGeoScale = orig.myGeoScale;
193 102351 : myCos = orig.myCos;
194 102351 : mySin = orig.mySin;
195 102351 : myUseInverseProjection = orig.myUseInverseProjection;
196 102351 : myFlatten = orig.myFlatten;
197 : #ifdef PROJ_API_FILE
198 102351 : 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 102351 : 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 102351 : 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 102351 : if (orig.myProjection != nullptr) {
223 28616 : initProj(myProjString);
224 : }
225 102351 : 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 102351 : 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 102351 : return *this;
241 : }
242 :
243 :
244 : bool
245 2456 : GeoConvHelper::init(OptionsCont& oc) {
246 2456 : std::string proj = "!"; // the default
247 2456 : double scale = oc.getFloat("proj.scale");
248 2456 : double rot = oc.getFloat("proj.rotate");
249 4912 : Position offset = Position(oc.getFloat("offset.x"), oc.getFloat("offset.y"), oc.getFloat("offset.z"));
250 4912 : bool inverse = oc.exists("proj.inverse") && oc.getBool("proj.inverse");
251 4912 : bool flatten = oc.exists("flatten") && oc.getBool("flatten");
252 :
253 4912 : if (oc.getBool("simple-projection")) {
254 : proj = "-";
255 : }
256 :
257 : #ifdef PROJ_API_FILE
258 : #ifdef PROJ_VERSION_MAJOR
259 : #ifdef WIN32
260 : if (getenv("SUMO_HOME") != nullptr) {
261 : const char* paths[] = { (std::string(getenv("SUMO_HOME")) + "/data/proj").c_str() };
262 : proj_context_set_search_paths(PJ_DEFAULT_CTX, 1, paths);
263 : }
264 : #endif
265 : #endif
266 2456 : if (oc.getBool("proj.inverse") && oc.getString("proj") == "!") {
267 0 : WRITE_ERROR(TL("Inverse projection works only with explicit proj parameters."));
268 0 : return false;
269 : }
270 4912 : unsigned numProjections = oc.getBool("simple-projection") + oc.getBool("proj.utm") + oc.getBool("proj.dhdn") + oc.getBool("proj.dhdnutm") + (oc.getString("proj").length() > 1);
271 2456 : if (numProjections > 1) {
272 0 : WRITE_ERROR(TL("The projection method needs to be uniquely defined."));
273 0 : return false;
274 : }
275 :
276 4912 : if (oc.getBool("proj.utm")) {
277 : proj = "UTM";
278 3696 : } else if (oc.getBool("proj.dhdn")) {
279 : proj = "DHDN";
280 3696 : } else if (oc.getBool("proj.dhdnutm")) {
281 : proj = "DHDN_UTM";
282 3696 : } else if (!oc.isDefault("proj")) {
283 4 : proj = oc.getString("proj");
284 : }
285 : #endif
286 2456 : myProcessing = GeoConvHelper(proj, offset, Boundary(), Boundary(), scale, rot, inverse, flatten);
287 2456 : myFinal = myProcessing;
288 : return true;
289 : }
290 :
291 :
292 : void
293 46928 : GeoConvHelper::init(const std::string& proj, const Position& offset, const Boundary& orig,
294 : const Boundary& conv, double scale) {
295 46928 : myProcessing = GeoConvHelper(proj, offset, orig, conv, scale);
296 46928 : myProcessing.resolveAbstractProjection();
297 46928 : myFinal = myProcessing;
298 46928 : }
299 :
300 : void
301 48283 : GeoConvHelper::resolveAbstractProjection() {
302 : #ifdef PROJ_API_FILE
303 48283 : if (myProjection == nullptr &&
304 34245 : myProjectionMethod != NONE && myProjectionMethod != SIMPLE) {
305 : const std::string origProj = myProjString;
306 : // try to initialized projection based on origBoundary
307 1 : Position tmp = myOrigBoundary.getCenter();
308 1 : x2cartesian(tmp, false);
309 1 : if (myProjection == nullptr) {
310 0 : WRITE_WARNING("Failed to intialized projection '" + origProj + "' based on origBoundary centered on '" + toString(myOrigBoundary.getCenter()) + "'");
311 0 : myProjectionMethod = NONE;
312 : }
313 : }
314 : #endif
315 48283 : }
316 :
317 : void
318 2494 : GeoConvHelper::addProjectionOptions(OptionsCont& oc) {
319 2494 : oc.addOptionSubTopic("Projection");
320 :
321 2494 : oc.doRegister("simple-projection", new Option_Bool(false));
322 4988 : oc.addSynonyme("simple-projection", "proj.simple", true);
323 4988 : oc.addDescription("simple-projection", "Projection", TL("Uses a simple method for projection"));
324 :
325 2494 : oc.doRegister("proj.scale", new Option_Float(1.0));
326 4988 : oc.addDescription("proj.scale", "Projection", TL("Scaling factor for input coordinates"));
327 :
328 2494 : oc.doRegister("proj.rotate", new Option_Float(0.0));
329 4988 : oc.addDescription("proj.rotate", "Projection", TL("Rotation (clockwise degrees) for input coordinates"));
330 :
331 : #ifdef PROJ_API_FILE
332 2494 : oc.doRegister("proj.utm", new Option_Bool(false));
333 4988 : oc.addDescription("proj.utm", "Projection", TL("Determine the UTM zone (for a universal transversal mercator projection based on the WGS84 ellipsoid)"));
334 :
335 2494 : oc.doRegister("proj.dhdn", new Option_Bool(false));
336 4988 : oc.addDescription("proj.dhdn", "Projection", "Determine the DHDN zone (for a transversal mercator projection based on the bessel ellipsoid, \"Gauss-Krueger\")");
337 :
338 4988 : oc.doRegister("proj", new Option_String("!"));
339 4988 : oc.addDescription("proj", "Projection", TL("Uses STR as proj.4 definition for projection"));
340 :
341 2494 : oc.doRegister("proj.inverse", new Option_Bool(false));
342 4988 : oc.addDescription("proj.inverse", "Projection", TL("Inverses projection"));
343 :
344 2494 : oc.doRegister("proj.dhdnutm", new Option_Bool(false));
345 4988 : oc.addDescription("proj.dhdnutm", "Projection", TL("Convert from Gauss-Krueger to UTM"));
346 : #endif // PROJ_API_FILE
347 2494 : }
348 :
349 :
350 : bool
351 972338 : GeoConvHelper::usingGeoProjection() const {
352 972338 : return myProjectionMethod != NONE;
353 : }
354 :
355 :
356 : bool
357 344 : GeoConvHelper::usingInverseGeoProjection() const {
358 344 : return myUseInverseProjection;
359 : }
360 :
361 :
362 : void
363 23173 : GeoConvHelper::cartesian2geo(Position& cartesian) const {
364 23173 : cartesian.sub(getOffsetBase());
365 23173 : if (myProjectionMethod == NONE) {
366 2908 : return;
367 : }
368 20265 : if (myProjectionMethod == SIMPLE) {
369 0 : const double y = cartesian.y() / 111136.;
370 0 : const double x = cartesian.x() / 111320. / cos(DEG2RAD(y));
371 : cartesian.set(x, y);
372 0 : return;
373 : }
374 : #ifdef PROJ_API_FILE
375 : #ifdef PROJ_VERSION_MAJOR
376 20265 : PJ_COORD c = proj_coord(cartesian.x(), cartesian.y(), cartesian.z(), 0);
377 20265 : c = proj_trans(myProjection, PJ_INV, c);
378 20265 : checkError(myProjection);
379 20265 : if (myRadians) {
380 20265 : cartesian.set(proj_todeg(c.lp.lam), proj_todeg(c.lp.phi));
381 : } else {
382 : cartesian.set(c.lp.lam, c.lp.phi);
383 : }
384 : #else
385 : projUV p;
386 : p.u = cartesian.x();
387 : p.v = cartesian.y();
388 : p = pj_inv(p, myProjection);
389 : //!!! check pj_errno
390 : p.u *= RAD_TO_DEG;
391 : p.v *= RAD_TO_DEG;
392 : cartesian.set((double) p.u, (double) p.v);
393 : #endif
394 : #endif
395 : }
396 :
397 :
398 : #ifdef PROJ_API_FILE
399 : #ifdef PROJ_VERSION_MAJOR
400 : bool
401 317834 : GeoConvHelper::checkError(projPJ projection) const {
402 317834 : const int err_no = proj_context_errno(PJ_DEFAULT_CTX);
403 : const char* err_string = "";
404 317834 : if (err_no != 0) {
405 : #if PROJ_VERSION_MAJOR > 7
406 0 : err_string = proj_context_errno_string(PJ_DEFAULT_CTX, err_no);
407 : #else
408 : err_string = proj_errno_string(err_no);
409 : #endif
410 : }
411 317834 : if (projection == nullptr) {
412 0 : if (err_no == 0) {
413 0 : WRITE_WARNING(TL("Failed to create transformation, reason unknown."));
414 : } else {
415 0 : WRITE_WARNINGF(TL("Failed to create transformation, %."), err_string);
416 : }
417 0 : return false;
418 : }
419 317834 : if (err_no != 0) {
420 0 : WRITE_WARNINGF(TL("Failed to transform, %."), err_string);
421 0 : return false;
422 : }
423 : return true;
424 : }
425 : #endif
426 : #endif
427 :
428 :
429 : bool
430 1112582 : GeoConvHelper::x2cartesian(Position& from, bool includeInBoundary) {
431 1112582 : if (includeInBoundary) {
432 734806 : myOrigBoundary.add(from);
433 : }
434 : // init projection parameter on first use
435 : #ifdef PROJ_API_FILE
436 1112582 : if (myProjection == nullptr) {
437 859958 : double x = from.x() * myGeoScale;
438 859958 : switch (myProjectionMethod) {
439 0 : case DHDN_UTM: {
440 0 : int zone = (int)((x - 500000.) / 1000000.);
441 0 : if (zone < 1 || zone > 5) {
442 0 : WRITE_WARNING("Attempt to initialize DHDN_UTM-projection on invalid longitude " + toString(x));
443 0 : return false;
444 : }
445 0 : myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
446 0 : " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
447 0 : " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
448 : #ifdef PROJ_VERSION_MAJOR
449 0 : myInverseProjection = proj_create(PJ_DEFAULT_CTX, myProjString.c_str());
450 0 : checkError(myInverseProjection);
451 0 : myGeoProjection = proj_create(PJ_DEFAULT_CTX, "+proj=latlong +datum=WGS84");
452 0 : checkError(myGeoProjection);
453 : #else
454 : myInverseProjection = pj_init_plus(myProjString.c_str());
455 : myGeoProjection = pj_init_plus("+proj=latlong +datum=WGS84");
456 : #endif
457 0 : x = ((x - 500000.) / 1000000.) * 3; // continues with UTM
458 : }
459 : FALLTHROUGH;
460 595 : case UTM: {
461 595 : int zone = (int)(x + 180) / 6 + 1;
462 1190 : myProjString = "+proj=utm +zone=" + toString(zone) +
463 595 : " +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
464 : #ifdef PROJ_VERSION_MAJOR
465 595 : myProjection = proj_create(PJ_DEFAULT_CTX, myProjString.c_str());
466 595 : checkError(myProjection);
467 : #else
468 : myProjection = pj_init_plus(myProjString.c_str());
469 : #endif
470 : }
471 595 : break;
472 0 : case DHDN: {
473 0 : int zone = (int)(x / 3);
474 0 : if (zone < 1 || zone > 5) {
475 0 : WRITE_WARNING("Attempt to initialize DHDN-projection on invalid longitude " + toString(x));
476 0 : return false;
477 : }
478 0 : myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
479 0 : " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
480 0 : " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
481 : #ifdef PROJ_VERSION_MAJOR
482 0 : myProjection = proj_create(PJ_DEFAULT_CTX, myProjString.c_str());
483 0 : checkError(myProjection);
484 : #else
485 : myProjection = pj_init_plus(myProjString.c_str());
486 : #endif
487 : }
488 : break;
489 : default:
490 : break;
491 : }
492 : }
493 1112582 : if (myInverseProjection != nullptr) {
494 : #ifdef PROJ_VERSION_MAJOR
495 0 : PJ_COORD c = proj_coord(from.x(), from.y(), from.z(), 0);
496 0 : c = proj_trans(myInverseProjection, PJ_INV, c);
497 0 : checkError(myInverseProjection);
498 0 : from.set(proj_todeg(c.lp.lam), proj_todeg(c.lp.phi));
499 : #else
500 : double x = from.x();
501 : double y = from.y();
502 : if (pj_transform(myInverseProjection, myGeoProjection, 1, 1, &x, &y, nullptr)) {
503 : WRITE_WARNINGF(TL("Could not transform (%,%)"), toString(x), toString(y));
504 : }
505 : from.set(double(x * RAD_TO_DEG), double(y * RAD_TO_DEG));
506 : #endif
507 : }
508 : #endif
509 : // perform conversion
510 1112582 : bool ok = x2cartesian_const(from);
511 1112582 : if (ok) {
512 1112581 : if (includeInBoundary) {
513 734805 : myConvBoundary.add(from);
514 : }
515 : }
516 : return ok;
517 : }
518 :
519 :
520 : bool
521 1113238 : GeoConvHelper::x2cartesian_const(Position& from) const {
522 1113238 : const double x2 = from.x() * myGeoScale;
523 1113238 : const double y2 = from.y() * myGeoScale;
524 1113238 : double x = x2 * myCos - y2 * mySin;
525 1113238 : double y = x2 * mySin + y2 * myCos;
526 1113238 : if (myProjectionMethod == NONE) {
527 : // do nothing
528 253889 : } else if (myUseInverseProjection) {
529 0 : cartesian2geo(from);
530 : } else {
531 253889 : if (x > 180.1 || x < -180.1) {
532 0 : WRITE_WARNING("Invalid longitude " + toString(x));
533 0 : return false;
534 : }
535 253889 : if (y > 90.1 || y < -90.1) {
536 2 : WRITE_WARNING("Invalid latitude " + toString(y));
537 1 : return false;
538 : }
539 : #ifdef PROJ_API_FILE
540 253888 : if (myProjection != nullptr) {
541 : #ifdef PROJ_VERSION_MAJOR
542 253874 : PJ_COORD c = proj_coord(myRadians ? proj_torad(x) : x, myRadians ? proj_torad(y) : y, from.z(), 0);
543 253874 : c = proj_trans(myProjection, PJ_FWD, c);
544 253874 : checkError(myProjection);
545 253874 : x = c.xy.x;
546 253874 : y = c.xy.y;
547 : #else
548 : projUV p;
549 : p.u = x * DEG_TO_RAD;
550 : p.v = y * DEG_TO_RAD;
551 : p = pj_fwd(p, myProjection);
552 : //!!! check pj_errno
553 : x = p.u;
554 : y = p.v;
555 : #endif
556 : }
557 : #endif
558 253888 : if (myProjectionMethod == SIMPLE) {
559 : // Sinusoidal projection (https://en.wikipedia.org/wiki/Sinusoidal_projection)
560 14 : x *= 111320. * cos(DEG2RAD(y));
561 14 : y *= 111136.;
562 : //!!! recheck whether the axes are mirrored
563 : }
564 : }
565 1113237 : if (x > std::numeric_limits<double>::max() ||
566 1113237 : y > std::numeric_limits<double>::max()) {
567 : return false;
568 : }
569 : from.set(x, y);
570 : from.add(myOffset);
571 1113237 : if (myFlatten) {
572 : from.setz(0);
573 : }
574 : return true;
575 : }
576 :
577 :
578 : void
579 1653 : GeoConvHelper::moveConvertedBy(double x, double y) {
580 : myOffset.add(x, y);
581 1653 : myConvBoundary.moveby(x, y);
582 1653 : }
583 :
584 :
585 : const Boundary&
586 5664 : GeoConvHelper::getOrigBoundary() const {
587 5664 : return myOrigBoundary;
588 : }
589 :
590 :
591 : const Boundary&
592 21143 : GeoConvHelper::getConvBoundary() const {
593 21143 : return myConvBoundary;
594 : }
595 :
596 :
597 : const Position
598 4436 : GeoConvHelper::getOffset() const {
599 4436 : return myOffset;
600 : }
601 :
602 :
603 : const Position
604 27767 : GeoConvHelper::getOffsetBase() const {
605 27767 : return myOffset;
606 : }
607 :
608 :
609 : const std::string&
610 4593 : GeoConvHelper::getProjString() const {
611 4593 : return myProjString;
612 : }
613 :
614 :
615 : void
616 2265 : GeoConvHelper::computeFinal(bool lefthand) {
617 2265 : if (myNumLoaded == 0) {
618 1163 : myFinal = myProcessing;
619 1163 : if (lefthand) {
620 : myFinal.myOffset.mul(1, -1);
621 : }
622 : } else {
623 1102 : if (lefthand) {
624 : myProcessing.myOffset.mul(1, -1);
625 : }
626 4408 : myFinal = GeoConvHelper(
627 : // prefer options over loaded location
628 2204 : myProcessing.usingGeoProjection() ? myProcessing.getProjString() : myLoaded.getProjString(),
629 : // let offset and boundary lead back to the original coords of the loaded data
630 2204 : myProcessing.getOffset() + myLoaded.getOffset(),
631 : myLoaded.getOrigBoundary(),
632 : // the new boundary (updated during loading)
633 1102 : myProcessing.getConvBoundary());
634 : }
635 1111 : if (lefthand) {
636 23 : myFinal.myConvBoundary.flipY();
637 : }
638 2265 : }
639 :
640 :
641 : void
642 1331 : GeoConvHelper::setLoaded(const GeoConvHelper& loaded) {
643 1331 : myNumLoaded++;
644 1331 : if (myNumLoaded > 1) {
645 26 : WRITE_WARNINGF(TL("Ignoring loaded location attribute nr. % for tracking of original location"), toString(myNumLoaded));
646 : } else {
647 1318 : myLoaded = loaded;
648 : }
649 1331 : }
650 :
651 :
652 : void
653 570 : GeoConvHelper::setLoadedPlain(const std::string& nodFile, const GeoConvHelper& loaded) {
654 1140 : myLoadedPlain[nodFile] = {loaded.getProjString(), loaded.getOffset()};
655 570 : }
656 :
657 :
658 : GeoConvHelper*
659 1416 : GeoConvHelper::getLoadedPlain(const std::string& plainFile, const std::string& suffix) {
660 4248 : std::string nodFile = StringUtils::replace(plainFile, suffix, ".nod.xml");
661 : auto it = myLoadedPlain.find(nodFile);
662 1416 : if (it != myLoadedPlain.end()) {
663 699 : return new GeoConvHelper(it->second.first, it->second.second, Boundary(), Boundary());
664 : } else {
665 : return nullptr;
666 : }
667 : }
668 :
669 :
670 : void
671 0 : GeoConvHelper::resetLoaded() {
672 0 : myNumLoaded = 0;
673 : myLoadedPlain.clear();
674 0 : }
675 :
676 :
677 : void
678 2340 : GeoConvHelper::writeLocation(OutputDevice& into) {
679 2340 : into.openTag(SUMO_TAG_LOCATION);
680 2340 : into.writeAttr(SUMO_ATTR_NET_OFFSET, myFinal.getOffsetBase());
681 2340 : into.writeAttr(SUMO_ATTR_CONV_BOUNDARY, myFinal.getConvBoundary());
682 2340 : if (myFinal.usingGeoProjection()) {
683 894 : into.setPrecision(gPrecisionGeo);
684 : }
685 2340 : into.writeAttr(SUMO_ATTR_ORIG_BOUNDARY, myFinal.getOrigBoundary());
686 2340 : if (myFinal.usingGeoProjection()) {
687 894 : into.setPrecision();
688 : }
689 2340 : into.writeAttr(SUMO_ATTR_ORIG_PROJ, StringUtils::escapeXML(myFinal.getProjString()));
690 2340 : into.closeTag();
691 2340 : into.lf();
692 2340 : }
693 :
694 :
695 : /****************************************************************************/
|
