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