Line data Source code
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 : /****************************************************************************/
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 <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 : // ===========================================================================
57 : NBHeightMapper NBHeightMapper::myInstance;
58 :
59 :
60 : // ===========================================================================
61 : // method definitions
62 : // ===========================================================================
63 2040 : NBHeightMapper::NBHeightMapper():
64 2040 : myRTree(&Triangle::addSelf) {
65 2040 : }
66 :
67 :
68 2040 : NBHeightMapper::~NBHeightMapper() {
69 2040 : clearData();
70 2040 : }
71 :
72 :
73 : const NBHeightMapper&
74 428910 : NBHeightMapper::get() {
75 428910 : return myInstance;
76 : }
77 :
78 :
79 : bool
80 429517 : NBHeightMapper::ready() const {
81 429517 : return myRasters.size() > 0 || myTriangles.size() > 0;
82 : }
83 :
84 :
85 : double
86 607 : NBHeightMapper::getZ(const Position& geo) const {
87 607 : if (!ready()) {
88 0 : WRITE_WARNING(TL("Cannot supply height since no height data was loaded"));
89 0 : return 0;
90 : }
91 857 : for (auto& item : myRasters) {
92 604 : const Boundary& boundary = item.boundary;
93 604 : int16_t* raster = item.raster;
94 : double result = -1e6;
95 604 : if (boundary.around(geo)) {
96 354 : const int xSize = item.xSize;
97 354 : const double normX = (geo.x() - boundary.xmin()) / mySizeOfPixel.x();
98 354 : const double normY = (geo.y() - boundary.ymax()) / mySizeOfPixel.y();
99 354 : PositionVector corners;
100 354 : corners.push_back(Position(floor(normX) + 0.5, floor(normY) + 0.5, raster[(int)normY * xSize + (int)normX]));
101 354 : if (normX - floor(normX) > 0.5) {
102 180 : corners.push_back(Position(floor(normX) + 1.5, floor(normY) + 0.5, raster[(int)normY * xSize + (int)normX + 1]));
103 : } else {
104 174 : corners.push_back(Position(floor(normX) - 0.5, floor(normY) + 0.5, raster[(int)normY * xSize + (int)normX - 1]));
105 : }
106 354 : if (normY - floor(normY) > 0.5 && ((int)normY + 1) < item.ySize) {
107 169 : corners.push_back(Position(floor(normX) + 0.5, floor(normY) + 1.5, raster[((int)normY + 1) * xSize + (int)normX]));
108 : } else {
109 185 : corners.push_back(Position(floor(normX) + 0.5, floor(normY) - 0.5, raster[((int)normY - 1) * xSize + (int)normX]));
110 : }
111 354 : result = Triangle(corners).getZ(Position(normX, normY));
112 354 : }
113 354 : 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 253 : minB[0] = (float)geo.x() - 0.00001f;
121 253 : minB[1] = (float)geo.y() - 0.00001f;
122 253 : maxB[0] = (float)geo.x() + 0.00001f;
123 253 : maxB[1] = (float)geo.y() + 0.00001f;
124 : QueryResult queryResult;
125 253 : int hits = myRTree.Search(minB, maxB, queryResult);
126 253 : Triangles result = queryResult.triangles;
127 : assert(hits == (int)result.size());
128 : UNUSED_PARAMETER(hits); // only used for assertion
129 :
130 254 : for (Triangles::iterator it = result.begin(); it != result.end(); it++) {
131 4 : const Triangle* triangle = *it;
132 4 : if (triangle->contains(geo)) {
133 3 : return triangle->getZ(geo);
134 : }
135 : }
136 500 : WRITE_WARNINGF(TL("Could not get height data for coordinate %"), toString(geo));
137 250 : return 0;
138 253 : }
139 :
140 :
141 : void
142 3 : NBHeightMapper::addTriangle(PositionVector corners) {
143 3 : Triangle* triangle = new Triangle(corners);
144 3 : myTriangles.push_back(triangle);
145 3 : Boundary b = corners.getBoxBoundary();
146 3 : const float cmin[2] = {(float) b.xmin(), (float) b.ymin()};
147 3 : const float cmax[2] = {(float) b.xmax(), (float) b.ymax()};
148 3 : myRTree.Insert(cmin, cmax, triangle);
149 3 : }
150 :
151 :
152 : void
153 1892 : NBHeightMapper::loadIfSet(OptionsCont& oc) {
154 3784 : if (oc.isSet("heightmap.geotiff")) {
155 : // parse file(s)
156 6 : std::vector<std::string> files = oc.getStringVector("heightmap.geotiff");
157 6 : for (std::vector<std::string>::const_iterator file = files.begin(); file != files.end(); ++file) {
158 9 : PROGRESS_BEGIN_MESSAGE("Parsing from GeoTIFF '" + *file + "'");
159 3 : int numFeatures = myInstance.loadTiff(*file);
160 6 : MsgHandler::getMessageInstance()->endProcessMsg(
161 6 : " done (parsed " + toString(numFeatures) +
162 9 : " features, Boundary: " + toString(myInstance.getBoundary()) + ").");
163 : }
164 3 : }
165 3784 : if (oc.isSet("heightmap.shapefiles")) {
166 : // parse file(s)
167 0 : std::vector<std::string> files = oc.getStringVector("heightmap.shapefiles");
168 0 : for (std::vector<std::string>::const_iterator file = files.begin(); file != files.end(); ++file) {
169 0 : PROGRESS_BEGIN_MESSAGE("Parsing from shape-file '" + *file + "'");
170 0 : int numFeatures = myInstance.loadShapeFile(*file);
171 0 : MsgHandler::getMessageInstance()->endProcessMsg(
172 0 : " done (parsed " + toString(numFeatures) +
173 0 : " features, Boundary: " + toString(myInstance.getBoundary()) + ").");
174 : }
175 0 : }
176 1892 : }
177 :
178 :
179 : int
180 0 : NBHeightMapper::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 0 : GDALAllRegister();
187 0 : GDALDataset* ds = (GDALDataset*)GDALOpenEx(file.c_str(), GDAL_OF_VECTOR | GA_ReadOnly, nullptr, nullptr, nullptr);
188 : #endif
189 0 : if (ds == nullptr) {
190 0 : throw ProcessError(TLF("Could not open shape file '%'.", file));
191 : }
192 :
193 : // begin file parsing
194 0 : OGRLayer* layer = ds->GetLayer(0);
195 0 : layer->ResetReading();
196 :
197 : // triangle coordinates are stored in WGS84 and later matched with network coordinates in WGS84
198 : // build coordinate transformation
199 0 : OGRSpatialReference* sr_src = layer->GetSpatialRef();
200 0 : OGRSpatialReference sr_dest;
201 0 : sr_dest.SetWellKnownGeogCS("WGS84");
202 0 : OGRCoordinateTransformation* toWGS84 = OGRCreateCoordinateTransformation(sr_src, &sr_dest);
203 0 : if (toWGS84 == nullptr) {
204 0 : WRITE_WARNING(TL("Could not create geocoordinates converter; check whether proj.4 is installed."));
205 : }
206 :
207 : int numFeatures = 0;
208 : OGRFeature* feature;
209 0 : layer->ResetReading();
210 0 : while ((feature = layer->GetNextFeature()) != nullptr) {
211 0 : OGRGeometry* geom = feature->GetGeometryRef();
212 : assert(geom != 0);
213 :
214 0 : OGRwkbGeometryType gtype = geom->getGeometryType();
215 0 : if (gtype == wkbPolygon) {
216 : assert(std::string(geom->getGeometryName()) == std::string("POLYGON"));
217 : // try transform to wgs84
218 0 : geom->transform(toWGS84);
219 0 : OGRLinearRing* cgeom = ((OGRPolygon*) geom)->getExteriorRing();
220 : // assume TIN with with 4 points and point0 == point3
221 : assert(cgeom->getNumPoints() == 4);
222 0 : PositionVector corners;
223 0 : for (int j = 0; j < 3; j++) {
224 0 : Position pos((double) cgeom->getX(j), (double) cgeom->getY(j), (double) cgeom->getZ(j));
225 0 : corners.push_back(pos);
226 0 : myBoundary.add(pos);
227 : }
228 0 : addTriangle(corners);
229 0 : numFeatures++;
230 0 : } else {
231 0 : 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 0 : OGRFeature::DestroyFeature(feature);
265 : }
266 : #if GDAL_VERSION_MAJOR < 2
267 : OGRDataSource::DestroyDataSource(ds);
268 : #else
269 0 : GDALClose(ds);
270 : #endif
271 0 : OCTDestroyCoordinateTransformation(reinterpret_cast<OGRCoordinateTransformationH>(toWGS84));
272 0 : OGRCleanupAll();
273 0 : 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 0 : }
280 :
281 :
282 : int
283 3 : NBHeightMapper::loadTiff(const std::string& file) {
284 : #ifdef HAVE_GDAL
285 3 : GDALAllRegister();
286 3 : GDALDataset* poDataset = (GDALDataset*)GDALOpen(file.c_str(), GA_ReadOnly);
287 3 : if (poDataset == 0) {
288 0 : WRITE_ERROR(TL("Cannot load GeoTIFF file."));
289 0 : return 0;
290 : }
291 3 : Boundary boundary;
292 3 : const int xSize = poDataset->GetRasterXSize();
293 3 : const int ySize = poDataset->GetRasterYSize();
294 : double adfGeoTransform[6];
295 3 : if (poDataset->GetGeoTransform(adfGeoTransform) == CE_None) {
296 3 : Position topLeft(adfGeoTransform[0], adfGeoTransform[3]);
297 3 : mySizeOfPixel.set(adfGeoTransform[1], adfGeoTransform[5]);
298 3 : const double horizontalSize = xSize * mySizeOfPixel.x();
299 3 : const double verticalSize = ySize * mySizeOfPixel.y();
300 3 : boundary.add(topLeft);
301 3 : boundary.add(topLeft.x() + horizontalSize, topLeft.y() + verticalSize);
302 : } else {
303 0 : WRITE_ERRORF(TL("Could not parse geo information from %."), file);
304 0 : return 0;
305 : }
306 3 : const int picSize = xSize * ySize;
307 3 : int16_t* raster = (int16_t*)CPLMalloc(sizeof(int16_t) * picSize);
308 : bool ok = true;
309 6 : for (int i = 1; i <= poDataset->GetRasterCount(); i++) {
310 3 : GDALRasterBand* poBand = poDataset->GetRasterBand(i);
311 3 : if (poBand->GetColorInterpretation() != GCI_GrayIndex) {
312 0 : WRITE_ERRORF(TL("Unknown color band in %."), file);
313 0 : clearData();
314 : ok = false;
315 : break;
316 : }
317 : assert(xSize == poBand->GetXSize() && ySize == poBand->GetYSize());
318 3 : if (poBand->RasterIO(GF_Read, 0, 0, xSize, ySize, raster, xSize, ySize, GDT_Int16, 0, 0) == CE_Failure) {
319 0 : WRITE_ERRORF(TL("Failure in reading %."), file);
320 0 : clearData();
321 : ok = false;
322 : break;
323 : }
324 : }
325 3 : double min = std::numeric_limits<double>::max();
326 3 : double max = -std::numeric_limits<double>::max();
327 45363 : for (int i = 0; i < picSize; i++) {
328 45360 : min = MIN2(min, (double)raster[i]);
329 45510 : max = MAX2(max, (double)raster[i]);
330 : }
331 3 : GDALClose(poDataset);
332 3 : if (ok) {
333 21 : 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 3 : rasterData.raster = raster;
338 : rasterData.boundary = boundary;
339 3 : rasterData.xSize = xSize;
340 3 : rasterData.ySize = ySize;
341 3 : 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 3 : }
352 :
353 :
354 : void
355 2041 : NBHeightMapper::clearData() {
356 2044 : for (Triangles::iterator it = myTriangles.begin(); it != myTriangles.end(); it++) {
357 6 : delete *it;
358 : }
359 : myTriangles.clear();
360 : #ifdef HAVE_GDAL
361 2044 : for (auto& item : myRasters) {
362 3 : CPLFree(item.raster);
363 : }
364 : myRasters.clear();
365 : #endif
366 2041 : myBoundary.reset();
367 2041 : }
368 :
369 :
370 : // ===========================================================================
371 : // Triangle member methods
372 : // ===========================================================================
373 357 : NBHeightMapper::Triangle::Triangle(const PositionVector& corners):
374 : myCorners(corners) {
375 : assert(myCorners.size() == 3);
376 : // @todo assert non-colinearity
377 357 : }
378 :
379 :
380 : void
381 5 : NBHeightMapper::Triangle::addSelf(const QueryResult& queryResult) const {
382 5 : queryResult.triangles.push_back(this);
383 5 : }
384 :
385 :
386 : bool
387 4 : NBHeightMapper::Triangle::contains(const Position& pos) const {
388 4 : return myCorners.around(pos);
389 : }
390 :
391 :
392 : double
393 357 : NBHeightMapper::Triangle::getZ(const Position& geo) const {
394 : // en.wikipedia.org/wiki/Line-plane_intersection
395 357 : Position p0 = myCorners.front();
396 : Position line(0, 0, 1);
397 : p0.sub(geo); // p0 - l0
398 357 : Position normal = normalVector();
399 357 : return p0.dotProduct(normal) / line.dotProduct(normal);
400 : }
401 :
402 :
403 : Position
404 357 : NBHeightMapper::Triangle::normalVector() const {
405 : // @todo maybe cache result to avoid multiple computations?
406 357 : Position side1 = myCorners[1] - myCorners[0];
407 357 : Position side2 = myCorners[2] - myCorners[0];
408 357 : return side1.crossProduct(side2);
409 : }
410 :
411 :
412 : /****************************************************************************/
|
