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-2013 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 #endif
65  myOffset(offset),
66  myGeoScale(pow(10, (double) - shift)),
67  myProjectionMethod(NONE),
68  myUseInverseProjection(inverse),
69  myOrigBoundary(orig),
70  myConvBoundary(conv) {
71  if (proj == "!") {
73  } else if (proj == "-") {
75  } else if (proj == "UTM") {
77  } else if (proj == "DHDN") {
79 #ifdef HAVE_PROJ
80  } else {
82  myProjection = pj_init_plus(proj.c_str());
83  if (myProjection == 0) {
84  // !!! check pj_errno
85  throw ProcessError("Could not build projection!");
86  }
87 #endif
88  }
89 }
90 
91 
93 #ifdef HAVE_PROJ
94  if (myProjection != 0) {
95  pj_free(myProjection);
96  }
97 #endif
98 }
99 
100 
103  myProjString = orig.myProjString;
104  myOffset = orig.myOffset;
108  myGeoScale = orig.myGeoScale;
110 #ifdef HAVE_PROJ
111  if (myProjection != 0) {
112  pj_free(myProjection);
113  myProjection = 0;
114  }
115  if (orig.myProjection != 0) {
116  myProjection = pj_init_plus(orig.myProjString.c_str());
117  }
118 #endif
119  return *this;
120 }
121 
122 
123 bool
125  std::string proj = "!"; // the default
126  int shift = oc.getInt("proj.scale");
127  Position offset = Position(oc.getFloat("offset.x"), oc.getFloat("offset.y"));
128  bool inverse = oc.exists("proj.inverse") && oc.getBool("proj.inverse");
129 
130  if (oc.getBool("simple-projection")) {
131  proj = "-";
132  }
133 
134 #ifdef HAVE_PROJ
135  if (oc.getBool("proj.inverse") && oc.getString("proj") == "!") {
136  WRITE_ERROR("Inverse projection works only with explicit proj parameters.");
137  return false;
138  }
139  unsigned numProjections = oc.getBool("simple-projection") + oc.getBool("proj.utm") + oc.getBool("proj.dhdn") + (oc.getString("proj").length() > 1);
140  if (numProjections > 1) {
141  WRITE_ERROR("The projection method needs to be uniquely defined.");
142  return false;
143  }
144 
145  if (oc.getBool("proj.utm")) {
146  proj = "UTM";
147  } else if (oc.getBool("proj.dhdn")) {
148  proj = "DHDN";
149  } else {
150  proj = oc.getString("proj");
151  }
152 #endif
153  myProcessing = GeoConvHelper(proj, offset, Boundary(), Boundary(), shift, inverse);
155  return true;
156 }
157 
158 
159 void
160 GeoConvHelper::init(const std::string& proj,
161  const Position& offset,
162  const Boundary& orig,
163  const Boundary& conv,
164  int shift) {
165  myProcessing = GeoConvHelper(proj, offset, orig, conv, shift);
167 }
168 
169 
170 void
172  oc.addOptionSubTopic("Projection");
173 
174  oc.doRegister("simple-projection", new Option_Bool(false));
175  oc.addSynonyme("simple-projection", "proj.simple", true);
176  oc.addDescription("simple-projection", "Projection", "Uses a simple method for projection");
177 
178  oc.doRegister("proj.scale", new Option_Integer(0));
179  oc.addSynonyme("proj.scale", "proj.shift", true);
180  oc.addDescription("proj.scale", "Projection", "Number of places to shift decimal point to right in geo-coordinates");
181 
182 #ifdef HAVE_PROJ
183  oc.doRegister("proj.utm", new Option_Bool(false));
184  oc.addDescription("proj.utm", "Projection", "Determine the UTM zone (for a universal transversal mercator projection based on the WGS84 ellipsoid)");
185 
186  oc.doRegister("proj.dhdn", new Option_Bool(false));
187  oc.addDescription("proj.dhdn", "Projection", "Determine the DHDN zone (for a transversal mercator projection based on the bessel ellipsoid)");
188 
189  oc.doRegister("proj", new Option_String("!"));
190  oc.addDescription("proj", "Projection", "Uses STR as proj.4 definition for projection");
191 
192  oc.doRegister("proj.inverse", new Option_Bool(false));
193  oc.addDescription("proj.inverse", "Projection", "Inverses projection");
194 #endif // HAVE_PROJ
195 }
196 
197 
198 bool
200  return myProjectionMethod != NONE;
201 }
202 
203 
204 bool
206  return myUseInverseProjection;
207 }
208 
209 
210 void
212  cartesian.sub(getOffsetBase());
213  if (myProjectionMethod == NONE) {
214  return;
215  }
216 #ifdef HAVE_PROJ
217  projUV p;
218  p.u = cartesian.x();
219  p.v = cartesian.y();
220  p = pj_inv(p, myProjection);
222  p.u *= RAD_TO_DEG;
223  p.v *= RAD_TO_DEG;
224  cartesian.set((SUMOReal) p.u, (SUMOReal) p.v);
225 #endif
226 }
227 
228 
229 bool
230 GeoConvHelper::x2cartesian(Position& from, bool includeInBoundary) {
231  myOrigBoundary.add(from);
232  // init projection parameter on first use
233 #ifdef HAVE_PROJ
234  if (myProjection == 0) {
235  const double x = from.x() * myGeoScale;
236  switch (myProjectionMethod) {
237  case UTM: {
238  int zone = (int)(x + 180) / 6 + 1;
239  myProjString = "+proj=utm +zone=" + toString(zone) +
240  " +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
241  myProjection = pj_init_plus(myProjString.c_str());
243  }
244  break;
245  case DHDN: {
246  int zone = (int)(x / 3);
247  if (zone < 1 || zone > 5) {
248  WRITE_WARNING("Attempt to initialize DHDN-projection on invalid longitude " + toString(x));
249  return false;
250  }
251  myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
252  " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
253  " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
254  myProjection = pj_init_plus(myProjString.c_str());
256  }
257  break;
258  default:
259  break;
260  }
261  }
262 #endif
263  // perform conversion
264  bool ok = x2cartesian_const(from);
265  if (ok) {
266  if (includeInBoundary) {
267  myConvBoundary.add(from);
268  }
269  }
270  return ok;
271 }
272 
273 
274 bool
276  double x = from.x() * myGeoScale;
277  double y = from.y() * myGeoScale;
278  if (myProjectionMethod == NONE) {
279  from.add(myOffset);
280  } else if (myUseInverseProjection) {
281  cartesian2geo(from);
282  } else {
283  if (x > 180.1 || x < -180.1 || y > 90.1 || y < -90.1) {
284  return false;
285  }
286 #ifdef HAVE_PROJ
287  if (myProjection != 0) {
288  projUV p;
289  p.u = x * DEG_TO_RAD;
290  p.v = y * DEG_TO_RAD;
291  p = pj_fwd(p, myProjection);
293  x = p.u;
294  y = p.v;
295  }
296 #endif
297  if (myProjectionMethod == SIMPLE) {
298  double ys = y;
299  x *= 111320. * cos(ys * PI / 180.0);
300  y *= 111136.;
301  from.set((SUMOReal)x, (SUMOReal)y);
303  from.add(myOffset);
304  }
305  }
308  return false;
309  }
310  if (myProjectionMethod != SIMPLE) {
311  from.set((SUMOReal)x, (SUMOReal)y);
312  from.add(myOffset);
313  }
314  return true;
315 }
316 
317 
318 void
320  myOffset.add(x, y);
321  myConvBoundary.moveby(x, y);
322 }
323 
324 
325 const Boundary&
327  return myOrigBoundary;
328 }
329 
330 
331 const Boundary&
333  return myConvBoundary;
334 }
335 
336 
337 const Position
339  return myOffset;
340 }
341 
342 
343 const Position
345  return myOffset;
346 }
347 
348 
349 const std::string&
351  return myProjString;
352 }
353 
354 
355 void
357  if (myNumLoaded == 0) {
359  } else {
361  // prefer options over loaded location
363  // let offset and boundary lead back to the original coords of the loaded data
366  // the new boundary (updated during loading)
368  }
369 }
370 
371 
372 void
374  myNumLoaded++;
375  if (myNumLoaded > 1) {
376  WRITE_WARNING("Ignoring loaded location attribute nr. " + toString(myNumLoaded) + " for tracking of original location");
377  } else {
378  myLoaded = loaded;
379  }
380 }
381 
382 
383 void
385  myNumLoaded = 0;
386 }
387 
388 
389 /****************************************************************************/
390