SUMO - Simulation of Urban MObility
GeoConvHelper.cpp
Go to the documentation of this file.
1 /****************************************************************************/
2 // Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.org/sumo
3 // Copyright (C) 2001-2017 German Aerospace Center (DLR) and others.
4 /****************************************************************************/
5 //
6 // This program and the accompanying materials
7 // are made available under the terms of the Eclipse Public License v2.0
8 // which accompanies this distribution, and is available at
9 // http://www.eclipse.org/legal/epl-v20.html
10 //
11 /****************************************************************************/
19 // static methods for processing the coordinates conversion for the current net
20 /****************************************************************************/
21 
22 
23 // ===========================================================================
24 // included modules
25 // ===========================================================================
26 #ifdef _MSC_VER
27 #include <windows_config.h>
28 #else
29 #include <config.h>
30 #endif
31 
32 #include <map>
33 #include <cmath>
34 #include <cassert>
35 #include <climits>
37 #include <utils/common/ToString.h>
38 #include <utils/geom/GeomHelper.h>
41 #include "GeoConvHelper.h"
42 
43 
44 // ===========================================================================
45 // static member variables
46 // ===========================================================================
47 
52 
53 // ===========================================================================
54 // method definitions
55 // ===========================================================================
56 
57 GeoConvHelper::GeoConvHelper(const std::string& proj, const Position& offset,
58  const Boundary& orig, const Boundary& conv, double scale, bool inverse):
59  myProjString(proj),
60 #ifdef HAVE_PROJ
61  myProjection(0),
62  myInverseProjection(0),
63  myGeoProjection(0),
64 #endif
65  myOffset(offset),
66  myGeoScale(scale),
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  } else if (proj == "DHDN_UTM") {
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  if (myInverseProjection != 0) {
100  pj_free(myInverseProjection);
101  }
102  if (myGeoProjection != 0) {
103  pj_free(myInverseProjection);
104  }
105 #endif
106 }
107 
108 
111  myProjString = orig.myProjString;
112  myOffset = orig.myOffset;
116  myGeoScale = orig.myGeoScale;
118 #ifdef HAVE_PROJ
119  if (myProjection != 0) {
120  pj_free(myProjection);
121  myProjection = 0;
122  }
123  if (myInverseProjection != 0) {
124  pj_free(myInverseProjection);
125  myInverseProjection = 0;
126  }
127  if (myGeoProjection != 0) {
128  pj_free(myGeoProjection);
129  myGeoProjection = 0;
130  }
131  if (orig.myProjection != 0) {
132  myProjection = pj_init_plus(orig.myProjString.c_str());
133  }
134  if (orig.myInverseProjection != 0) {
135  myInverseProjection = pj_init_plus(pj_get_def(orig.myInverseProjection, 0));
136  }
137  if (orig.myGeoProjection != 0) {
138  myGeoProjection = pj_init_plus(pj_get_def(orig.myGeoProjection, 0));
139  }
140 #endif
141  return *this;
142 }
143 
144 
145 bool
147  std::string proj = "!"; // the default
148  double scale = oc.getFloat("proj.scale");
149  Position offset = Position(oc.getFloat("offset.x"), oc.getFloat("offset.y"));
150  bool inverse = oc.exists("proj.inverse") && oc.getBool("proj.inverse");
151 
152  if (oc.getBool("simple-projection")) {
153  proj = "-";
154  }
155 
156 #ifdef HAVE_PROJ
157  if (oc.getBool("proj.inverse") && oc.getString("proj") == "!") {
158  WRITE_ERROR("Inverse projection works only with explicit proj parameters.");
159  return false;
160  }
161  unsigned numProjections = oc.getBool("simple-projection") + oc.getBool("proj.utm") + oc.getBool("proj.dhdn") + oc.getBool("proj.dhdnutm") + (oc.getString("proj").length() > 1);
162  if (numProjections > 1) {
163  WRITE_ERROR("The projection method needs to be uniquely defined.");
164  return false;
165  }
166 
167  if (oc.getBool("proj.utm")) {
168  proj = "UTM";
169  } else if (oc.getBool("proj.dhdn")) {
170  proj = "DHDN";
171  } else if (oc.getBool("proj.dhdnutm")) {
172  proj = "DHDN_UTM";
173  } else if (!oc.isDefault("proj")) {
174  proj = oc.getString("proj");
175  }
176 #endif
177  myProcessing = GeoConvHelper(proj, offset, Boundary(), Boundary(), scale, inverse);
179  return true;
180 }
181 
182 
183 void
184 GeoConvHelper::init(const std::string& proj, const Position& offset, const Boundary& orig,
185  const Boundary& conv, double scale) {
186  myProcessing = GeoConvHelper(proj, offset, orig, conv, scale);
188 }
189 
190 
191 void
193  oc.addOptionSubTopic("Projection");
194 
195  oc.doRegister("simple-projection", new Option_Bool(false));
196  oc.addSynonyme("simple-projection", "proj.simple", true);
197  oc.addDescription("simple-projection", "Projection", "Uses a simple method for projection");
198 
199  oc.doRegister("proj.scale", new Option_Float(1.0));
200  oc.addDescription("proj.scale", "Projection", "Scaling factor for input coordinates");
201 
202 #ifdef HAVE_PROJ
203  oc.doRegister("proj.utm", new Option_Bool(false));
204  oc.addDescription("proj.utm", "Projection", "Determine the UTM zone (for a universal transversal mercator projection based on the WGS84 ellipsoid)");
205 
206  oc.doRegister("proj.dhdn", new Option_Bool(false));
207  oc.addDescription("proj.dhdn", "Projection", "Determine the DHDN zone (for a transversal mercator projection based on the bessel ellipsoid, \"Gauss-Krueger\")");
208 
209  oc.doRegister("proj", new Option_String("!"));
210  oc.addDescription("proj", "Projection", "Uses STR as proj.4 definition for projection");
211 
212  oc.doRegister("proj.inverse", new Option_Bool(false));
213  oc.addDescription("proj.inverse", "Projection", "Inverses projection");
214 
215  oc.doRegister("proj.dhdnutm", new Option_Bool(false));
216  oc.addDescription("proj.dhdnutm", "Projection", "Convert from Gauss-Krueger to UTM");
217 #endif // HAVE_PROJ
218 }
219 
220 
221 bool
223  return myProjectionMethod != NONE;
224 }
225 
226 
227 bool
229  return myUseInverseProjection;
230 }
231 
232 
233 void
235  cartesian.sub(getOffsetBase());
236  if (myProjectionMethod == NONE) {
237  return;
238  }
239  if (myProjectionMethod == SIMPLE) {
240  const double y = cartesian.y() / 111136.;
241  const double x = cartesian.x() / 111320. / cos(DEG2RAD(y));
242  cartesian.set(x, y);
243  return;
244  }
245 #ifdef HAVE_PROJ
246  projUV p;
247  p.u = cartesian.x();
248  p.v = cartesian.y();
249  p = pj_inv(p, myProjection);
251  p.u *= RAD_TO_DEG;
252  p.v *= RAD_TO_DEG;
253  cartesian.set((double) p.u, (double) p.v);
254 #endif
255 }
256 
257 
258 bool
259 GeoConvHelper::x2cartesian(Position& from, bool includeInBoundary) {
260  if (includeInBoundary) {
261  myOrigBoundary.add(from);
262  }
263  // init projection parameter on first use
264 #ifdef HAVE_PROJ
265  if (myProjection == 0) {
266  double x = from.x() * myGeoScale;
267  switch (myProjectionMethod) {
268  case DHDN_UTM: {
269  int zone = (int)((x - 500000.) / 1000000.);
270  if (zone < 1 || zone > 5) {
271  WRITE_WARNING("Attempt to initialize DHDN_UTM-projection on invalid longitude " + toString(x));
272  return false;
273  }
274  myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
275  " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
276  " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
277  myInverseProjection = pj_init_plus(myProjString.c_str());
278  myGeoProjection = pj_init_plus("+proj=latlong +datum=WGS84");
280  x = ((x - 500000.) / 1000000.) * 3; // continues with UTM
281  }
282  case UTM: {
283  int zone = (int)(x + 180) / 6 + 1;
284  myProjString = "+proj=utm +zone=" + toString(zone) +
285  " +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
286  myProjection = pj_init_plus(myProjString.c_str());
288  }
289  break;
290  case DHDN: {
291  int zone = (int)(x / 3);
292  if (zone < 1 || zone > 5) {
293  WRITE_WARNING("Attempt to initialize DHDN-projection on invalid longitude " + toString(x));
294  return false;
295  }
296  myProjString = "+proj=tmerc +lat_0=0 +lon_0=" + toString(3 * zone) +
297  " +k=1 +x_0=" + toString(zone * 1000000 + 500000) +
298  " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
299  myProjection = pj_init_plus(myProjString.c_str());
301  }
302  break;
303  default:
304  break;
305  }
306  }
307  if (myInverseProjection != 0) {
308  double x = from.x();
309  double y = from.y();
310  if (pj_transform(myInverseProjection, myGeoProjection, 1, 1, &x, &y, NULL)) {
311  WRITE_WARNING("Could not transform (" + toString(x) + "," + toString(y) + ")");
312  }
313  from.set(double(x * RAD_TO_DEG), double(y * RAD_TO_DEG));
314  }
315 #endif
316  // perform conversion
317  bool ok = x2cartesian_const(from);
318  if (ok) {
319  if (includeInBoundary) {
320  myConvBoundary.add(from);
321  }
322  }
323  return ok;
324 }
325 
326 
327 bool
329  double x = from.x() * myGeoScale;
330  double y = from.y() * myGeoScale;
331  if (myProjectionMethod == NONE) {
332  from.add(myOffset);
333  } else if (myUseInverseProjection) {
334  cartesian2geo(from);
335  } else {
336  if (x > 180.1 || x < -180.1) {
337  WRITE_WARNING("Invalid longitude " + toString(x));
338  return false;
339  }
340  if (y > 90.1 || y < -90.1) {
341  WRITE_WARNING("Invalid latitude " + toString(y));
342  return false;
343  }
344 #ifdef HAVE_PROJ
345  if (myProjection != 0) {
346  projUV p;
347  p.u = x * DEG_TO_RAD;
348  p.v = y * DEG_TO_RAD;
349  p = pj_fwd(p, myProjection);
351  x = p.u;
352  y = p.v;
353  }
354 #endif
355  if (myProjectionMethod == SIMPLE) {
356  x *= 111320. * cos(DEG2RAD(y));
357  y *= 111136.;
359  }
360  }
361  if (x > std::numeric_limits<double>::max() ||
362  y > std::numeric_limits<double>::max()) {
363  return false;
364  }
365  from.set(x, y);
366  from.add(myOffset);
367  return true;
368 }
369 
370 
371 void
372 GeoConvHelper::moveConvertedBy(double x, double y) {
373  myOffset.add(x, y);
374  myConvBoundary.moveby(x, y);
375 }
376 
377 
378 const Boundary&
380  return myOrigBoundary;
381 }
382 
383 
384 const Boundary&
386  return myConvBoundary;
387 }
388 
389 
390 const Position
392  return myOffset;
393 }
394 
395 
396 const Position
398  return myOffset;
399 }
400 
401 
402 const std::string&
404  return myProjString;
405 }
406 
407 
408 void
410  if (myNumLoaded == 0) {
412  if (lefthand) {
413  myFinal.myOffset.mul(1, -1);
414  }
415  } else {
416  if (lefthand) {
417  myProcessing.myOffset.mul(1, -1);
418  }
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) {
430  }
431 }
432 
433 
434 void
436  myNumLoaded++;
437  if (myNumLoaded > 1) {
438  WRITE_WARNING("Ignoring loaded location attribute nr. " + toString(myNumLoaded) + " for tracking of original location");
439  } else {
440  myLoaded = loaded;
441  }
442 }
443 
444 
445 void
447  myNumLoaded = 0;
448 }
449 
450 
451 void
456  if (myFinal.usingGeoProjection()) {
458  }
460  if (myFinal.usingGeoProjection()) {
461  into.setPrecision();
462  }
464  into.closeTag();
465  into.lf();
466 }
467 
468 
469 
470 /****************************************************************************/
471 
void doRegister(const std::string &name, Option *v)
Adds an option under the given name.
Definition: OptionsCont.cpp:81
OutputDevice & writeAttr(const SumoXMLAttr attr, const T &val)
writes a named attribute
Definition: OutputDevice.h:260
static void writeLocation(OutputDevice &into)
writes the location element
static GeoConvHelper myProcessing
coordinate transformation to use for input conversion and processing
~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:132
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)
Converts the given coordinate into a cartesian and optionally update myConvBoundary.
double y() const
Returns the y-position.
Definition: Position.h:67
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:350
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:62
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:92
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:47
#define WRITE_WARNING(msg)
Definition: MsgHandler.h:199
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.
static methods for processing the coordinates conversion for the current net
Definition: GeoConvHelper.h:59
std::string toString(const T &t, std::streamsize accuracy=gPrecision)
Definition: ToString.h:55
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:45
int gPrecisionGeo
Definition: StdDefs.cpp:30
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:44
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:205
static int myNumLoaded
the numer of coordinate transformations loaded from location elements
void flipY()
flips ymin and ymax
Definition: Boundary.cpp:323
A storage for options typed value containers)
Definition: OptionsCont.h:98
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:70
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:112
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:85
#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:238
GeoConvHelper & operator=(const GeoConvHelper &)
make assignment operator private
void sub(double dx, double dy)
Substracts the given position from this one.
Definition: Position.h:152