37#pragma warning(disable: 4435 5219 5220)
40#pragma GCC diagnostic push
41#pragma GCC diagnostic ignored "-Wpedantic"
43#include <ogrsf_frmts.h>
47#pragma GCC diagnostic pop
88 WRITE_WARNING(
TL(
"Cannot supply height since no height data was loaded"));
92 const Boundary& boundary = item.boundary;
93 int16_t* raster = item.raster;
95 if (boundary.
around(geo)) {
96 const int xSize = item.xSize;
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]));
104 corners.push_back(
Position(floor(normX) - 0.5, floor(normY) + 0.5, raster[(int)normY * xSize + (
int)normX - 1]));
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]));
109 corners.push_back(
Position(floor(normX) + 0.5, floor(normY) - 0.5, raster[((int)normY - 1) * xSize + (int)normX]));
113 if (result > -1e5 && result < 1e5) {
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;
125 int hits =
myRTree.Search(minB, maxB, queryResult);
127 assert(hits == (
int)result.size());
130 for (Triangles::iterator it = result.begin(); it != result.end(); it++) {
133 return triangle->
getZ(geo);
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);
154 if (oc.
isSet(
"heightmap.geotiff")) {
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) {
161 " done (parsed " +
toString(numFeatures) +
165 if (oc.
isSet(
"heightmap.shapefiles")) {
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) {
172 " done (parsed " +
toString(numFeatures) +
182#if GDAL_VERSION_MAJOR < 2
184 OGRDataSource* ds = OGRSFDriverRegistrar::Open(file.c_str(), FALSE);
187 GDALDataset* ds = (GDALDataset*)GDALOpenEx(file.c_str(), GDAL_OF_VECTOR | GA_ReadOnly,
nullptr,
nullptr,
nullptr);
194 OGRLayer* layer = ds->GetLayer(0);
195 layer->ResetReading();
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."));
209 layer->ResetReading();
210 while ((feature = layer->GetNextFeature()) !=
nullptr) {
211 OGRGeometry* geom = feature->GetGeometryRef();
214 OGRwkbGeometryType gtype = geom->getGeometryType();
215 if (gtype == wkbPolygon) {
216 assert(std::string(geom->getGeometryName()) == std::string(
"POLYGON"));
218 geom->transform(toWGS84);
219 OGRLinearRing* cgeom = ((OGRPolygon*) geom)->getExteriorRing();
221 assert(cgeom->getNumPoints() == 4);
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);
231 WRITE_WARNINGF(
TL(
"Ignored heightmap feature type %"), geom->getGeometryName());
264 OGRFeature::DestroyFeature(feature);
266#if GDAL_VERSION_MAJOR < 2
267 OGRDataSource::DestroyDataSource(ds);
271 OCTDestroyCoordinateTransformation(
reinterpret_cast<OGRCoordinateTransformationH
>(toWGS84));
276 WRITE_ERROR(
TL(
"Cannot load shape file since SUMO was compiled without GDAL support."));
286 GDALDataset* poDataset = (GDALDataset*)GDALOpen(file.c_str(), GA_ReadOnly);
287 if (poDataset == 0) {
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]);
300 boundary.
add(topLeft);
301 boundary.
add(topLeft.
x() + horizontalSize, topLeft.
y() + verticalSize);
303 WRITE_ERRORF(
TL(
"Could not parse geo information from %."), file);
306 const int picSize = xSize * ySize;
307 int16_t* raster = (int16_t*)CPLMalloc(
sizeof(int16_t) * picSize);
309 for (
int i = 1; i <= poDataset->GetRasterCount(); i++) {
310 GDALRasterBand* poBand = poDataset->GetRasterBand(i);
311 if (poBand->GetColorInterpretation() != GCI_GrayIndex) {
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) {
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]);
331 GDALClose(poDataset);
334 +
" for geo boundary [" +
toString(boundary)
337 rasterData.
raster = raster;
339 rasterData.
xSize = xSize;
340 rasterData.
ySize = ySize;
348 WRITE_ERROR(
TL(
"Cannot load GeoTIFF file since SUMO was compiled without GDAL support."));
362 CPLFree(item.raster);
388 return myCorners.around(pos);
406 Position side1 = myCorners[1] - myCorners[0];
407 Position side2 = myCorners[2] - myCorners[0];
#define WRITE_WARNINGF(...)
#define WRITE_ERRORF(...)
#define WRITE_MESSAGE(msg)
#define WRITE_WARNING(msg)
#define PROGRESS_BEGIN_MESSAGE(msg)
#define UNUSED_PARAMETER(x)
std::string toString(const T &t, std::streamsize accuracy=gPrecision)
A class that stores a 2D geometrical boundary.
void add(double x, double y, double z=0)
Makes the boundary include the given coordinate.
double ymin() const
Returns minimum y-coordinate.
void reset()
Resets the boundary.
double xmin() const
Returns minimum x-coordinate.
bool around(const Position &p, double offset=0) const
Returns whether the boundary contains the given coordinate.
double ymax() const
Returns maximum y-coordinate.
double xmax() const
Returns maximum x-coordinate.
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
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)
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.
double dotProduct(const Position &pos) const
returns the dot product (scalar product) between this point and the second one
void set(double x, double y)
set positions x and y
void sub(double dx, double dy)
Subtracts the given position from this one.
double x() const
Returns the x-position.
Position crossProduct(const Position &pos)
returns the cross product between this point and the second one
double y() const
Returns the y-position.
Boundary getBoxBoundary() const
Returns a boundary enclosing this list of lines.