SUMO - Simulation of Urban MObility
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GeoConvHelper.cpp
Go to the documentation of this file.
1 /****************************************************************************/
9 // static methods for processing the coordinates conversion for the current net
10 /****************************************************************************/
11 // SUMO, Simulation of Urban MObility; see http://sumo-sim.org/
12 // Copyright (C) 2001-2014 DLR (http://www.dlr.de/) and contributors
13 /****************************************************************************/
14 //
15 // This file is part of SUMO.
16 // SUMO is free software: you can redistribute it and/or modify
17 // it under the terms of the GNU General Public License as published by
18 // the Free Software Foundation, either version 3 of the License, or
19 // (at your option) any later version.
20 //
21 /****************************************************************************/
22 
23 
24 // ===========================================================================
25 // included modules
26 // ===========================================================================
27 #ifdef _MSC_VER
28 #include <windows_config.h>
29 #else
30 #include <config.h>
31 #endif
32 
33 #include <map>
34 #include <cmath>
35 #include <cassert>
36 #include <climits>
38 #include <utils/common/ToString.h>
39 #include <utils/geom/GeomHelper.h>
41 #include "GeoConvHelper.h"
42 
43 #ifdef CHECK_MEMORY_LEAKS
44 #include <foreign/nvwa/debug_new.h>
45 #endif // CHECK_MEMORY_LEAKS
46 
47 
48 // ===========================================================================
49 // static member variables
50 // ===========================================================================
55 
56 // ===========================================================================
57 // method definitions
58 // ===========================================================================
59 GeoConvHelper::GeoConvHelper(const std::string& proj, const Position& offset,
60  const Boundary& orig, const Boundary& conv, int shift, bool inverse):
61  myProjString(proj),
62 #ifdef HAVE_PROJ
63  myProjection(0),
64  myInverseProjection(0),
65  myGeoProjection(0),
66 #endif
67  myOffset(offset),
68  myGeoScale(pow(10, (double) - shift)),
69  myProjectionMethod(NONE),
70  myUseInverseProjection(inverse),
71  myOrigBoundary(orig),
72  myConvBoundary(conv) {
73  if (proj == "!") {
75  } else if (proj == "-") {
77  } else if (proj == "UTM") {
79  } else if (proj == "DHDN") {
81  } else if (proj == "DHDN_UTM") {
83 #ifdef HAVE_PROJ
84  } else {
86  myProjection = pj_init_plus(proj.c_str());
87  if (myProjection == 0) {
88  // !!! check pj_errno
89  throw ProcessError("Could not build projection!");
90  }
91 #endif
92  }
93 }
94 
95 
97 #ifdef HAVE_PROJ
98  if (myProjection != 0) {
99  pj_free(myProjection);
100  }
101  if (myInverseProjection != 0) {
102  pj_free(myInverseProjection);
103  }
104  if (myGeoProjection != 0) {
105  pj_free(myInverseProjection);
106  }
107 #endif
108 }
109 
110 
113  myProjString = orig.myProjString;
114  myOffset = orig.myOffset;
118  myGeoScale = orig.myGeoScale;
120 #ifdef HAVE_PROJ
121  if (myProjection != 0) {
122  pj_free(myProjection);
123  myProjection = 0;
124  }
125  if (myInverseProjection != 0) {
126  pj_free(myInverseProjection);
127  myInverseProjection = 0;
128  }
129  if (myGeoProjection != 0) {
130  pj_free(myGeoProjection);
131  myGeoProjection = 0;
132  }
133  if (orig.myProjection != 0) {
134  myProjection = pj_init_plus(orig.myProjString.c_str());
135  }
136  if (orig.myInverseProjection != 0) {
137  myInverseProjection = pj_init_plus(pj_get_def(orig.myInverseProjection, 0));
138  }
139  if (orig.myGeoProjection != 0) {
140  myGeoProjection = pj_init_plus(pj_get_def(orig.myGeoProjection, 0));
141  }
142 #endif
143  return *this;
144 }
145 
146 
147 bool
149  std::string proj = "!"; // the default
150  int shift = oc.getInt("proj.scale");
151  Position offset = Position(oc.getFloat("offset.x"), oc.getFloat("offset.y"));
152  bool inverse = oc.exists("proj.inverse") && oc.getBool("proj.inverse");
153 
154  if (oc.getBool("simple-projection")) {
155  proj = "-";
156  }
157 
158 #ifdef HAVE_PROJ
159  if (oc.getBool("proj.inverse") && oc.getString("proj") == "!") {
160  WRITE_ERROR("Inverse projection works only with explicit proj parameters.");
161  return false;
162  }
163  unsigned numProjections = oc.getBool("simple-projection") + oc.getBool("proj.utm") + oc.getBool("proj.dhdn") + oc.getBool("proj.dhdnutm") + (oc.getString("proj").length() > 1);
164  if (numProjections > 1) {
165  WRITE_ERROR("The projection method needs to be uniquely defined.");
166  return false;
167  }
168 
169  if (oc.getBool("proj.utm")) {
170  proj = "UTM";
171  } else if (oc.getBool("proj.dhdn")) {
172  proj = "DHDN";
173  } else if (oc.getBool("proj.dhdnutm")) {
174  proj = "DHDN_UTM";
175  } else {
176  proj = oc.getString("proj");
177  }
178 #endif
179  myProcessing = GeoConvHelper(proj, offset, Boundary(), Boundary(), shift, inverse);
181  return true;
182 }
183 
184 
185 void
186 GeoConvHelper::init(const std::string& proj,
187  const Position& offset,
188  const Boundary& orig,
189  const Boundary& conv,
190  int shift) {
191  myProcessing = GeoConvHelper(proj, offset, orig, conv, shift);
193 }
194 
195 
196 void
198  oc.addOptionSubTopic("Projection");
199 
200  oc.doRegister("simple-projection", new Option_Bool(false));
201  oc.addSynonyme("simple-projection", "proj.simple", true);
202  oc.addDescription("simple-projection", "Projection", "Uses a simple method for projection");
203 
204  oc.doRegister("proj.scale", new Option_Integer(0));
205  oc.addSynonyme("proj.scale", "proj.shift", true);
206  oc.addDescription("proj.scale", "Projection", "Number of places to shift decimal point to right in geo-coordinates");
207 
208 #ifdef HAVE_PROJ
209  oc.doRegister("proj.utm", new Option_Bool(false));
210  oc.addDescription("proj.utm", "Projection", "Determine the UTM zone (for a universal transversal mercator projection based on the WGS84 ellipsoid)");
211 
212  oc.doRegister("proj.dhdn", new Option_Bool(false));
213  oc.addDescription("proj.dhdn", "Projection", "Determine the DHDN zone (for a transversal mercator projection based on the bessel ellipsoid, \"Gauss-Krueger\")");
214 
215  oc.doRegister("proj", new Option_String("!"));
216  oc.addDescription("proj", "Projection", "Uses STR as proj.4 definition for projection");
217 
218  oc.doRegister("proj.inverse", new Option_Bool(false));
219  oc.addDescription("proj.inverse", "Projection", "Inverses projection");
220 
221  oc.doRegister("proj.dhdnutm", new Option_Bool(false));
222  oc.addDescription("proj.dhdnutm", "Projection", "Convert from Gauss-Krueger to UTM");
223 #endif // HAVE_PROJ
224 }
225 
226 
227 bool
229  return myProjectionMethod != NONE;
230 }
231 
232 
233 bool
235  return myUseInverseProjection;
236 }
237 
238 
239 void
241  cartesian.sub(getOffsetBase());
242  if (myProjectionMethod == NONE) {
243  return;
244  }
245 #ifdef HAVE_PROJ
246  projUV p;
247  p.u = cartesian.x();
248  p.v = cartesian.y();
249  p = pj_inv(p, myProjection);
251  p.u *= RAD_TO_DEG;
252  p.v *= RAD_TO_DEG;
253  cartesian.set((SUMOReal) p.u, (SUMOReal) p.v);
254 #endif
255 }
256 
257 
258 bool
259 GeoConvHelper::x2cartesian(Position& from, bool includeInBoundary) {
260  if (includeInBoundary) {
261  myOrigBoundary.add(from);
262  }
263  // init projection parameter on first use
264 #ifdef HAVE_PROJ
265  if (myProjection == 0) {
266  double x = from.x() * myGeoScale;
267  switch (myProjectionMethod) {
268  case DHDN_UTM: {
269  int zone = (int)((x - 500000.) / 1000000.);
270  if (zone < 1 || zone > 5) {
271  WRITE_WARNING("Attempt to initialize DHDN_UTM-projection on invalid longitude " + toString(x));
272  return false;
273  }
274  myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
275  " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
276  " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
277  myInverseProjection = pj_init_plus(myProjString.c_str());
278  myGeoProjection = pj_init_plus("+proj=latlong +datum=WGS84");
280  x = ((x - 500000.) / 1000000.) * 3; // continues with UTM
281  }
282  case UTM: {
283  int zone = (int)(x + 180) / 6 + 1;
284  myProjString = "+proj=utm +zone=" + toString(zone) +
285  " +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
286  myProjection = pj_init_plus(myProjString.c_str());
288  }
289  break;
290  case DHDN: {
291  int zone = (int)(x / 3);
292  if (zone < 1 || zone > 5) {
293  WRITE_WARNING("Attempt to initialize DHDN-projection on invalid longitude " + toString(x));
294  return false;
295  }
296  myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
297  " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
298  " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
299  myProjection = pj_init_plus(myProjString.c_str());
301  }
302  break;
303  default:
304  break;
305  }
306  }
307  if (myInverseProjection != 0) {
308  double x = from.x();
309  double y = from.y();
310  if (pj_transform(myInverseProjection, myGeoProjection, 1, 1, &x, &y, NULL)) {
311  WRITE_WARNING("Could not transform (" + toString(x) + "," + toString(y) + ")");
312  }
313  from.set(SUMOReal(x * RAD_TO_DEG), SUMOReal(y * RAD_TO_DEG));
314  }
315 #endif
316  // perform conversion
317  bool ok = x2cartesian_const(from);
318  if (ok) {
319  if (includeInBoundary) {
320  myConvBoundary.add(from);
321  }
322  }
323  return ok;
324 }
325 
326 
327 bool
329  double x = from.x() * myGeoScale;
330  double y = from.y() * myGeoScale;
331  if (myProjectionMethod == NONE) {
332  from.add(myOffset);
333  } else if (myUseInverseProjection) {
334  cartesian2geo(from);
335  } else {
336  if (x > 180.1 || x < -180.1 || y > 90.1 || y < -90.1) {
337  return false;
338  }
339 #ifdef HAVE_PROJ
340  if (myProjection != 0) {
341  projUV p;
342  p.u = x * DEG_TO_RAD;
343  p.v = y * DEG_TO_RAD;
344  p = pj_fwd(p, myProjection);
346  x = p.u;
347  y = p.v;
348  }
349 #endif
350  if (myProjectionMethod == SIMPLE) {
351  double ys = y;
352  x *= 111320. * cos(DEG2RAD(ys));
353  y *= 111136.;
354  from.set((SUMOReal)x, (SUMOReal)y);
356  from.add(myOffset);
357  }
358  }
361  return false;
362  }
363  if (myProjectionMethod != SIMPLE) {
364  from.set((SUMOReal)x, (SUMOReal)y);
365  from.add(myOffset);
366  }
367  return true;
368 }
369 
370 
371 void
373  myOffset.add(x, y);
374  myConvBoundary.moveby(x, y);
375 }
376 
377 
378 const Boundary&
380  return myOrigBoundary;
381 }
382 
383 
384 const Boundary&
386  return myConvBoundary;
387 }
388 
389 
390 const Position
392  return myOffset;
393 }
394 
395 
396 const Position
398  return myOffset;
399 }
400 
401 
402 const std::string&
404  return myProjString;
405 }
406 
407 
408 void
410  if (myNumLoaded == 0) {
412  } else {
414  // prefer options over loaded location
416  // let offset and boundary lead back to the original coords of the loaded data
419  // the new boundary (updated during loading)
421  }
422 }
423 
424 
425 void
427  myNumLoaded++;
428  if (myNumLoaded > 1) {
429  WRITE_WARNING("Ignoring loaded location attribute nr. " + toString(myNumLoaded) + " for tracking of original location");
430  } else {
431  myLoaded = loaded;
432  }
433 }
434 
435 
436 void
438  myNumLoaded = 0;
439 }
440 
441 
442 /****************************************************************************/
443