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