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.sourceforge.net/
12 // Copyright (C) 2001-2012 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>
36 #include <utils/common/ToString.h>
37 #include <utils/geom/GeomHelper.h>
39 #include "GeoConvHelper.h"
40 
41 #ifdef CHECK_MEMORY_LEAKS
42 #include <foreign/nvwa/debug_new.h>
43 #endif // CHECK_MEMORY_LEAKS
44 
45 
46 // ===========================================================================
47 // static member variables
48 // ===========================================================================
53 
54 // ===========================================================================
55 // method definitions
56 // ===========================================================================
57 GeoConvHelper::GeoConvHelper(const std::string& proj, const Position& offset,
58  const Boundary& orig, const Boundary& conv, int shift, bool inverse, bool baseFound):
59  myProjString(proj),
60 #ifdef HAVE_PROJ
61  myProjection(0),
62 #endif
63  myOffset(offset),
64  myGeoScale(pow(10, (double) - shift)),
65  myProjectionMethod(NONE),
66  myUseInverseProjection(inverse),
67  myBaseFound(baseFound),
68  myBaseX(0),
69  myBaseY(0),
70  myOrigBoundary(orig),
71  myConvBoundary(conv) {
72  if (proj == "!") {
74  } else if (proj == "-") {
76  } else if (proj == "UTM") {
78  } else if (proj == "DHDN") {
80 #ifdef HAVE_PROJ
81  } else {
83  myProjection = pj_init_plus(proj.c_str());
84  if (myProjection == 0) {
85  // !!! check pj_errno
86  throw ProcessError("Could not build projection!");
87  }
88 #endif
89  }
90 }
91 
92 
94 #ifdef HAVE_PROJ
95  if (myProjection != 0) {
96  pj_free(myProjection);
97  }
98 #endif
99 }
100 
101 
104  myProjString = orig.myProjString;
105  myOffset = orig.myOffset;
109  myGeoScale = orig.myGeoScale;
111  myBaseFound = orig.myBaseFound;
112  myBaseX = orig.myBaseX;
113  myBaseY = orig.myBaseY;
114 #ifdef HAVE_PROJ
115  if (myProjection != 0) {
116  pj_free(myProjection);
117  myProjection = 0;
118  }
119  if (orig.myProjection != 0) {
120  myProjection = pj_init_plus(orig.myProjString.c_str());
121  }
122 #endif
123  return *this;
124 }
125 
126 
127 bool
129  std::string proj = "!"; // the default
130  int shift = oc.getInt("proj.scale");
131  Position offset = Position(oc.getFloat("offset.x"), oc.getFloat("offset.y"));
132  bool inverse = oc.exists("proj.inverse") && oc.getBool("proj.inverse");
133  bool baseFound = !oc.exists("offset.disable-normalization") ||
134  oc.getBool("offset.disable-normalization") ||
135  !oc.isDefault("offset.x") ||
136  !oc.isDefault("offset.y");
137 
138  if (oc.getBool("simple-projection")) {
139  proj = "-";
140  }
141 
142 #ifdef HAVE_PROJ
143  if (oc.getBool("proj.inverse") && oc.getString("proj") == "!") {
144  WRITE_ERROR("Inverse projection works only with explicit proj parameters.");
145  return false;
146  }
147  unsigned numProjections = oc.getBool("simple-projection") + oc.getBool("proj.utm") + oc.getBool("proj.dhdn") + (oc.getString("proj").length() > 1);
148  if (numProjections > 1) {
149  WRITE_ERROR("The projection method needs to be uniquely defined.");
150  return false;
151  }
152 
153  if (oc.getBool("proj.utm")) {
154  proj = "UTM";
155  } else if (oc.getBool("proj.dhdn")) {
156  proj = "DHDN";
157  } else {
158  proj = oc.getString("proj");
159  }
160 #endif
161  myProcessing = GeoConvHelper(proj, offset, Boundary(), Boundary(), shift, inverse, baseFound);
163  return true;
164 }
165 
166 
167 void
168 GeoConvHelper::init(const std::string& proj,
169  const Position& offset,
170  const Boundary& orig,
171  const Boundary& conv) {
172  myProcessing = GeoConvHelper(proj, offset, orig, conv);
174 }
175 
176 
177 void
179  oc.addOptionSubTopic("Projection");
180 
181  oc.doRegister("simple-projection", new Option_Bool(false));
182  oc.addSynonyme("simple-projection", "proj.simple", true);
183  oc.addDescription("simple-projection", "Projection", "Uses a simple method for projection");
184 
185  oc.doRegister("proj.scale", new Option_Integer(0));
186  oc.addSynonyme("proj.scale", "proj.shift", true);
187  oc.addDescription("proj.scale", "Projection", "Number of places to shift decimal point to right in geo-coordinates");
188 
189 #ifdef HAVE_PROJ
190  oc.doRegister("proj.utm", new Option_Bool(false));
191  oc.addDescription("proj.utm", "Projection", "Determine the UTM zone (for a universal transversal mercator projection based on the WGS84 ellipsoid)");
192 
193  oc.doRegister("proj.dhdn", new Option_Bool(false));
194  oc.addDescription("proj.dhdn", "Projection", "Determine the DHDN zone (for a transversal mercator projection based on the bessel ellipsoid)");
195 
196  oc.doRegister("proj", new Option_String("!"));
197  oc.addDescription("proj", "Projection", "Uses STR as proj.4 definition for projection");
198 
199  oc.doRegister("proj.inverse", new Option_Bool(false));
200  oc.addDescription("proj.inverse", "Projection", "Inverses projection");
201 #endif // HAVE_PROJ
202 }
203 
204 
205 bool
207  return myProjectionMethod != NONE;
208 }
209 
210 
211 bool
213  return myUseInverseProjection;
214 }
215 
216 
217 void
219  cartesian.sub(getOffsetBase());
220  if (myProjectionMethod == NONE) {
221  return;
222  }
223 #ifdef HAVE_PROJ
224  projUV p;
225  p.u = cartesian.x();
226  p.v = cartesian.y();
227  p = pj_inv(p, myProjection);
229  p.u *= RAD_TO_DEG;
230  p.v *= RAD_TO_DEG;
231  cartesian.set((SUMOReal) p.u, (SUMOReal) p.v);
232 #endif
233 }
234 
235 
236 bool
237 GeoConvHelper::x2cartesian(Position& from, bool includeInBoundary) {
238  myOrigBoundary.add(from);
239  // init projection parameter on first use
240 #ifdef HAVE_PROJ
241  if (myProjection == 0) {
242  const double x = from.x() * myGeoScale;
243  switch (myProjectionMethod) {
244  case UTM: {
245  int zone = (int)(x + 180) / 6 + 1;
246  myProjString = "+proj=utm +zone=" + toString(zone) +
247  " +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
248  myProjection = pj_init_plus(myProjString.c_str());
250  }
251  break;
252  case DHDN: {
253  int zone = (int)(x / 3);
254  if (zone < 1 || zone > 5) {
255  WRITE_WARNING("Attempt to initialize DHDN-projection on invalid longitude " + toString(x));
256  return false;
257  }
258  myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
259  " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
260  " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
261  myProjection = pj_init_plus(myProjString.c_str());
263  }
264  break;
265  default:
266  break;
267  }
268  }
269 #endif
270  // perform conversion
271  bool ok = x2cartesian_const(from);
272  if (ok) {
273  if (!myBaseFound) {
274  // avoid very large coordinates to reduce numerical errors
275  if (myProjectionMethod == SIMPLE || from.x() > 100000 || from.y() > 100000) {
276  myBaseX = from.x();
277  myBaseY = from.y();
278  from.set(0, 0);
279  }
280  myBaseFound = true;
281  }
282 
283  if (includeInBoundary) {
284  myConvBoundary.add(from);
285  }
286  }
287  return ok;
288 }
289 
290 
291 bool
293  double x = from.x() * myGeoScale;
294  double y = from.y() * myGeoScale;
295  if (myProjectionMethod == NONE) {
296  from.add(myOffset);
297  } else if (myUseInverseProjection) {
298  cartesian2geo(from);
299  } else {
300  if (x > 180.1 || x < -180.1 || y > 90.1 || y < -90.1) {
301  return false;
302  }
303 #ifdef HAVE_PROJ
304  if (myProjection != 0) {
305  projUV p;
306  p.u = x * DEG_TO_RAD;
307  p.v = y * DEG_TO_RAD;
308  p = pj_fwd(p, myProjection);
310  x = p.u;
311  y = p.v;
312  }
313 #endif
314  if (myProjectionMethod == SIMPLE) {
315  double ys = y;
316  x -= myBaseX;
317  y -= myBaseY;
318  x *= 111320. * cos(ys * PI / 180.0);
319  y *= 111136.;
320  from.set((SUMOReal)x, (SUMOReal)y);
322  from.add(myOffset);
323  }
324  }
325  if (myProjectionMethod != SIMPLE) {
326  x -= myBaseX;
327  y -= myBaseY;
328  from.set((SUMOReal)x, (SUMOReal)y);
329  from.add(myOffset);
330  }
331  return true;
332 }
333 
334 
335 void
337  myOffset.add(x, y);
338  myConvBoundary.moveby(x, y);
339 }
340 
341 
342 const Boundary&
344  return myOrigBoundary;
345 }
346 
347 
348 const Boundary&
350  return myConvBoundary;
351 }
352 
353 
354 const Position
356  return myOffset;
357 }
358 
359 
360 const Position
362  return Position(myOffset.x() - myBaseX, myOffset.y() - myBaseY);
363 }
364 
365 
366 const std::string&
368  return myProjString;
369 }
370 
371 
372 void
374  if (myNumLoaded == 0) {
376  } else {
378  // prefer options over loaded location
380  // let offset and boundary lead back to the original coords of the loaded data
383  // the new boundary (updated during loading)
385  }
386 }
387 
388 
389 void
391  myNumLoaded++;
392  if (myNumLoaded > 1) {
393  WRITE_WARNING("Ignoring loaded location attribute nr. " + toString(myNumLoaded) + " for tracking of original location");
394  } else {
395  myLoaded = loaded;
396  }
397 }
398 
399 
400 void
402  myNumLoaded = 0;
403 }
404 
405 
406 /****************************************************************************/
407