SUMO - Simulation of Urban MObility
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
NBHeightMapper.cpp
Go to the documentation of this file.
1 /****************************************************************************/
9 // Set z-values for all network positions based on data from a height map
10 /****************************************************************************/
11 // SUMO, Simulation of Urban MObility; see http://sumo.dlr.de/
12 // Copyright (C) 2011-2015 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 <string>
35 #include <utils/common/ToString.h>
38 #include <utils/geom/GeomHelper.h>
39 #include "NBHeightMapper.h"
41 #include <utils/common/RGBColor.h>
43 
44 #ifdef HAVE_GDAL
45 #include <ogrsf_frmts.h>
46 #include <ogr_api.h>
47 #include <gdal_priv.h>
48 #endif
49 
50 #ifdef CHECK_MEMORY_LEAKS
51 #include <foreign/nvwa/debug_new.h>
52 #endif // CHECK_MEMORY_LEAKS
53 
54 // ===========================================================================
55 // static members
56 // ===========================================================================
58 
59 // ===========================================================================
60 // method definitions
61 // ===========================================================================
62 
63 
65  myRTree(&Triangle::addSelf), myRaster(0)
66 { }
67 
68 
70  clearData();
71 }
72 
73 
74 const NBHeightMapper&
76  return Singleton;
77 }
78 
79 
80 bool
82  return myRaster != 0 || myTriangles.size() > 0;
83 }
84 
85 
87 NBHeightMapper::getZ(const Position& geo) const {
88  if (!ready()) {
89  WRITE_WARNING("Cannot supply height since no height data was loaded");
90  return 0;
91  }
92  if (myRaster != 0) {
93  SUMOReal result = -1e6;
94  if (myBoundary.around(geo)) {
95  const int xSize = int((myBoundary.xmax() - myBoundary.xmin()) / mySizeOfPixel.x() + .5);
96  const int iX = int((geo.x() - myBoundary.xmin()) / mySizeOfPixel.x() + .5);
97  const int iY = int((geo.y() - myBoundary.ymax()) / mySizeOfPixel.y() + .5);
98  result = myRaster[iY * xSize + iX];
99  }
100  if (result > -1e5 && result < 1e5) {
101  return result;
102  }
103  }
104  // coordinates in degrees hence a small search window
105  float minB[2];
106  float maxB[2];
107  minB[0] = (float)geo.x() - 0.00001f;
108  minB[1] = (float)geo.y() - 0.00001f;
109  maxB[0] = (float)geo.x() + 0.00001f;
110  maxB[1] = (float)geo.y() + 0.00001f;
111  QueryResult queryResult;
112  int hits = myRTree.Search(minB, maxB, queryResult);
113  Triangles result = queryResult.triangles;
114  assert(hits == (int)result.size());
115  UNUSED_PARAMETER(hits); // only used for assertion
116 
117  for (Triangles::iterator it = result.begin(); it != result.end(); it++) {
118  const Triangle* triangle = *it;
119  if (triangle->contains(geo)) {
120  return triangle->getZ(geo);
121  }
122  }
123  WRITE_WARNING("Could not get height data for coordinate " + toString(geo));
124  return 0;
125 }
126 
127 
128 void
130  Triangle* triangle = new Triangle(corners);
131  myTriangles.push_back(triangle);
132  Boundary b = corners.getBoxBoundary();
133  const float cmin[2] = {(float) b.xmin(), (float) b.ymin()};
134  const float cmax[2] = {(float) b.xmax(), (float) b.ymax()};
135  myRTree.Insert(cmin, cmax, triangle);
136 }
137 
138 
139 void
141  if (oc.isSet("heightmap.geotiff")) {
142  // parse file(s)
143  std::vector<std::string> files = oc.getStringVector("heightmap.geotiff");
144  for (std::vector<std::string>::const_iterator file = files.begin(); file != files.end(); ++file) {
145  PROGRESS_BEGIN_MESSAGE("Parsing from GeoTIFF '" + *file + "'");
146  int numFeatures = Singleton.loadTiff(*file);
148  " done (parsed " + toString(numFeatures) +
149  " features, Boundary: " + toString(Singleton.getBoundary()) + ").");
150  }
151  }
152  if (oc.isSet("heightmap.shapefiles")) {
153  // parse file(s)
154  std::vector<std::string> files = oc.getStringVector("heightmap.shapefiles");
155  for (std::vector<std::string>::const_iterator file = files.begin(); file != files.end(); ++file) {
156  PROGRESS_BEGIN_MESSAGE("Parsing from shape-file '" + *file + "'");
157  int numFeatures = Singleton.loadShapeFile(*file);
159  " done (parsed " + toString(numFeatures) +
160  " features, Boundary: " + toString(Singleton.getBoundary()) + ").");
161  }
162  }
163 }
164 
165 
166 int
167 NBHeightMapper::loadShapeFile(const std::string& file) {
168 #ifdef HAVE_GDAL
169  OGRRegisterAll();
170  OGRDataSource* ds = OGRSFDriverRegistrar::Open(file.c_str(), FALSE);
171  if (ds == NULL) {
172  throw ProcessError("Could not open shape file '" + file + "'.");
173  }
174 
175  // begin file parsing
176  OGRLayer* layer = ds->GetLayer(0);
177  layer->ResetReading();
178 
179  // triangle coordinates are stored in WGS84 and later matched with network coordinates in WGS84
180  // build coordinate transformation
181  OGRSpatialReference* sr_src = layer->GetSpatialRef();
182  OGRSpatialReference sr_dest;
183  sr_dest.SetWellKnownGeogCS("WGS84");
184  OGRCoordinateTransformation* toWGS84 = OGRCreateCoordinateTransformation(sr_src, &sr_dest);
185  if (toWGS84 == 0) {
186  WRITE_WARNING("Could not create geocoordinates converter; check whether proj.4 is installed.");
187  }
188 
189  int numFeatures = 0;
190  OGRFeature* feature;
191  layer->ResetReading();
192  while ((feature = layer->GetNextFeature()) != NULL) {
193  OGRGeometry* geom = feature->GetGeometryRef();
194  assert(geom != 0);
195 
196  // @todo gracefull handling of shapefiles with unexpected contents or any error handling for that matter
197  assert(std::string(geom->getGeometryName()) == std::string("POLYGON"));
198  // try transform to wgs84
199  geom->transform(toWGS84);
200  OGRLinearRing* cgeom = ((OGRPolygon*) geom)->getExteriorRing();
201  // assume TIN with with 4 points and point0 == point3
202  assert(cgeom->getNumPoints() == 4);
203  PositionVector corners;
204  for (int j = 0; j < 3; j++) {
205  Position pos((SUMOReal) cgeom->getX(j), (SUMOReal) cgeom->getY(j), (SUMOReal) cgeom->getZ(j));
206  corners.push_back(pos);
207  myBoundary.add(pos);
208  }
209  addTriangle(corners);
210  numFeatures++;
211 
212  /*
213  OGRwkbGeometryType gtype = geom->getGeometryType();
214  switch (gtype) {
215  case wkbPolygon: {
216  break;
217  }
218  case wkbPoint: {
219  WRITE_WARNING("got wkbPoint");
220  break;
221  }
222  case wkbLineString: {
223  WRITE_WARNING("got wkbLineString");
224  break;
225  }
226  case wkbMultiPoint: {
227  WRITE_WARNING("got wkbMultiPoint");
228  break;
229  }
230  case wkbMultiLineString: {
231  WRITE_WARNING("got wkbMultiLineString");
232  break;
233  }
234  case wkbMultiPolygon: {
235  WRITE_WARNING("got wkbMultiPolygon");
236  break;
237  }
238  default:
239  WRITE_WARNING("Unsupported shape type occured");
240  break;
241  }
242  */
243  OGRFeature::DestroyFeature(feature);
244  }
245  OGRDataSource::DestroyDataSource(ds);
246  OCTDestroyCoordinateTransformation(toWGS84);
247  OGRCleanupAll();
248  return numFeatures;
249 #else
250  WRITE_ERROR("Cannot load shape file since SUMO was compiled without GDAL support.");
251  return 0;
252 #endif
253 }
254 
255 
256 int
257 NBHeightMapper::loadTiff(const std::string& file) {
258 #ifdef HAVE_GDAL
259  GDALAllRegister();
260  GDALDataset* poDataset = (GDALDataset*)GDALOpen(file.c_str(), GA_ReadOnly);
261  if (poDataset == 0) {
262  WRITE_ERROR("Cannot load GeoTIFF file.");
263  return 0;
264  }
265  const int xSize = poDataset->GetRasterXSize();
266  const int ySize = poDataset->GetRasterYSize();
267  double adfGeoTransform[6];
268  if (poDataset->GetGeoTransform(adfGeoTransform) == CE_None) {
269  Position topLeft(adfGeoTransform[0], adfGeoTransform[3]);
270  mySizeOfPixel.set(adfGeoTransform[1], adfGeoTransform[5]);
271  const double horizontalSize = xSize * mySizeOfPixel.x();
272  const double verticalSize = ySize * mySizeOfPixel.y();
273  myBoundary.add(topLeft);
274  myBoundary.add(topLeft.x() + horizontalSize, topLeft.y() + verticalSize);
275  } else {
276  WRITE_ERROR("Could not parse geo information from " + file + ".");
277  return 0;
278  }
279  const int picSize = xSize * ySize;
280  myRaster = (int16_t*)CPLMalloc(sizeof(int16_t) * picSize);
281  for (int i = 1; i <= poDataset->GetRasterCount(); i++) {
282  GDALRasterBand* poBand = poDataset->GetRasterBand(i);
283  if (poBand->GetColorInterpretation() != GCI_GrayIndex) {
284  WRITE_ERROR("Unknown color band in " + file + ".");
285  clearData();
286  break;
287  }
288  if (poBand->GetRasterDataType() != GDT_Int16) {
289  WRITE_ERROR("Unknown data type in " + file + ".");
290  clearData();
291  break;
292  }
293  assert(xSize == poBand->GetXSize() && ySize == poBand->GetYSize());
294  if (poBand->RasterIO(GF_Read, 0, 0, xSize, ySize, myRaster, xSize, ySize, GDT_Int16, 0, 0) == CE_Failure) {
295  WRITE_ERROR("Failure in reading " + file + ".");
296  clearData();
297  break;
298  }
299  }
300  GDALClose(poDataset);
301  return picSize;
302 #else
303  WRITE_ERROR("Cannot load GeoTIFF file since SUMO was compiled without GDAL support.");
304  return 0;
305 #endif
306 }
307 
308 
309 void
311  for (Triangles::iterator it = myTriangles.begin(); it != myTriangles.end(); it++) {
312  delete *it;
313  }
314  myTriangles.clear();
315 #ifdef HAVE_GDAL
316  if (myRaster != 0) {
317  CPLFree(myRaster);
318  myRaster = 0;
319  }
320 #endif
321  myBoundary.reset();
322 }
323 
324 
325 // ===========================================================================
326 // Triangle member methods
327 // ===========================================================================
329  myCorners(corners) {
330  assert(myCorners.size() == 3);
331  // @todo assert non-colinearity
332 }
333 
334 
335 void
337  queryResult.triangles.push_back(this);
338 }
339 
340 
341 bool
343  return myCorners.around(pos);
344 }
345 
346 
347 SUMOReal
349  // en.wikipedia.org/wiki/Line-plane_intersection
350  Position p0 = myCorners.front();
351  Position line(0, 0, 1);
352  p0.sub(geo); // p0 - l0
353  Position normal = normalVector();
354  return p0.dotProduct(normal) / line.dotProduct(normal);
355 }
356 
357 
358 Position
360  // @todo maybe cache result to avoid multiple computations?
361  Position side1 = myCorners[1] - myCorners[0];
362  Position side2 = myCorners[2] - myCorners[0];
363  return side1.crossProduct(side2);
364 }
365 
366 
367 /****************************************************************************/
368 
void sub(SUMOReal dx, SUMOReal dy)
Substracts the given position from this one.
Definition: Position.h:139
std::vector< std::string > getStringVector(const std::string &name) const
Returns the list of string-vector-value of the named option (only for Option_String) ...
SUMOReal getZ(const Position &geo) const
returns the projection of the give geoCoordinate (WGS84) onto triangle plane
Triangles myTriangles
bool ready() const
returns whether the NBHeightMapper has data
bool around(const Position &p, SUMOReal offset=0) const
Returns whether the boundary contains the given coordinate.
Definition: Boundary.cpp:148
Position normalVector() const
returns the normal vector for this triangles plane
int loadShapeFile(const std::string &file)
load height data from Arcgis-shape file and returns the number of parsed features ...
SUMOReal getZ(const Position &geo) const
returns height for the given geo coordinate (WGS84)
SUMOReal ymin() const
Returns minimum y-coordinate.
Definition: Boundary.cpp:124
SUMOReal xmin() const
Returns minimum x-coordinate.
Definition: Boundary.cpp:112
Triangle(const PositionVector &corners)
void addTriangle(PositionVector corners)
adds one triangles worth of height data
void clearData()
clears loaded data
static NBHeightMapper Singleton
the singleton instance
NBHeightMapper()
private constructor and destructor (Singleton)
SUMOReal x() const
Returns the x-position.
Definition: Position.h:63
#define UNUSED_PARAMETER(x)
Definition: StdDefs.h:39
SUMOReal xmax() const
Returns maximum x-coordinate.
Definition: Boundary.cpp:118
A class that stores a 2D geometrical boundary.
Definition: Boundary.h:48
#define WRITE_WARNING(msg)
Definition: MsgHandler.h:200
int loadTiff(const std::string &file)
load height data from GeoTIFF file and returns the number of non void pixels
static const NBHeightMapper & get()
return the singleton instance (maybe 0)
class for cirumventing the const-restriction of RTree::Search-context
static void loadIfSet(OptionsCont &oc)
loads heigh map data if any loading options are set
TRIANGLE_RTREE_QUAL myRTree
The RTree for spatial queries.
Position crossProduct(const Position &pos)
returns the cross product between this point and the second one
Definition: Position.h:250
A point in 2D or 3D with translation and scaling methods.
Definition: Position.h:46
A list of positions.
std::vector< const Triangle * > Triangles
#define PROGRESS_BEGIN_MESSAGE(msg)
Definition: MsgHandler.h:202
static MsgHandler * getMessageInstance()
Returns the instance to add normal messages to.
Definition: MsgHandler.cpp:62
int16_t * myRaster
raster height information in m
bool contains(const Position &pos) const
checks whether pos lies within triangle (only checks x,y)
std::string toString(const T &t, std::streamsize accuracy=OUTPUT_ACCURACY)
Definition: ToString.h:53
Position mySizeOfPixel
dimensions of one pixel in raster data
Boundary myBoundary
convex boundary of all known triangles;
#define WRITE_ERROR(msg)
Definition: MsgHandler.h:205
void add(SUMOReal x, SUMOReal y)
Makes the boundary include the given coordinate.
Definition: Boundary.cpp:76
void reset()
Resets the boundary.
Definition: Boundary.cpp:66
SUMOReal dotProduct(const Position &pos)
returns the dot product (scalar product) between this point and the second one
Definition: Position.h:258
SUMOReal y() const
Returns the y-position.
Definition: Position.h:68
PositionVector myCorners
the corners of the triangle
A storage for options typed value containers)
Definition: OptionsCont.h:108
void set(SUMOReal x, SUMOReal y)
Definition: Position.h:78
void addSelf(const QueryResult &queryResult) const
callback for RTree search
#define SUMOReal
Definition: config.h:218
SUMOReal ymax() const
Returns maximum y-coordinate.
Definition: Boundary.cpp:130
const Boundary & getBoundary()
returns the convex boundary of all known triangles
Boundary getBoxBoundary() const
Returns a boundary enclosing this list of lines.
Set z-values for all network positions based on data from a height map.
bool isSet(const std::string &name, bool failOnNonExistant=true) const
Returns the information whether the named option is set.
void endProcessMsg(std::string msg)
Ends a process information.
Definition: MsgHandler.cpp:131