Line data Source code
1 : /****************************************************************************/
2 : // Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.dev/sumo
3 : // Copyright (C) 2011-2026 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 : /****************************************************************************/
14 : /// @file NBHeightMapper.cpp
15 : /// @author Jakob Erdmann
16 : /// @author Laura Bieker
17 : /// @author Michael Behrisch
18 : /// @date Sept 2011
19 : ///
20 : // Set z-values for all network positions based on data from a height map
21 : /****************************************************************************/
22 : #include <config.h>
23 :
24 : #include <string>
25 : #include <utils/common/MsgHandler.h>
26 : #include <utils/common/ToString.h>
27 : #include <utils/common/StringUtils.h>
28 : #include <utils/options/OptionsCont.h>
29 : #include <utils/geom/GeomHelper.h>
30 : #include "NBHeightMapper.h"
31 : #include <utils/geom/GeoConvHelper.h>
32 : #include <utils/common/RGBColor.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 <gdal_version.h>
44 : #include <ogrsf_frmts.h>
45 : #include <ogr_api.h>
46 : #include <gdal_priv.h>
47 : #if __GNUC__ > 3
48 : #pragma GCC diagnostic pop
49 : #endif
50 : #ifdef _MSC_VER
51 : #pragma warning(pop)
52 : #endif
53 : #endif
54 :
55 : // ===========================================================================
56 : // static members
57 : // ===========================================================================
58 : NBHeightMapper NBHeightMapper::myInstance;
59 :
60 :
61 : // ===========================================================================
62 : // method definitions
63 : // ===========================================================================
64 2589 : NBHeightMapper::NBHeightMapper():
65 2589 : myRTree(&Triangle::addSelf) {
66 2589 : }
67 :
68 :
69 2589 : NBHeightMapper::~NBHeightMapper() {
70 2589 : clearData();
71 2589 : }
72 :
73 :
74 : const NBHeightMapper&
75 1010275 : NBHeightMapper::get() {
76 1010275 : return myInstance;
77 : }
78 :
79 :
80 : bool
81 1011810 : NBHeightMapper::ready() const {
82 1011810 : return myRasters.size() > 0 || myTriangles.size() > 0;
83 : }
84 :
85 :
86 : double
87 1535 : NBHeightMapper::getZ(const Position& geo) const {
88 1535 : if (!ready()) {
89 0 : WRITE_WARNING(TL("Cannot supply height since no height data was loaded"));
90 0 : return 0;
91 : }
92 2383 : for (auto& item : myRasters) {
93 1532 : const Boundary& boundary = item.boundary;
94 1532 : float* raster = item.raster;
95 : double result = -1e6;
96 :
97 1532 : double x = geo.x();
98 1532 : double y = geo.y();
99 :
100 : #ifdef HAVE_GDAL
101 : // Transform geo coordinates to the coordinate system of this
102 : // raster image for lookup in its raster, if applicable.
103 1532 : if (item.transform != nullptr) {
104 : // Since the input coordinates are always WGS84 (they may be
105 : // transformed to it in NBNetBuilder::transformCoordinate), and
106 : // WGS84 uses latitude-longitude order (y-x), we have to swap the
107 : // input coordinates here.
108 : std::swap(x, y);
109 :
110 1532 : item.transform->Transform(1, &x, &y);
111 : }
112 : #endif
113 :
114 1532 : if (boundary.around2D(x, y)) {
115 684 : const int xSize = item.xSize;
116 684 : const double normX = (x - boundary.xmin()) / mySizeOfPixel.x();
117 684 : const double normY = (y - boundary.ymax()) / mySizeOfPixel.y();
118 684 : PositionVector corners;
119 684 : corners.push_back(Position(floor(normX) + 0.5, floor(normY) + 0.5, raster[(int)normY * xSize + (int)normX]));
120 684 : if (normX - floor(normX) > 0.5) {
121 294 : corners.push_back(Position(floor(normX) + 1.5, floor(normY) + 0.5, raster[(int)normY * xSize + (int)normX + 1]));
122 : } else {
123 390 : corners.push_back(Position(floor(normX) - 0.5, floor(normY) + 0.5, raster[(int)normY * xSize + (int)normX - 1]));
124 : }
125 684 : if (normY - floor(normY) > 0.5 && ((int)normY + 1) < item.ySize) {
126 325 : corners.push_back(Position(floor(normX) + 0.5, floor(normY) + 1.5, raster[((int)normY + 1) * xSize + (int)normX]));
127 : } else {
128 359 : corners.push_back(Position(floor(normX) + 0.5, floor(normY) - 0.5, raster[((int)normY - 1) * xSize + (int)normX]));
129 : }
130 684 : result = Triangle(corners).getZ(Position(normX, normY));
131 684 : }
132 684 : if (result > -1e5 && result < 1e5) {
133 684 : return result;
134 : }
135 : }
136 : // coordinates in degrees hence a small search window
137 : float minB[2];
138 : float maxB[2];
139 851 : minB[0] = (float)geo.x() - 0.00001f;
140 851 : minB[1] = (float)geo.y() - 0.00001f;
141 851 : maxB[0] = (float)geo.x() + 0.00001f;
142 851 : maxB[1] = (float)geo.y() + 0.00001f;
143 : QueryResult queryResult;
144 851 : int hits = myRTree.Search(minB, maxB, queryResult);
145 851 : Triangles result = queryResult.triangles;
146 : assert(hits == (int)result.size());
147 : UNUSED_PARAMETER(hits); // only used for assertion
148 :
149 852 : for (Triangles::iterator it = result.begin(); it != result.end(); it++) {
150 4 : const Triangle* triangle = *it;
151 4 : if (triangle->contains(geo)) {
152 3 : return triangle->getZ(geo);
153 : }
154 : }
155 1696 : WRITE_WARNINGF(TL("Could not get height data for coordinate %"), toString(geo));
156 848 : return 0;
157 851 : }
158 :
159 :
160 : void
161 3 : NBHeightMapper::addTriangle(PositionVector corners) {
162 3 : Triangle* triangle = new Triangle(corners);
163 3 : myTriangles.push_back(triangle);
164 3 : Boundary b = corners.getBoxBoundary();
165 3 : const float cmin[2] = {(float) b.xmin(), (float) b.ymin()};
166 3 : const float cmax[2] = {(float) b.xmax(), (float) b.ymax()};
167 3 : myRTree.Insert(cmin, cmax, triangle);
168 3 : }
169 :
170 :
171 : void
172 2437 : NBHeightMapper::loadIfSet(OptionsCont& oc) {
173 4874 : if (oc.isSet("heightmap.geotiff")) {
174 : // parse file(s)
175 10 : std::vector<std::string> files = oc.getStringVector("heightmap.geotiff");
176 10 : for (std::vector<std::string>::const_iterator file = files.begin(); file != files.end(); ++file) {
177 15 : PROGRESS_BEGIN_MESSAGE("Parsing from GeoTIFF '" + *file + "'");
178 5 : int numFeatures = myInstance.loadTiff(*file);
179 10 : MsgHandler::getMessageInstance()->endProcessMsg(
180 10 : " done (parsed " + toString(numFeatures) +
181 15 : " features, Boundary: " + toString(myInstance.getBoundary()) + ").");
182 : }
183 5 : }
184 4874 : if (oc.isSet("heightmap.shapefiles")) {
185 : // parse file(s)
186 0 : std::vector<std::string> files = oc.getStringVector("heightmap.shapefiles");
187 0 : for (std::vector<std::string>::const_iterator file = files.begin(); file != files.end(); ++file) {
188 0 : PROGRESS_BEGIN_MESSAGE("Parsing from shape-file '" + *file + "'");
189 0 : int numFeatures = myInstance.loadShapeFile(*file);
190 0 : MsgHandler::getMessageInstance()->endProcessMsg(
191 0 : " done (parsed " + toString(numFeatures) +
192 0 : " features, Boundary: " + toString(myInstance.getBoundary()) + ").");
193 : }
194 0 : }
195 2437 : }
196 :
197 :
198 : int
199 0 : NBHeightMapper::loadShapeFile(const std::string& file) {
200 : #ifdef HAVE_GDAL
201 : #if GDAL_VERSION_MAJOR < 2
202 : OGRRegisterAll();
203 : OGRDataSource* ds = OGRSFDriverRegistrar::Open(file.c_str(), FALSE);
204 : #else
205 0 : GDALAllRegister();
206 0 : GDALDataset* ds = (GDALDataset*)GDALOpenEx(file.c_str(), GDAL_OF_VECTOR | GA_ReadOnly, nullptr, nullptr, nullptr);
207 : #endif
208 0 : if (ds == nullptr) {
209 0 : throw ProcessError(TLF("Could not open shape file '%'.", file));
210 : }
211 :
212 : // begin file parsing
213 0 : OGRLayer* layer = ds->GetLayer(0);
214 0 : layer->ResetReading();
215 :
216 : // triangle coordinates are stored in WGS84 and later matched with network coordinates in WGS84
217 : // build coordinate transformation
218 : #if GDAL_VERSION_MAJOR < 3
219 : OGRSpatialReference* sr_src = layer->GetSpatialRef();
220 : #else
221 0 : const OGRSpatialReference* sr_src = layer->GetSpatialRef();
222 : #endif
223 0 : OGRSpatialReference sr_dest;
224 0 : sr_dest.SetWellKnownGeogCS("WGS84");
225 0 : OGRCoordinateTransformation* toWGS84 = OGRCreateCoordinateTransformation(sr_src, &sr_dest);
226 0 : if (toWGS84 == nullptr) {
227 0 : WRITE_WARNING(TL("Could not create geocoordinates converter; check whether proj.4 is installed."));
228 : }
229 :
230 : int numFeatures = 0;
231 : OGRFeature* feature;
232 0 : layer->ResetReading();
233 0 : while ((feature = layer->GetNextFeature()) != nullptr) {
234 0 : OGRGeometry* geom = feature->GetGeometryRef();
235 : assert(geom != 0);
236 :
237 0 : OGRwkbGeometryType gtype = geom->getGeometryType();
238 0 : if (gtype == wkbPolygon) {
239 : assert(std::string(geom->getGeometryName()) == std::string("POLYGON"));
240 : // try transform to wgs84
241 0 : geom->transform(toWGS84);
242 0 : OGRLinearRing* cgeom = ((OGRPolygon*) geom)->getExteriorRing();
243 : // assume TIN with with 4 points and point0 == point3
244 : assert(cgeom->getNumPoints() == 4);
245 0 : PositionVector corners;
246 0 : for (int j = 0; j < 3; j++) {
247 0 : Position pos((double) cgeom->getX(j), (double) cgeom->getY(j), (double) cgeom->getZ(j));
248 0 : corners.push_back(pos);
249 0 : myBoundary.add(pos);
250 : }
251 0 : addTriangle(corners);
252 0 : numFeatures++;
253 0 : } else {
254 0 : WRITE_WARNINGF(TL("Ignored heightmap feature type %"), geom->getGeometryName());
255 : }
256 :
257 : /*
258 : switch (gtype) {
259 : case wkbPolygon: {
260 : break;
261 : }
262 : case wkbPoint: {
263 : WRITE_WARNING(TL("got wkbPoint"));
264 : break;
265 : }
266 : case wkbLineString: {
267 : WRITE_WARNING(TL("got wkbLineString"));
268 : break;
269 : }
270 : case wkbMultiPoint: {
271 : WRITE_WARNING(TL("got wkbMultiPoint"));
272 : break;
273 : }
274 : case wkbMultiLineString: {
275 : WRITE_WARNING(TL("got wkbMultiLineString"));
276 : break;
277 : }
278 : case wkbMultiPolygon: {
279 : WRITE_WARNING(TL("got wkbMultiPolygon"));
280 : break;
281 : }
282 : default:
283 : WRITE_WARNING(TL("Unsupported shape type occurred"));
284 : break;
285 : }
286 : */
287 0 : OGRFeature::DestroyFeature(feature);
288 : }
289 : #if GDAL_VERSION_MAJOR < 2
290 : OGRDataSource::DestroyDataSource(ds);
291 : #else
292 0 : GDALClose(ds);
293 : #endif
294 0 : OCTDestroyCoordinateTransformation(reinterpret_cast<OGRCoordinateTransformationH>(toWGS84));
295 0 : OGRCleanupAll();
296 0 : return numFeatures;
297 : #else
298 : UNUSED_PARAMETER(file);
299 : WRITE_ERROR(TL("Cannot load shape file since SUMO was compiled without GDAL support."));
300 : return 0;
301 : #endif
302 0 : }
303 :
304 :
305 : int
306 5 : NBHeightMapper::loadTiff(const std::string& file) {
307 : #ifdef HAVE_GDAL
308 5 : GDALAllRegister();
309 5 : GDALDataset* poDataset = (GDALDataset*)GDALOpen(file.c_str(), GA_ReadOnly);
310 5 : if (poDataset == 0) {
311 0 : WRITE_ERROR(TL("Cannot load GeoTIFF file."));
312 0 : return 0;
313 : }
314 5 : Boundary boundary;
315 5 : const int xSize = poDataset->GetRasterXSize();
316 5 : const int ySize = poDataset->GetRasterYSize();
317 : double adfGeoTransform[6];
318 5 : if (poDataset->GetGeoTransform(adfGeoTransform) == CE_None) {
319 5 : Position topLeft(adfGeoTransform[0], adfGeoTransform[3]);
320 5 : mySizeOfPixel.set(adfGeoTransform[1], adfGeoTransform[5]);
321 5 : const double horizontalSize = xSize * mySizeOfPixel.x();
322 5 : const double verticalSize = ySize * mySizeOfPixel.y();
323 5 : boundary.add(topLeft);
324 5 : boundary.add(topLeft.x() + horizontalSize, topLeft.y() + verticalSize);
325 : } else {
326 0 : WRITE_ERRORF(TL("Could not parse geo information from %."), file);
327 0 : return 0;
328 : }
329 5 : const int picSize = xSize * ySize;
330 5 : float* raster = (float*)CPLMalloc(sizeof(float) * picSize);
331 : bool ok = true;
332 10 : for (int i = 1; i <= poDataset->GetRasterCount(); i++) {
333 5 : GDALRasterBand* poBand = poDataset->GetRasterBand(i);
334 5 : if (poBand->GetColorInterpretation() != GCI_GrayIndex) {
335 0 : WRITE_ERRORF(TL("Unknown color band in %."), file);
336 0 : clearData();
337 : ok = false;
338 : break;
339 : }
340 : assert(xSize == poBand->GetXSize() && ySize == poBand->GetYSize());
341 5 : if (poBand->RasterIO(GF_Read, 0, 0, xSize, ySize, raster, xSize, ySize, GDT_Float32, 0, 0) == CE_Failure) {
342 0 : WRITE_ERRORF(TL("Failure in reading %."), file);
343 0 : clearData();
344 : ok = false;
345 : break;
346 : }
347 : }
348 5 : double min = std::numeric_limits<double>::max();
349 5 : double max = -std::numeric_limits<double>::max();
350 75605 : for (int i = 0; i < picSize; i++) {
351 75600 : min = MIN2(min, (double)raster[i]);
352 75850 : max = MAX2(max, (double)raster[i]);
353 : }
354 :
355 : // Make a copy, GDALClose will destroy the original
356 : #if GDAL_VERSION_MAJOR < 3
357 : OGRSpatialReference spatialRef;
358 : char* wkt = const_cast<char*>(poDataset->GetProjectionRef());
359 : spatialRef.importFromWkt(&wkt);
360 : #else
361 5 : OGRSpatialReference spatialRef(*poDataset->GetSpatialRef());
362 : #endif
363 5 : GDALClose(poDataset);
364 5 : if (ok) {
365 35 : WRITE_MESSAGE("Read geotiff heightmap with size " + toString(xSize) + "," + toString(ySize)
366 : + " for geo boundary [" + toString(boundary)
367 : + "] with elevation range [" + toString(min) + "," + toString(max) + "].");
368 5 : OGRSpatialReference wgs;
369 5 : wgs.SetWellKnownGeogCS("WGS84");
370 5 : myRasters.push_back(RasterData{raster, boundary, xSize, ySize, OGRCreateCoordinateTransformation(&wgs, &spatialRef)});
371 : return picSize;
372 5 : }
373 : return 0;
374 : #else
375 : UNUSED_PARAMETER(file);
376 : WRITE_ERROR(TL("Cannot load GeoTIFF file since SUMO was compiled without GDAL support."));
377 : return 0;
378 : #endif
379 5 : }
380 :
381 :
382 : void
383 2590 : NBHeightMapper::clearData() {
384 2593 : for (Triangles::iterator it = myTriangles.begin(); it != myTriangles.end(); it++) {
385 6 : delete *it;
386 : }
387 : myTriangles.clear();
388 : #ifdef HAVE_GDAL
389 2595 : for (auto& item : myRasters) {
390 5 : CPLFree(item.raster);
391 5 : if (item.transform != nullptr) {
392 5 : delete item.transform;
393 : }
394 : }
395 : myRasters.clear();
396 : #endif
397 2590 : myBoundary.reset();
398 2590 : }
399 :
400 :
401 : // ===========================================================================
402 : // Triangle member methods
403 : // ===========================================================================
404 687 : NBHeightMapper::Triangle::Triangle(const PositionVector& corners):
405 : myCorners(corners) {
406 : assert(myCorners.size() == 3);
407 : // @todo assert non-colinearity
408 687 : }
409 :
410 :
411 : void
412 5 : NBHeightMapper::Triangle::addSelf(const QueryResult& queryResult) const {
413 5 : queryResult.triangles.push_back(this);
414 5 : }
415 :
416 :
417 : bool
418 4 : NBHeightMapper::Triangle::contains(const Position& pos) const {
419 4 : return myCorners.around(pos);
420 : }
421 :
422 :
423 : double
424 687 : NBHeightMapper::Triangle::getZ(const Position& geo) const {
425 : // en.wikipedia.org/wiki/Line-plane_intersection
426 687 : Position p0 = myCorners.front();
427 : Position line(0, 0, 1);
428 : p0.sub(geo); // p0 - l0
429 687 : Position normal = normalVector();
430 687 : return p0.dotProduct(normal) / line.dotProduct(normal);
431 : }
432 :
433 :
434 : Position
435 687 : NBHeightMapper::Triangle::normalVector() const {
436 : // @todo maybe cache result to avoid multiple computations?
437 687 : Position side1 = myCorners[1] - myCorners[0];
438 687 : Position side2 = myCorners[2] - myCorners[0];
439 687 : return side1.crossProduct(side2);
440 : }
441 :
442 :
443 : /****************************************************************************/
|
