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-2017 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 
45 // ===========================================================================
46 // static member variables
47 // ===========================================================================
52 
53 // ===========================================================================
54 // method definitions
55 // ===========================================================================
56 GeoConvHelper::GeoConvHelper(const std::string& proj, const Position& offset,
57  const Boundary& orig, const Boundary& conv, int shift, bool inverse):
58  myProjString(proj),
59 #ifdef HAVE_PROJ
60  myProjection(0),
61  myInverseProjection(0),
62  myGeoProjection(0),
63 #endif
64  myOffset(offset),
65  myGeoScale(pow(10, (double) - shift)),
66  myProjectionMethod(NONE),
67  myUseInverseProjection(inverse),
68  myOrigBoundary(orig),
69  myConvBoundary(conv) {
70  if (proj == "!") {
72  } else if (proj == "-") {
74  } else if (proj == "UTM") {
76  } else if (proj == "DHDN") {
78  } else if (proj == "DHDN_UTM") {
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  if (myInverseProjection != 0) {
99  pj_free(myInverseProjection);
100  }
101  if (myGeoProjection != 0) {
102  pj_free(myInverseProjection);
103  }
104 #endif
105 }
106 
107 
110  myProjString = orig.myProjString;
111  myOffset = orig.myOffset;
115  myGeoScale = orig.myGeoScale;
117 #ifdef HAVE_PROJ
118  if (myProjection != 0) {
119  pj_free(myProjection);
120  myProjection = 0;
121  }
122  if (myInverseProjection != 0) {
123  pj_free(myInverseProjection);
124  myInverseProjection = 0;
125  }
126  if (myGeoProjection != 0) {
127  pj_free(myGeoProjection);
128  myGeoProjection = 0;
129  }
130  if (orig.myProjection != 0) {
131  myProjection = pj_init_plus(orig.myProjString.c_str());
132  }
133  if (orig.myInverseProjection != 0) {
134  myInverseProjection = pj_init_plus(pj_get_def(orig.myInverseProjection, 0));
135  }
136  if (orig.myGeoProjection != 0) {
137  myGeoProjection = pj_init_plus(pj_get_def(orig.myGeoProjection, 0));
138  }
139 #endif
140  return *this;
141 }
142 
143 
144 bool
146  std::string proj = "!"; // the default
147  int shift = oc.getInt("proj.scale");
148  Position offset = Position(oc.getFloat("offset.x"), oc.getFloat("offset.y"));
149  bool inverse = oc.exists("proj.inverse") && oc.getBool("proj.inverse");
150 
151  if (oc.getBool("simple-projection")) {
152  proj = "-";
153  }
154 
155 #ifdef HAVE_PROJ
156  if (oc.getBool("proj.inverse") && oc.getString("proj") == "!") {
157  WRITE_ERROR("Inverse projection works only with explicit proj parameters.");
158  return false;
159  }
160  unsigned numProjections = oc.getBool("simple-projection") + oc.getBool("proj.utm") + oc.getBool("proj.dhdn") + oc.getBool("proj.dhdnutm") + (oc.getString("proj").length() > 1);
161  if (numProjections > 1) {
162  WRITE_ERROR("The projection method needs to be uniquely defined.");
163  return false;
164  }
165 
166  if (oc.getBool("proj.utm")) {
167  proj = "UTM";
168  } else if (oc.getBool("proj.dhdn")) {
169  proj = "DHDN";
170  } else if (oc.getBool("proj.dhdnutm")) {
171  proj = "DHDN_UTM";
172  } else if (!oc.isDefault("proj")) {
173  proj = oc.getString("proj");
174  }
175 #endif
176  myProcessing = GeoConvHelper(proj, offset, Boundary(), Boundary(), shift, inverse);
178  return true;
179 }
180 
181 
182 void
183 GeoConvHelper::init(const std::string& proj,
184  const Position& offset,
185  const Boundary& orig,
186  const Boundary& conv,
187  int shift) {
188  myProcessing = GeoConvHelper(proj, offset, orig, conv, shift);
190 }
191 
192 
193 void
195  oc.addOptionSubTopic("Projection");
196 
197  oc.doRegister("simple-projection", new Option_Bool(false));
198  oc.addSynonyme("simple-projection", "proj.simple", true);
199  oc.addDescription("simple-projection", "Projection", "Uses a simple method for projection");
200 
201  oc.doRegister("proj.scale", new Option_Integer(0));
202  oc.addSynonyme("proj.scale", "proj.shift", true);
203  oc.addDescription("proj.scale", "Projection", "Number of places to shift decimal point to right in geo-coordinates");
204 
205 #ifdef HAVE_PROJ
206  oc.doRegister("proj.utm", new Option_Bool(false));
207  oc.addDescription("proj.utm", "Projection", "Determine the UTM zone (for a universal transversal mercator projection based on the WGS84 ellipsoid)");
208 
209  oc.doRegister("proj.dhdn", new Option_Bool(false));
210  oc.addDescription("proj.dhdn", "Projection", "Determine the DHDN zone (for a transversal mercator projection based on the bessel ellipsoid, \"Gauss-Krueger\")");
211 
212  oc.doRegister("proj", new Option_String("!"));
213  oc.addDescription("proj", "Projection", "Uses STR as proj.4 definition for projection");
214 
215  oc.doRegister("proj.inverse", new Option_Bool(false));
216  oc.addDescription("proj.inverse", "Projection", "Inverses projection");
217 
218  oc.doRegister("proj.dhdnutm", new Option_Bool(false));
219  oc.addDescription("proj.dhdnutm", "Projection", "Convert from Gauss-Krueger to UTM");
220 #endif // HAVE_PROJ
221 }
222 
223 
224 bool
226  return myProjectionMethod != NONE;
227 }
228 
229 
230 bool
232  return myUseInverseProjection;
233 }
234 
235 
236 void
238  cartesian.sub(getOffsetBase());
239  if (myProjectionMethod == NONE) {
240  return;
241  }
242 #ifdef HAVE_PROJ
243  projUV p;
244  p.u = cartesian.x();
245  p.v = cartesian.y();
246  p = pj_inv(p, myProjection);
248  p.u *= RAD_TO_DEG;
249  p.v *= RAD_TO_DEG;
250  cartesian.set((double) p.u, (double) p.v);
251 #endif
252 }
253 
254 
255 bool
256 GeoConvHelper::x2cartesian(Position& from, bool includeInBoundary) {
257  if (includeInBoundary) {
258  myOrigBoundary.add(from);
259  }
260  // init projection parameter on first use
261 #ifdef HAVE_PROJ
262  if (myProjection == 0) {
263  double x = from.x() * myGeoScale;
264  switch (myProjectionMethod) {
265  case DHDN_UTM: {
266  int zone = (int)((x - 500000.) / 1000000.);
267  if (zone < 1 || zone > 5) {
268  WRITE_WARNING("Attempt to initialize DHDN_UTM-projection on invalid longitude " + toString(x));
269  return false;
270  }
271  myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
272  " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
273  " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
274  myInverseProjection = pj_init_plus(myProjString.c_str());
275  myGeoProjection = pj_init_plus("+proj=latlong +datum=WGS84");
277  x = ((x - 500000.) / 1000000.) * 3; // continues with UTM
278  }
279  case UTM: {
280  int zone = (int)(x + 180) / 6 + 1;
281  myProjString = "+proj=utm +zone=" + toString(zone) +
282  " +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
283  myProjection = pj_init_plus(myProjString.c_str());
285  }
286  break;
287  case DHDN: {
288  int zone = (int)(x / 3);
289  if (zone < 1 || zone > 5) {
290  WRITE_WARNING("Attempt to initialize DHDN-projection on invalid longitude " + toString(x));
291  return false;
292  }
293  myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
294  " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
295  " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
296  myProjection = pj_init_plus(myProjString.c_str());
298  }
299  break;
300  default:
301  break;
302  }
303  }
304  if (myInverseProjection != 0) {
305  double x = from.x();
306  double y = from.y();
307  if (pj_transform(myInverseProjection, myGeoProjection, 1, 1, &x, &y, NULL)) {
308  WRITE_WARNING("Could not transform (" + toString(x) + "," + toString(y) + ")");
309  }
310  from.set(double(x * RAD_TO_DEG), double(y * RAD_TO_DEG));
311  }
312 #endif
313  // perform conversion
314  bool ok = x2cartesian_const(from);
315  if (ok) {
316  if (includeInBoundary) {
317  myConvBoundary.add(from);
318  }
319  }
320  return ok;
321 }
322 
323 
324 bool
326  double x = from.x() * myGeoScale;
327  double y = from.y() * myGeoScale;
328  if (myProjectionMethod == NONE) {
329  from.add(myOffset);
330  } else if (myUseInverseProjection) {
331  cartesian2geo(from);
332  } else {
333  if (x > 180.1 || x < -180.1) {
334  WRITE_WARNING("Invalid longitude " + toString(x));
335  return false;
336  }
337  if (y > 90.1 || y < -90.1) {
338  WRITE_WARNING("Invalid latitude " + toString(y));
339  return false;
340  }
341 #ifdef HAVE_PROJ
342  if (myProjection != 0) {
343  projUV p;
344  p.u = x * DEG_TO_RAD;
345  p.v = y * DEG_TO_RAD;
346  p = pj_fwd(p, myProjection);
348  x = p.u;
349  y = p.v;
350  }
351 #endif
352  if (myProjectionMethod == SIMPLE) {
353  double ys = y;
354  x *= 111320. * cos(DEG2RAD(ys));
355  y *= 111136.;
356  from.set((double)x, (double)y);
358  from.add(myOffset);
359  }
360  }
363  return false;
364  }
365  if (myProjectionMethod != SIMPLE) {
366  from.set((double)x, (double)y);
367  from.add(myOffset);
368  }
369  return true;
370 }
371 
372 
373 void
374 GeoConvHelper::moveConvertedBy(double x, double y) {
375  myOffset.add(x, y);
376  myConvBoundary.moveby(x, y);
377 }
378 
379 
380 const Boundary&
382  return myOrigBoundary;
383 }
384 
385 
386 const Boundary&
388  return myConvBoundary;
389 }
390 
391 
392 const Position
394  return myOffset;
395 }
396 
397 
398 const Position
400  return myOffset;
401 }
402 
403 
404 const std::string&
406  return myProjString;
407 }
408 
409 
410 void
412  if (myNumLoaded == 0) {
414  } else {
416  // prefer options over loaded location
418  // let offset and boundary lead back to the original coords of the loaded data
421  // the new boundary (updated during loading)
423  }
424  if (lefthand) {
425  myFinal.myOffset.mul(1, -1);
427  }
428 }
429 
430 
431 void
433  myNumLoaded++;
434  if (myNumLoaded > 1) {
435  WRITE_WARNING("Ignoring loaded location attribute nr. " + toString(myNumLoaded) + " for tracking of original location");
436  } else {
437  myLoaded = loaded;
438  }
439 }
440 
441 
442 void
444  myNumLoaded = 0;
445 }
446 
447 
448 void
453  if (myFinal.usingGeoProjection()) {
455  }
457  if (myFinal.usingGeoProjection()) {
458  into.setPrecision();
459  }
461  into.closeTag();
462  into.lf();
463 }
464 
465 
466 
467 /****************************************************************************/
468 
void doRegister(const std::string &name, Option *v)
Adds an option under the given name.
Definition: OptionsCont.cpp:82
OutputDevice & writeAttr(const SumoXMLAttr attr, const T &val)
writes a named attribute
Definition: OutputDevice.h:256
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:133
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...
bool x2cartesian(Position &from, bool includeInBoundary=true)
double y() const
Returns the y-position.
Definition: Position.h:68
bool usingGeoProjection() const
Returns whether a transformation from geo to metric coordinates will be performed.
void moveby(double x, double y, double z=0)
Moves the boundary by the given amount.
Definition: Boundary.cpp:283
Boundary myOrigBoundary
The boundary before conversion (x2cartesian)
static void setLoaded(const GeoConvHelper &loaded)
sets the coordinate transformation loaded from a location element
double x() const
Returns the x-position.
Definition: Position.h:63
void setPrecision(int precision=gPrecision)
Sets the precison or resets it to default.
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.
void set(double x, double y)
set positions x and y
Definition: Position.h:93
bool myUseInverseProjection
Information whether inverse projection shall be used.
void moveConvertedBy(double x, double y)
Shifts the converted boundary by the given amounts.
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
std::string toString(const T &t, std::streamsize accuracy=gPrecision)
Definition: ToString.h:56
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
int gPrecisionGeo
Definition: StdDefs.cpp:31
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.
void addOptionSubTopic(const std::string &topic)
Adds an option subtopic.
#define DEG2RAD(x)
Definition: GeomHelper.h:45
static void addProjectionOptions(OptionsCont &oc)
Adds projection options to the given container.
GeoConvHelper(OptionsCont &oc)
Constructor based on the stored options.
double getFloat(const std::string &name) const
Returns the double-value of the named option (only for Option_Float)
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 flipY()
flips ymin and ymax
Definition: Boundary.cpp:256
A storage for options typed value containers)
Definition: OptionsCont.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.
void mul(double val)
Multiplies both positions with the given value.
Definition: Position.h:113
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.
void add(double x, double y, double z=0)
Makes the boundary include the given coordinate.
Definition: Boundary.cpp:86
#define HAVE_PROJ
Definition: config.h:74
OutputDevice & openTag(const std::string &xmlElement)
Opens an XML tag.
void lf()
writes a line feed if applicable
Definition: OutputDevice.h:234
GeoConvHelper & operator=(const GeoConvHelper &)
assignment operator.
void sub(double dx, double dy)
Substracts the given position from this one.
Definition: Position.h:153