SUMO - Simulation of Urban MObility
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.dlr.de/
12 // Copyright (C) 2001-2016 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>
42 #include "GeoConvHelper.h"
43 
44 #ifdef CHECK_MEMORY_LEAKS
45 #include <foreign/nvwa/debug_new.h>
46 #endif // CHECK_MEMORY_LEAKS
47 
48 
49 // ===========================================================================
50 // static member variables
51 // ===========================================================================
56 
57 // ===========================================================================
58 // method definitions
59 // ===========================================================================
60 GeoConvHelper::GeoConvHelper(const std::string& proj, const Position& offset,
61  const Boundary& orig, const Boundary& conv, int shift, bool inverse):
62  myProjString(proj),
63 #ifdef HAVE_PROJ
64  myProjection(0),
65  myInverseProjection(0),
66  myGeoProjection(0),
67 #endif
68  myOffset(offset),
69  myGeoScale(pow(10, (double) - shift)),
70  myProjectionMethod(NONE),
71  myUseInverseProjection(inverse),
72  myOrigBoundary(orig),
73  myConvBoundary(conv) {
74  if (proj == "!") {
76  } else if (proj == "-") {
78  } else if (proj == "UTM") {
80  } else if (proj == "DHDN") {
82  } else if (proj == "DHDN_UTM") {
84 #ifdef HAVE_PROJ
85  } else {
87  myProjection = pj_init_plus(proj.c_str());
88  if (myProjection == 0) {
89  // !!! check pj_errno
90  throw ProcessError("Could not build projection!");
91  }
92 #endif
93  }
94 }
95 
96 
98 #ifdef HAVE_PROJ
99  if (myProjection != 0) {
100  pj_free(myProjection);
101  }
102  if (myInverseProjection != 0) {
103  pj_free(myInverseProjection);
104  }
105  if (myGeoProjection != 0) {
106  pj_free(myInverseProjection);
107  }
108 #endif
109 }
110 
111 
114  myProjString = orig.myProjString;
115  myOffset = orig.myOffset;
119  myGeoScale = orig.myGeoScale;
121 #ifdef HAVE_PROJ
122  if (myProjection != 0) {
123  pj_free(myProjection);
124  myProjection = 0;
125  }
126  if (myInverseProjection != 0) {
127  pj_free(myInverseProjection);
128  myInverseProjection = 0;
129  }
130  if (myGeoProjection != 0) {
131  pj_free(myGeoProjection);
132  myGeoProjection = 0;
133  }
134  if (orig.myProjection != 0) {
135  myProjection = pj_init_plus(orig.myProjString.c_str());
136  }
137  if (orig.myInverseProjection != 0) {
138  myInverseProjection = pj_init_plus(pj_get_def(orig.myInverseProjection, 0));
139  }
140  if (orig.myGeoProjection != 0) {
141  myGeoProjection = pj_init_plus(pj_get_def(orig.myGeoProjection, 0));
142  }
143 #endif
144  return *this;
145 }
146 
147 
148 bool
150  std::string proj = "!"; // the default
151  int shift = oc.getInt("proj.scale");
152  Position offset = Position(oc.getFloat("offset.x"), oc.getFloat("offset.y"));
153  bool inverse = oc.exists("proj.inverse") && oc.getBool("proj.inverse");
154 
155  if (oc.getBool("simple-projection")) {
156  proj = "-";
157  }
158 
159 #ifdef HAVE_PROJ
160  if (oc.getBool("proj.inverse") && oc.getString("proj") == "!") {
161  WRITE_ERROR("Inverse projection works only with explicit proj parameters.");
162  return false;
163  }
164  unsigned numProjections = oc.getBool("simple-projection") + oc.getBool("proj.utm") + oc.getBool("proj.dhdn") + oc.getBool("proj.dhdnutm") + (oc.getString("proj").length() > 1);
165  if (numProjections > 1) {
166  WRITE_ERROR("The projection method needs to be uniquely defined.");
167  return false;
168  }
169 
170  if (oc.getBool("proj.utm")) {
171  proj = "UTM";
172  } else if (oc.getBool("proj.dhdn")) {
173  proj = "DHDN";
174  } else if (oc.getBool("proj.dhdnutm")) {
175  proj = "DHDN_UTM";
176  } else if (!oc.isDefault("proj")) {
177  proj = oc.getString("proj");
178  }
179 #endif
180  myProcessing = GeoConvHelper(proj, offset, Boundary(), Boundary(), shift, inverse);
182  return true;
183 }
184 
185 
186 void
187 GeoConvHelper::init(const std::string& proj,
188  const Position& offset,
189  const Boundary& orig,
190  const Boundary& conv,
191  int shift) {
192  myProcessing = GeoConvHelper(proj, offset, orig, conv, shift);
194 }
195 
196 
197 void
199  oc.addOptionSubTopic("Projection");
200 
201  oc.doRegister("simple-projection", new Option_Bool(false));
202  oc.addSynonyme("simple-projection", "proj.simple", true);
203  oc.addDescription("simple-projection", "Projection", "Uses a simple method for projection");
204 
205  oc.doRegister("proj.scale", new Option_Integer(0));
206  oc.addSynonyme("proj.scale", "proj.shift", true);
207  oc.addDescription("proj.scale", "Projection", "Number of places to shift decimal point to right in geo-coordinates");
208 
209 #ifdef HAVE_PROJ
210  oc.doRegister("proj.utm", new Option_Bool(false));
211  oc.addDescription("proj.utm", "Projection", "Determine the UTM zone (for a universal transversal mercator projection based on the WGS84 ellipsoid)");
212 
213  oc.doRegister("proj.dhdn", new Option_Bool(false));
214  oc.addDescription("proj.dhdn", "Projection", "Determine the DHDN zone (for a transversal mercator projection based on the bessel ellipsoid, \"Gauss-Krueger\")");
215 
216  oc.doRegister("proj", new Option_String("!"));
217  oc.addDescription("proj", "Projection", "Uses STR as proj.4 definition for projection");
218 
219  oc.doRegister("proj.inverse", new Option_Bool(false));
220  oc.addDescription("proj.inverse", "Projection", "Inverses projection");
221 
222  oc.doRegister("proj.dhdnutm", new Option_Bool(false));
223  oc.addDescription("proj.dhdnutm", "Projection", "Convert from Gauss-Krueger to UTM");
224 #endif // HAVE_PROJ
225 }
226 
227 
228 bool
230  return myProjectionMethod != NONE;
231 }
232 
233 
234 bool
236  return myUseInverseProjection;
237 }
238 
239 
240 void
242  cartesian.sub(getOffsetBase());
243  if (myProjectionMethod == NONE) {
244  return;
245  }
246 #ifdef HAVE_PROJ
247  projUV p;
248  p.u = cartesian.x();
249  p.v = cartesian.y();
250  p = pj_inv(p, myProjection);
252  p.u *= RAD_TO_DEG;
253  p.v *= RAD_TO_DEG;
254  cartesian.set((SUMOReal) p.u, (SUMOReal) p.v);
255 #endif
256 }
257 
258 
259 bool
260 GeoConvHelper::x2cartesian(Position& from, bool includeInBoundary) {
261  if (includeInBoundary) {
262  myOrigBoundary.add(from);
263  }
264  // init projection parameter on first use
265 #ifdef HAVE_PROJ
266  if (myProjection == 0) {
267  double x = from.x() * myGeoScale;
268  switch (myProjectionMethod) {
269  case DHDN_UTM: {
270  int zone = (int)((x - 500000.) / 1000000.);
271  if (zone < 1 || zone > 5) {
272  WRITE_WARNING("Attempt to initialize DHDN_UTM-projection on invalid longitude " + toString(x));
273  return false;
274  }
275  myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
276  " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
277  " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
278  myInverseProjection = pj_init_plus(myProjString.c_str());
279  myGeoProjection = pj_init_plus("+proj=latlong +datum=WGS84");
281  x = ((x - 500000.) / 1000000.) * 3; // continues with UTM
282  }
283  case UTM: {
284  int zone = (int)(x + 180) / 6 + 1;
285  myProjString = "+proj=utm +zone=" + toString(zone) +
286  " +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
287  myProjection = pj_init_plus(myProjString.c_str());
289  }
290  break;
291  case DHDN: {
292  int zone = (int)(x / 3);
293  if (zone < 1 || zone > 5) {
294  WRITE_WARNING("Attempt to initialize DHDN-projection on invalid longitude " + toString(x));
295  return false;
296  }
297  myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
298  " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
299  " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
300  myProjection = pj_init_plus(myProjString.c_str());
302  }
303  break;
304  default:
305  break;
306  }
307  }
308  if (myInverseProjection != 0) {
309  double x = from.x();
310  double y = from.y();
311  if (pj_transform(myInverseProjection, myGeoProjection, 1, 1, &x, &y, NULL)) {
312  WRITE_WARNING("Could not transform (" + toString(x) + "," + toString(y) + ")");
313  }
314  from.set(SUMOReal(x * RAD_TO_DEG), SUMOReal(y * RAD_TO_DEG));
315  }
316 #endif
317  // perform conversion
318  bool ok = x2cartesian_const(from);
319  if (ok) {
320  if (includeInBoundary) {
321  myConvBoundary.add(from);
322  }
323  }
324  return ok;
325 }
326 
327 
328 bool
330  double x = from.x() * myGeoScale;
331  double y = from.y() * myGeoScale;
332  if (myProjectionMethod == NONE) {
333  from.add(myOffset);
334  } else if (myUseInverseProjection) {
335  cartesian2geo(from);
336  } else {
337  if (x > 180.1 || x < -180.1) {
338  WRITE_WARNING("Invalid longitude " + toString(x));
339  return false;
340  }
341  if (y > 90.1 || y < -90.1) {
342  WRITE_WARNING("Invalid latitude " + toString(y));
343  return false;
344  }
345 #ifdef HAVE_PROJ
346  if (myProjection != 0) {
347  projUV p;
348  p.u = x * DEG_TO_RAD;
349  p.v = y * DEG_TO_RAD;
350  p = pj_fwd(p, myProjection);
352  x = p.u;
353  y = p.v;
354  }
355 #endif
356  if (myProjectionMethod == SIMPLE) {
357  double ys = y;
358  x *= 111320. * cos(DEG2RAD(ys));
359  y *= 111136.;
360  from.set((SUMOReal)x, (SUMOReal)y);
362  from.add(myOffset);
363  }
364  }
367  return false;
368  }
369  if (myProjectionMethod != SIMPLE) {
370  from.set((SUMOReal)x, (SUMOReal)y);
371  from.add(myOffset);
372  }
373  return true;
374 }
375 
376 
377 void
379  myOffset.add(x, y);
380  myConvBoundary.moveby(x, y);
381 }
382 
383 
384 const Boundary&
386  return myOrigBoundary;
387 }
388 
389 
390 const Boundary&
392  return myConvBoundary;
393 }
394 
395 
396 const Position
398  return myOffset;
399 }
400 
401 
402 const Position
404  return myOffset;
405 }
406 
407 
408 const std::string&
410  return myProjString;
411 }
412 
413 
414 void
416  if (myNumLoaded == 0) {
418  } else {
420  // prefer options over loaded location
422  // let offset and boundary lead back to the original coords of the loaded data
425  // the new boundary (updated during loading)
427  }
428  if (lefthand) {
429  myFinal.myOffset.mul(1, -1);
431  }
432 }
433 
434 
435 void
437  myNumLoaded++;
438  if (myNumLoaded > 1) {
439  WRITE_WARNING("Ignoring loaded location attribute nr. " + toString(myNumLoaded) + " for tracking of original location");
440  } else {
441  myLoaded = loaded;
442  }
443 }
444 
445 
446 void
448  myNumLoaded = 0;
449 }
450 
451 
452 void
457  if (myFinal.usingGeoProjection()) {
459  }
461  if (myFinal.usingGeoProjection()) {
462  into.setPrecision();
463  }
465  into.closeTag();
466  into.lf();
467 }
468 
469 
470 
471 /****************************************************************************/
472 
void sub(SUMOReal dx, SUMOReal dy)
Substracts the given position from this one.
Definition: Position.h:139
void doRegister(const std::string &name, Option *v)
Adds an option under the given name.
Definition: OptionsCont.cpp:86
OutputDevice & writeAttr(const SumoXMLAttr attr, const T &val)
writes a named attribute
Definition: OutputDevice.h:257
static void writeLocation(OutputDevice &into)
writes the location element
static GeoConvHelper myProcessing
coordinate transformation to use for input conversion and processing
int getInt(const std::string &name) const
Returns the int-value of the named option (only for Option_Integer)
~GeoConvHelper()
Destructor.
const Boundary & getConvBoundary() const
Returns the converted boundary.
void add(const Position &pos)
Adds the given position to this one.
Definition: Position.h:119
Position myOffset
The offset to apply.
static void computeFinal(bool lefthand=false)
compute the location attributes which will be used for output based on the loaded location data...
#define GEO_OUTPUT_ACCURACY
Definition: config.h:16
bool x2cartesian(Position &from, bool includeInBoundary=true)
bool usingGeoProjection() const
Returns whether a transformation from geo to metric coordinates will be performed.
Boundary myOrigBoundary
The boundary before conversion (x2cartesian)
void setPrecision(int precision=OUTPUT_ACCURACY)
Sets the precison or resets it to default.
static void setLoaded(const GeoConvHelper &loaded)
sets the coordinate transformation loaded from a location element
const Boundary & getOrigBoundary() const
Returns the original boundary.
static GeoConvHelper myLoaded
coordinate transformation loaded from a location element
bool getBool(const std::string &name) const
Returns the boolean-value of the named option (only for Option_Bool)
const std::string & getProjString() const
Returns the network offset.
static void resetLoaded()
resets loaded location elements
bool isDefault(const std::string &name) const
Returns the information whether the named option has still the default value.
bool myUseInverseProjection
Information whether inverse projection shall be used.
A class that stores a 2D geometrical boundary.
Definition: Boundary.h:48
#define WRITE_WARNING(msg)
Definition: MsgHandler.h:200
void addSynonyme(const std::string &name1, const std::string &name2, bool isDeprecated=false)
Adds a synonyme for an options name (any order)
void cartesian2geo(Position &cartesian) const
Converts the given cartesian (shifted) position to its geo (lat/long) representation.
double myGeoScale
The scaling to apply to geo-coordinates.
#define max(a, b)
Definition: polyfonts.c:65
static methods for processing the coordinates conversion for the current net
Definition: GeoConvHelper.h:60
static GeoConvHelper myFinal
coordinate transformation to use for writing the location element and for tracking the original coord...
A point in 2D or 3D with translation and scaling methods.
Definition: Position.h:46
ProjectionMethod myProjectionMethod
Information whether no projection shall be done.
std::string myProjString
A proj options string describing the proj.4-projection to use.
std::string getString(const std::string &name) const
Returns the string-value of the named option (only for Option_String)
bool exists(const std::string &name) const
Returns the information whether the named option is known.
static bool init(OptionsCont &oc)
Initialises the processing and the final instance using the given options.
SUMOReal x() const
Returns the x-position.
Definition: Position.h:63
void addOptionSubTopic(const std::string &topic)
Adds an option subtopic.
#define DEG2RAD(x)
Definition: GeomHelper.h:45
std::string toString(const T &t, std::streamsize accuracy=OUTPUT_ACCURACY)
Definition: ToString.h:55
static void addProjectionOptions(OptionsCont &oc)
Adds projection options to the given container.
GeoConvHelper(OptionsCont &oc)
Constructor based on the stored options.
Boundary myConvBoundary
The boundary after conversion (x2cartesian)
#define WRITE_ERROR(msg)
Definition: MsgHandler.h:206
static int myNumLoaded
the numer of coordinate transformations loaded from location elements
An integer-option.
Definition: Option.h:313
void add(SUMOReal x, SUMOReal y, SUMOReal z=0)
Makes the boundary include the given coordinate.
Definition: Boundary.cpp:90
void flipY()
flips ymin and ymax
Definition: Boundary.cpp:254
A storage for options typed value containers)
Definition: OptionsCont.h:99
void set(SUMOReal x, SUMOReal y)
Definition: Position.h:78
void mul(SUMOReal val)
Multiplies both positions with the given value.
Definition: Position.h:99
bool usingInverseGeoProjection() const
Returns the information whether an inverse transformation will happen.
const Position getOffsetBase() const
Returns the network base.
Static storage of an output device and its base (abstract) implementation.
Definition: OutputDevice.h:71
bool x2cartesian_const(Position &from) const
Converts the given coordinate into a cartesian using the previous initialisation. ...
bool closeTag()
Closes the most recently opened tag.
#define SUMOReal
Definition: config.h:213
const Position getOffset() const
Returns the network offset.
void addDescription(const std::string &name, const std::string &subtopic, const std::string &description)
Adds a description for an option.
#define HAVE_PROJ
Definition: config.h:83
SUMOReal y() const
Returns the y-position.
Definition: Position.h:68
void moveby(SUMOReal x, SUMOReal y, SUMOReal z=0)
Moves the boundary by the given amount.
Definition: Boundary.cpp:281
OutputDevice & openTag(const std::string &xmlElement)
Opens an XML tag.
void lf()
writes a line feed if applicable
Definition: OutputDevice.h:235
void moveConvertedBy(SUMOReal x, SUMOReal y)
Shifts the converted boundary by the given amounts.
GeoConvHelper & operator=(const GeoConvHelper &)
assignment operator.
SUMOReal getFloat(const std::string &name) const
Returns the SUMOReal-value of the named option (only for Option_Float)