Eclipse SUMO - Simulation of Urban MObility
Loading...
Searching...
No Matches
NBHeightMapper.cpp
Go to the documentation of this file.
1/****************************************************************************/
2// Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.dev/sumo
3// Copyright (C) 2011-2024 German Aerospace Center (DLR) and others.
4// This program and the accompanying materials are made available under the
5// terms of the Eclipse Public License 2.0 which is available at
6// https://www.eclipse.org/legal/epl-2.0/
7// This Source Code may also be made available under the following Secondary
8// Licenses when the conditions for such availability set forth in the Eclipse
9// Public License 2.0 are satisfied: GNU General Public License, version 2
10// or later which is available at
11// https://www.gnu.org/licenses/old-licenses/gpl-2.0-standalone.html
12// SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later
13/****************************************************************************/
20// Set z-values for all network positions based on data from a height map
21/****************************************************************************/
22#include <config.h>
23
24#include <string>
30#include "NBHeightMapper.h"
33
34#ifdef HAVE_GDAL
35#ifdef _MSC_VER
36#pragma warning(push)
37#pragma warning(disable: 4435 5219 5220)
38#endif
39#if __GNUC__ > 3
40#pragma GCC diagnostic push
41#pragma GCC diagnostic ignored "-Wpedantic"
42#endif
43#include <ogrsf_frmts.h>
44#include <ogr_api.h>
45#include <gdal_priv.h>
46#if __GNUC__ > 3
47#pragma GCC diagnostic pop
48#endif
49#ifdef _MSC_VER
50#pragma warning(pop)
51#endif
52#endif
53
54// ===========================================================================
55// static members
56// ===========================================================================
58
59
60// ===========================================================================
61// method definitions
62// ===========================================================================
64 myRTree(&Triangle::addSelf) {
65}
66
67
71
72
73const NBHeightMapper&
77
78
79bool
81 return myRasters.size() > 0 || myTriangles.size() > 0;
82}
83
84
85double
86NBHeightMapper::getZ(const Position& geo) const {
87 if (!ready()) {
88 WRITE_WARNING(TL("Cannot supply height since no height data was loaded"));
89 return 0;
90 }
91 for (auto& item : myRasters) {
92 const Boundary& boundary = item.boundary;
93 int16_t* raster = item.raster;
94 double result = -1e6;
95 if (boundary.around(geo)) {
96 const int xSize = item.xSize;
97 const double normX = (geo.x() - boundary.xmin()) / mySizeOfPixel.x();
98 const double normY = (geo.y() - boundary.ymax()) / mySizeOfPixel.y();
99 PositionVector corners;
100 corners.push_back(Position(floor(normX) + 0.5, floor(normY) + 0.5, raster[(int)normY * xSize + (int)normX]));
101 if (normX - floor(normX) > 0.5) {
102 corners.push_back(Position(floor(normX) + 1.5, floor(normY) + 0.5, raster[(int)normY * xSize + (int)normX + 1]));
103 } else {
104 corners.push_back(Position(floor(normX) - 0.5, floor(normY) + 0.5, raster[(int)normY * xSize + (int)normX - 1]));
105 }
106 if (normY - floor(normY) > 0.5 && ((int)normY + 1) < item.ySize) {
107 corners.push_back(Position(floor(normX) + 0.5, floor(normY) + 1.5, raster[((int)normY + 1) * xSize + (int)normX]));
108 } else {
109 corners.push_back(Position(floor(normX) + 0.5, floor(normY) - 0.5, raster[((int)normY - 1) * xSize + (int)normX]));
110 }
111 result = Triangle(corners).getZ(Position(normX, normY));
112 }
113 if (result > -1e5 && result < 1e5) {
114 return result;
115 }
116 }
117 // coordinates in degrees hence a small search window
118 float minB[2];
119 float maxB[2];
120 minB[0] = (float)geo.x() - 0.00001f;
121 minB[1] = (float)geo.y() - 0.00001f;
122 maxB[0] = (float)geo.x() + 0.00001f;
123 maxB[1] = (float)geo.y() + 0.00001f;
124 QueryResult queryResult;
125 int hits = myRTree.Search(minB, maxB, queryResult);
126 Triangles result = queryResult.triangles;
127 assert(hits == (int)result.size());
128 UNUSED_PARAMETER(hits); // only used for assertion
129
130 for (Triangles::iterator it = result.begin(); it != result.end(); it++) {
131 const Triangle* triangle = *it;
132 if (triangle->contains(geo)) {
133 return triangle->getZ(geo);
134 }
135 }
136 WRITE_WARNINGF(TL("Could not get height data for coordinate %"), toString(geo));
137 return 0;
138}
139
140
141void
143 Triangle* triangle = new Triangle(corners);
144 myTriangles.push_back(triangle);
145 Boundary b = corners.getBoxBoundary();
146 const float cmin[2] = {(float) b.xmin(), (float) b.ymin()};
147 const float cmax[2] = {(float) b.xmax(), (float) b.ymax()};
148 myRTree.Insert(cmin, cmax, triangle);
149}
150
151
152void
154 if (oc.isSet("heightmap.geotiff")) {
155 // parse file(s)
156 std::vector<std::string> files = oc.getStringVector("heightmap.geotiff");
157 for (std::vector<std::string>::const_iterator file = files.begin(); file != files.end(); ++file) {
158 PROGRESS_BEGIN_MESSAGE("Parsing from GeoTIFF '" + *file + "'");
159 int numFeatures = myInstance.loadTiff(*file);
161 " done (parsed " + toString(numFeatures) +
162 " features, Boundary: " + toString(myInstance.getBoundary()) + ").");
163 }
164 }
165 if (oc.isSet("heightmap.shapefiles")) {
166 // parse file(s)
167 std::vector<std::string> files = oc.getStringVector("heightmap.shapefiles");
168 for (std::vector<std::string>::const_iterator file = files.begin(); file != files.end(); ++file) {
169 PROGRESS_BEGIN_MESSAGE("Parsing from shape-file '" + *file + "'");
170 int numFeatures = myInstance.loadShapeFile(*file);
172 " done (parsed " + toString(numFeatures) +
173 " features, Boundary: " + toString(myInstance.getBoundary()) + ").");
174 }
175 }
176}
177
178
179int
180NBHeightMapper::loadShapeFile(const std::string& file) {
181#ifdef HAVE_GDAL
182#if GDAL_VERSION_MAJOR < 2
183 OGRRegisterAll();
184 OGRDataSource* ds = OGRSFDriverRegistrar::Open(file.c_str(), FALSE);
185#else
186 GDALAllRegister();
187 GDALDataset* ds = (GDALDataset*)GDALOpenEx(file.c_str(), GDAL_OF_VECTOR | GA_ReadOnly, nullptr, nullptr, nullptr);
188#endif
189 if (ds == nullptr) {
190 throw ProcessError(TLF("Could not open shape file '%'.", file));
191 }
192
193 // begin file parsing
194 OGRLayer* layer = ds->GetLayer(0);
195 layer->ResetReading();
196
197 // triangle coordinates are stored in WGS84 and later matched with network coordinates in WGS84
198 // build coordinate transformation
199 OGRSpatialReference* sr_src = layer->GetSpatialRef();
200 OGRSpatialReference sr_dest;
201 sr_dest.SetWellKnownGeogCS("WGS84");
202 OGRCoordinateTransformation* toWGS84 = OGRCreateCoordinateTransformation(sr_src, &sr_dest);
203 if (toWGS84 == nullptr) {
204 WRITE_WARNING(TL("Could not create geocoordinates converter; check whether proj.4 is installed."));
205 }
206
207 int numFeatures = 0;
208 OGRFeature* feature;
209 layer->ResetReading();
210 while ((feature = layer->GetNextFeature()) != nullptr) {
211 OGRGeometry* geom = feature->GetGeometryRef();
212 assert(geom != 0);
213
214 OGRwkbGeometryType gtype = geom->getGeometryType();
215 if (gtype == wkbPolygon) {
216 assert(std::string(geom->getGeometryName()) == std::string("POLYGON"));
217 // try transform to wgs84
218 geom->transform(toWGS84);
219 OGRLinearRing* cgeom = ((OGRPolygon*) geom)->getExteriorRing();
220 // assume TIN with with 4 points and point0 == point3
221 assert(cgeom->getNumPoints() == 4);
222 PositionVector corners;
223 for (int j = 0; j < 3; j++) {
224 Position pos((double) cgeom->getX(j), (double) cgeom->getY(j), (double) cgeom->getZ(j));
225 corners.push_back(pos);
226 myBoundary.add(pos);
227 }
228 addTriangle(corners);
229 numFeatures++;
230 } else {
231 WRITE_WARNINGF(TL("Ignored heightmap feature type %"), geom->getGeometryName());
232 }
233
234 /*
235 switch (gtype) {
236 case wkbPolygon: {
237 break;
238 }
239 case wkbPoint: {
240 WRITE_WARNING(TL("got wkbPoint"));
241 break;
242 }
243 case wkbLineString: {
244 WRITE_WARNING(TL("got wkbLineString"));
245 break;
246 }
247 case wkbMultiPoint: {
248 WRITE_WARNING(TL("got wkbMultiPoint"));
249 break;
250 }
251 case wkbMultiLineString: {
252 WRITE_WARNING(TL("got wkbMultiLineString"));
253 break;
254 }
255 case wkbMultiPolygon: {
256 WRITE_WARNING(TL("got wkbMultiPolygon"));
257 break;
258 }
259 default:
260 WRITE_WARNING(TL("Unsupported shape type occurred"));
261 break;
262 }
263 */
264 OGRFeature::DestroyFeature(feature);
265 }
266#if GDAL_VERSION_MAJOR < 2
267 OGRDataSource::DestroyDataSource(ds);
268#else
269 GDALClose(ds);
270#endif
271 OCTDestroyCoordinateTransformation(reinterpret_cast<OGRCoordinateTransformationH>(toWGS84));
272 OGRCleanupAll();
273 return numFeatures;
274#else
275 UNUSED_PARAMETER(file);
276 WRITE_ERROR(TL("Cannot load shape file since SUMO was compiled without GDAL support."));
277 return 0;
278#endif
279}
280
281
282int
283NBHeightMapper::loadTiff(const std::string& file) {
284#ifdef HAVE_GDAL
285 GDALAllRegister();
286 GDALDataset* poDataset = (GDALDataset*)GDALOpen(file.c_str(), GA_ReadOnly);
287 if (poDataset == 0) {
288 WRITE_ERROR(TL("Cannot load GeoTIFF file."));
289 return 0;
290 }
291 Boundary boundary;
292 const int xSize = poDataset->GetRasterXSize();
293 const int ySize = poDataset->GetRasterYSize();
294 double adfGeoTransform[6];
295 if (poDataset->GetGeoTransform(adfGeoTransform) == CE_None) {
296 Position topLeft(adfGeoTransform[0], adfGeoTransform[3]);
297 mySizeOfPixel.set(adfGeoTransform[1], adfGeoTransform[5]);
298 const double horizontalSize = xSize * mySizeOfPixel.x();
299 const double verticalSize = ySize * mySizeOfPixel.y();
300 boundary.add(topLeft);
301 boundary.add(topLeft.x() + horizontalSize, topLeft.y() + verticalSize);
302 } else {
303 WRITE_ERRORF(TL("Could not parse geo information from %."), file);
304 return 0;
305 }
306 const int picSize = xSize * ySize;
307 int16_t* raster = (int16_t*)CPLMalloc(sizeof(int16_t) * picSize);
308 bool ok = true;
309 for (int i = 1; i <= poDataset->GetRasterCount(); i++) {
310 GDALRasterBand* poBand = poDataset->GetRasterBand(i);
311 if (poBand->GetColorInterpretation() != GCI_GrayIndex) {
312 WRITE_ERRORF(TL("Unknown color band in %."), file);
313 clearData();
314 ok = false;
315 break;
316 }
317 assert(xSize == poBand->GetXSize() && ySize == poBand->GetYSize());
318 if (poBand->RasterIO(GF_Read, 0, 0, xSize, ySize, raster, xSize, ySize, GDT_Int16, 0, 0) == CE_Failure) {
319 WRITE_ERRORF(TL("Failure in reading %."), file);
320 clearData();
321 ok = false;
322 break;
323 }
324 }
325 double min = std::numeric_limits<double>::max();
326 double max = -std::numeric_limits<double>::max();
327 for (int i = 0; i < picSize; i++) {
328 min = MIN2(min, (double)raster[i]);
329 max = MAX2(max, (double)raster[i]);
330 }
331 GDALClose(poDataset);
332 if (ok) {
333 WRITE_MESSAGE("Read geotiff heightmap with size " + toString(xSize) + "," + toString(ySize)
334 + " for geo boundary [" + toString(boundary)
335 + "] with elevation range [" + toString(min) + "," + toString(max) + "].");
336 RasterData rasterData;
337 rasterData.raster = raster;
338 rasterData.boundary = boundary;
339 rasterData.xSize = xSize;
340 rasterData.ySize = ySize;
341 myRasters.push_back(rasterData);
342 return picSize;
343 } else {
344 return 0;
345 }
346#else
347 UNUSED_PARAMETER(file);
348 WRITE_ERROR(TL("Cannot load GeoTIFF file since SUMO was compiled without GDAL support."));
349 return 0;
350#endif
351}
352
353
354void
356 for (Triangles::iterator it = myTriangles.begin(); it != myTriangles.end(); it++) {
357 delete *it;
358 }
359 myTriangles.clear();
360#ifdef HAVE_GDAL
361 for (auto& item : myRasters) {
362 CPLFree(item.raster);
363 }
364 myRasters.clear();
365#endif
367}
368
369
370// ===========================================================================
371// Triangle member methods
372// ===========================================================================
374 myCorners(corners) {
375 assert(myCorners.size() == 3);
376 // @todo assert non-colinearity
377}
378
379
380void
382 queryResult.triangles.push_back(this);
383}
384
385
386bool
388 return myCorners.around(pos);
389}
390
391
392double
394 // en.wikipedia.org/wiki/Line-plane_intersection
395 Position p0 = myCorners.front();
396 Position line(0, 0, 1);
397 p0.sub(geo); // p0 - l0
398 Position normal = normalVector();
399 return p0.dotProduct(normal) / line.dotProduct(normal);
400}
401
402
405 // @todo maybe cache result to avoid multiple computations?
406 Position side1 = myCorners[1] - myCorners[0];
407 Position side2 = myCorners[2] - myCorners[0];
408 return side1.crossProduct(side2);
409}
410
411
412/****************************************************************************/
#define WRITE_WARNINGF(...)
Definition MsgHandler.h:296
#define WRITE_ERRORF(...)
Definition MsgHandler.h:305
#define WRITE_MESSAGE(msg)
Definition MsgHandler.h:297
#define WRITE_ERROR(msg)
Definition MsgHandler.h:304
#define WRITE_WARNING(msg)
Definition MsgHandler.h:295
#define TL(string)
Definition MsgHandler.h:315
#define TLF(string,...)
Definition MsgHandler.h:317
#define PROGRESS_BEGIN_MESSAGE(msg)
Definition MsgHandler.h:299
#define UNUSED_PARAMETER(x)
Definition StdDefs.h:30
T MIN2(T a, T b)
Definition StdDefs.h:76
T MAX2(T a, T b)
Definition StdDefs.h:82
std::string toString(const T &t, std::streamsize accuracy=gPrecision)
Definition ToString.h:46
A class that stores a 2D geometrical boundary.
Definition Boundary.h:39
void add(double x, double y, double z=0)
Makes the boundary include the given coordinate.
Definition Boundary.cpp:78
double ymin() const
Returns minimum y-coordinate.
Definition Boundary.cpp:130
void reset()
Resets the boundary.
Definition Boundary.cpp:66
double xmin() const
Returns minimum x-coordinate.
Definition Boundary.cpp:118
bool around(const Position &p, double offset=0) const
Returns whether the boundary contains the given coordinate.
Definition Boundary.cpp:172
double ymax() const
Returns maximum y-coordinate.
Definition Boundary.cpp:136
double xmax() const
Returns maximum x-coordinate.
Definition Boundary.cpp:124
virtual void endProcessMsg(std::string msg)
Ends a process information.
static MsgHandler * getMessageInstance()
Returns the instance to add normal messages to.
class for cirumventing the const-restriction of RTree::Search-context
Position normalVector() const
returns the normal vector for this triangles plane
double getZ(const Position &geo) const
returns the projection of the give geoCoordinate (WGS84) onto triangle plane
PositionVector myCorners
the corners of the triangle
Triangle(const PositionVector &corners)
void addSelf(const QueryResult &queryResult) const
callback for RTree search
bool contains(const Position &pos) const
checks whether pos lies within triangle (only checks x,y)
Set z-values for all network positions based on data from a height map.
std::vector< const Triangle * > Triangles
Triangles myTriangles
double getZ(const Position &geo) const
returns height for the given geo coordinate (WGS84)
static const NBHeightMapper & get()
return the singleton instance (maybe 0)
int loadShapeFile(const std::string &file)
load height data from Arcgis-shape file and returns the number of parsed features
static NBHeightMapper myInstance
the singleton instance
bool ready() const
returns whether the NBHeightMapper has data
std::vector< RasterData > myRasters
raster height information in m for all loaded files
NBHeightMapper()
private constructor and destructor (Singleton)
const Boundary & getBoundary()
returns the convex boundary of all known triangles
static void loadIfSet(OptionsCont &oc)
loads height map data if any loading options are set
void clearData()
clears loaded data
Position mySizeOfPixel
dimensions of one pixel in raster data
void addTriangle(PositionVector corners)
adds one triangles worth of height data
Boundary myBoundary
convex boundary of all known triangles;
TRIANGLE_RTREE_QUAL myRTree
The RTree for spatial queries.
int loadTiff(const std::string &file)
load height data from GeoTIFF file and returns the number of non void pixels
A storage for options typed value containers)
Definition OptionsCont.h:89
bool isSet(const std::string &name, bool failOnNonExistant=true) const
Returns the information whether the named option is set.
const StringVector & getStringVector(const std::string &name) const
Returns the list of string-value of the named option (only for Option_StringVector)
A point in 2D or 3D with translation and scaling methods.
Definition Position.h:37
double dotProduct(const Position &pos) const
returns the dot product (scalar product) between this point and the second one
Definition Position.h:304
void set(double x, double y)
set positions x and y
Definition Position.h:85
void sub(double dx, double dy)
Subtracts the given position from this one.
Definition Position.h:152
double x() const
Returns the x-position.
Definition Position.h:55
Position crossProduct(const Position &pos)
returns the cross product between this point and the second one
Definition Position.h:296
double y() const
Returns the y-position.
Definition Position.h:60
A list of positions.
Boundary getBoxBoundary() const
Returns a boundary enclosing this list of lines.