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-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  if (includeInBoundary) {
232  myOrigBoundary.add(from);
233  }
234  // init projection parameter on first use
235 #ifdef HAVE_PROJ
236  if (myProjection == 0) {
237  const double x = from.x() * myGeoScale;
238  switch (myProjectionMethod) {
239  case UTM: {
240  int zone = (int)(x + 180) / 6 + 1;
241  myProjString = "+proj=utm +zone=" + toString(zone) +
242  " +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
243  myProjection = pj_init_plus(myProjString.c_str());
245  }
246  break;
247  case DHDN: {
248  int zone = (int)(x / 3);
249  if (zone < 1 || zone > 5) {
250  WRITE_WARNING("Attempt to initialize DHDN-projection on invalid longitude " + toString(x));
251  return false;
252  }
253  myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
254  " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
255  " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
256  myProjection = pj_init_plus(myProjString.c_str());
258  }
259  break;
260  default:
261  break;
262  }
263  }
264 #endif
265  // perform conversion
266  bool ok = x2cartesian_const(from);
267  if (ok) {
268  if (includeInBoundary) {
269  myConvBoundary.add(from);
270  }
271  }
272  return ok;
273 }
274 
275 
276 bool
278  double x = from.x() * myGeoScale;
279  double y = from.y() * myGeoScale;
280  if (myProjectionMethod == NONE) {
281  from.add(myOffset);
282  } else if (myUseInverseProjection) {
283  cartesian2geo(from);
284  } else {
285  if (x > 180.1 || x < -180.1 || y > 90.1 || y < -90.1) {
286  return false;
287  }
288 #ifdef HAVE_PROJ
289  if (myProjection != 0) {
290  projUV p;
291  p.u = x * DEG_TO_RAD;
292  p.v = y * DEG_TO_RAD;
293  p = pj_fwd(p, myProjection);
295  x = p.u;
296  y = p.v;
297  }
298 #endif
299  if (myProjectionMethod == SIMPLE) {
300  double ys = y;
301  x *= 111320. * cos(DEG2RAD(ys));
302  y *= 111136.;
303  from.set((SUMOReal)x, (SUMOReal)y);
305  from.add(myOffset);
306  }
307  }
310  return false;
311  }
312  if (myProjectionMethod != SIMPLE) {
313  from.set((SUMOReal)x, (SUMOReal)y);
314  from.add(myOffset);
315  }
316  return true;
317 }
318 
319 
320 void
322  myOffset.add(x, y);
323  myConvBoundary.moveby(x, y);
324 }
325 
326 
327 const Boundary&
329  return myOrigBoundary;
330 }
331 
332 
333 const Boundary&
335  return myConvBoundary;
336 }
337 
338 
339 const Position
341  return myOffset;
342 }
343 
344 
345 const Position
347  return myOffset;
348 }
349 
350 
351 const std::string&
353  return myProjString;
354 }
355 
356 
357 void
359  if (myNumLoaded == 0) {
361  } else {
363  // prefer options over loaded location
365  // let offset and boundary lead back to the original coords of the loaded data
368  // the new boundary (updated during loading)
370  }
371 }
372 
373 
374 void
376  myNumLoaded++;
377  if (myNumLoaded > 1) {
378  WRITE_WARNING("Ignoring loaded location attribute nr. " + toString(myNumLoaded) + " for tracking of original location");
379  } else {
380  myLoaded = loaded;
381  }
382 }
383 
384 
385 void
387  myNumLoaded = 0;
388 }
389 
390 
391 /****************************************************************************/
392