SUMO - Simulation of Urban MObility
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>
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 #if GDAL_VERSION_MAJOR < 2
170  OGRRegisterAll();
171  OGRDataSource* ds = OGRSFDriverRegistrar::Open(file.c_str(), FALSE);
172 #else
173  GDALAllRegister();
174  GDALDataset* ds = (GDALDataset*) GDALOpen(file.c_str(), GA_ReadOnly);
175 #endif
176  if (ds == NULL) {
177  throw ProcessError("Could not open shape file '" + file + "'.");
178  }
179 
180  // begin file parsing
181  OGRLayer* layer = ds->GetLayer(0);
182  layer->ResetReading();
183 
184  // triangle coordinates are stored in WGS84 and later matched with network coordinates in WGS84
185  // build coordinate transformation
186  OGRSpatialReference* sr_src = layer->GetSpatialRef();
187  OGRSpatialReference sr_dest;
188  sr_dest.SetWellKnownGeogCS("WGS84");
189  OGRCoordinateTransformation* toWGS84 = OGRCreateCoordinateTransformation(sr_src, &sr_dest);
190  if (toWGS84 == 0) {
191  WRITE_WARNING("Could not create geocoordinates converter; check whether proj.4 is installed.");
192  }
193 
194  int numFeatures = 0;
195  OGRFeature* feature;
196  layer->ResetReading();
197  while ((feature = layer->GetNextFeature()) != NULL) {
198  OGRGeometry* geom = feature->GetGeometryRef();
199  assert(geom != 0);
200 
201  // @todo gracefull handling of shapefiles with unexpected contents or any error handling for that matter
202  assert(std::string(geom->getGeometryName()) == std::string("POLYGON"));
203  // try transform to wgs84
204  geom->transform(toWGS84);
205  OGRLinearRing* cgeom = ((OGRPolygon*) geom)->getExteriorRing();
206  // assume TIN with with 4 points and point0 == point3
207  assert(cgeom->getNumPoints() == 4);
208  PositionVector corners;
209  for (int j = 0; j < 3; j++) {
210  Position pos((SUMOReal) cgeom->getX(j), (SUMOReal) cgeom->getY(j), (SUMOReal) cgeom->getZ(j));
211  corners.push_back(pos);
212  myBoundary.add(pos);
213  }
214  addTriangle(corners);
215  numFeatures++;
216 
217  /*
218  OGRwkbGeometryType gtype = geom->getGeometryType();
219  switch (gtype) {
220  case wkbPolygon: {
221  break;
222  }
223  case wkbPoint: {
224  WRITE_WARNING("got wkbPoint");
225  break;
226  }
227  case wkbLineString: {
228  WRITE_WARNING("got wkbLineString");
229  break;
230  }
231  case wkbMultiPoint: {
232  WRITE_WARNING("got wkbMultiPoint");
233  break;
234  }
235  case wkbMultiLineString: {
236  WRITE_WARNING("got wkbMultiLineString");
237  break;
238  }
239  case wkbMultiPolygon: {
240  WRITE_WARNING("got wkbMultiPolygon");
241  break;
242  }
243  default:
244  WRITE_WARNING("Unsupported shape type occured");
245  break;
246  }
247  */
248  OGRFeature::DestroyFeature(feature);
249  }
250 #if GDAL_VERSION_MAJOR < 2
251  OGRDataSource::DestroyDataSource(ds);
252 #else
253  GDALClose(ds);
254 #endif
255  OCTDestroyCoordinateTransformation(toWGS84);
256  OGRCleanupAll();
257  return numFeatures;
258 #else
259  WRITE_ERROR("Cannot load shape file since SUMO was compiled without GDAL support.");
260  return 0;
261 #endif
262 }
263 
264 
265 int
266 NBHeightMapper::loadTiff(const std::string& file) {
267 #ifdef HAVE_GDAL
268  GDALAllRegister();
269  GDALDataset* poDataset = (GDALDataset*)GDALOpen(file.c_str(), GA_ReadOnly);
270  if (poDataset == 0) {
271  WRITE_ERROR("Cannot load GeoTIFF file.");
272  return 0;
273  }
274  const int xSize = poDataset->GetRasterXSize();
275  const int ySize = poDataset->GetRasterYSize();
276  double adfGeoTransform[6];
277  if (poDataset->GetGeoTransform(adfGeoTransform) == CE_None) {
278  Position topLeft(adfGeoTransform[0], adfGeoTransform[3]);
279  mySizeOfPixel.set(adfGeoTransform[1], adfGeoTransform[5]);
280  const double horizontalSize = xSize * mySizeOfPixel.x();
281  const double verticalSize = ySize * mySizeOfPixel.y();
282  myBoundary.add(topLeft);
283  myBoundary.add(topLeft.x() + horizontalSize, topLeft.y() + verticalSize);
284  } else {
285  WRITE_ERROR("Could not parse geo information from " + file + ".");
286  return 0;
287  }
288  const int picSize = xSize * ySize;
289  myRaster = (int16_t*)CPLMalloc(sizeof(int16_t) * picSize);
290  for (int i = 1; i <= poDataset->GetRasterCount(); i++) {
291  GDALRasterBand* poBand = poDataset->GetRasterBand(i);
292  if (poBand->GetColorInterpretation() != GCI_GrayIndex) {
293  WRITE_ERROR("Unknown color band in " + file + ".");
294  clearData();
295  break;
296  }
297  if (poBand->GetRasterDataType() != GDT_Int16) {
298  WRITE_ERROR("Unknown data type in " + file + ".");
299  clearData();
300  break;
301  }
302  assert(xSize == poBand->GetXSize() && ySize == poBand->GetYSize());
303  if (poBand->RasterIO(GF_Read, 0, 0, xSize, ySize, myRaster, xSize, ySize, GDT_Int16, 0, 0) == CE_Failure) {
304  WRITE_ERROR("Failure in reading " + file + ".");
305  clearData();
306  break;
307  }
308  }
309  GDALClose(poDataset);
310  return picSize;
311 #else
312  WRITE_ERROR("Cannot load GeoTIFF file since SUMO was compiled without GDAL support.");
313  return 0;
314 #endif
315 }
316 
317 
318 void
320  for (Triangles::iterator it = myTriangles.begin(); it != myTriangles.end(); it++) {
321  delete *it;
322  }
323  myTriangles.clear();
324 #ifdef HAVE_GDAL
325  if (myRaster != 0) {
326  CPLFree(myRaster);
327  myRaster = 0;
328  }
329 #endif
330  myBoundary.reset();
331 }
332 
333 
334 // ===========================================================================
335 // Triangle member methods
336 // ===========================================================================
338  myCorners(corners) {
339  assert(myCorners.size() == 3);
340  // @todo assert non-colinearity
341 }
342 
343 
344 void
346  queryResult.triangles.push_back(this);
347 }
348 
349 
350 bool
352  return myCorners.around(pos);
353 }
354 
355 
356 SUMOReal
358  // en.wikipedia.org/wiki/Line-plane_intersection
359  Position p0 = myCorners.front();
360  Position line(0, 0, 1);
361  p0.sub(geo); // p0 - l0
362  Position normal = normalVector();
363  return p0.dotProduct(normal) / line.dotProduct(normal);
364 }
365 
366 
367 Position
369  // @todo maybe cache result to avoid multiple computations?
370  Position side1 = myCorners[1] - myCorners[0];
371  Position side2 = myCorners[2] - myCorners[0];
372  return side1.crossProduct(side2);
373 }
374 
375 
376 /****************************************************************************/
377 
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
bool around(const Position &p, SUMOReal offset=0) const
Returns the information whether the position vector describes a polygon lying around the given point ...
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: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:54
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 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: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: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:213
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