36 #pragma GCC diagnostic push
37 #pragma GCC diagnostic ignored "-Wpedantic"
39 #include <ogrsf_frmts.h>
41 #include <gdal_priv.h>
43 #pragma GCC diagnostic pop
81 WRITE_WARNING(
"Cannot supply height since no height data was loaded");
85 const Boundary& boundary = item.boundary;
86 int16_t* raster = item.raster;
88 if (boundary.
around(geo)) {
89 const int xSize = item.xSize;
93 corners.push_back(
Position(floor(normX) + 0.5, floor(normY) + 0.5, raster[(int)normY * xSize + (
int)normX]));
94 if (normX - floor(normX) > 0.5) {
95 corners.push_back(
Position(floor(normX) + 1.5, floor(normY) + 0.5, raster[(int)normY * xSize + (
int)normX + 1]));
97 corners.push_back(
Position(floor(normX) - 0.5, floor(normY) + 0.5, raster[(int)normY * xSize + (
int)normX - 1]));
99 if (normY - floor(normY) > 0.5 && ((int)normY + 1) < item.ySize) {
100 corners.push_back(
Position(floor(normX) + 0.5, floor(normY) + 1.5, raster[((int)normY + 1) * xSize + (int)normX]));
102 corners.push_back(
Position(floor(normX) + 0.5, floor(normY) - 0.5, raster[((int)normY - 1) * xSize + (int)normX]));
106 if (result > -1e5 && result < 1e5) {
113 minB[0] = (float)geo.
x() - 0.00001f;
114 minB[1] = (float)geo.
y() - 0.00001f;
115 maxB[0] = (float)geo.
x() + 0.00001f;
116 maxB[1] = (float)geo.
y() + 0.00001f;
118 int hits =
myRTree.Search(minB, maxB, queryResult);
120 assert(hits == (
int)result.size());
123 for (Triangles::iterator it = result.begin(); it != result.end(); it++) {
126 return triangle->
getZ(geo);
139 const float cmin[2] = {(float) b.
xmin(), (float) b.
ymin()};
140 const float cmax[2] = {(float) b.
xmax(), (float) b.
ymax()};
141 myRTree.Insert(cmin, cmax, triangle);
147 if (oc.
isSet(
"heightmap.geotiff")) {
149 std::vector<std::string> files = oc.
getStringVector(
"heightmap.geotiff");
150 for (std::vector<std::string>::const_iterator file = files.begin(); file != files.end(); ++file) {
154 " done (parsed " +
toString(numFeatures) +
158 if (oc.
isSet(
"heightmap.shapefiles")) {
160 std::vector<std::string> files = oc.
getStringVector(
"heightmap.shapefiles");
161 for (std::vector<std::string>::const_iterator file = files.begin(); file != files.end(); ++file) {
165 " done (parsed " +
toString(numFeatures) +
175 #if GDAL_VERSION_MAJOR < 2
177 OGRDataSource* ds = OGRSFDriverRegistrar::Open(file.c_str(), FALSE);
180 GDALDataset* ds = (GDALDataset*)GDALOpenEx(file.c_str(), GDAL_OF_VECTOR | GA_ReadOnly,
nullptr,
nullptr,
nullptr);
183 throw ProcessError(
"Could not open shape file '" + file +
"'.");
187 OGRLayer* layer = ds->GetLayer(0);
188 layer->ResetReading();
192 OGRSpatialReference* sr_src = layer->GetSpatialRef();
193 OGRSpatialReference sr_dest;
194 sr_dest.SetWellKnownGeogCS(
"WGS84");
195 OGRCoordinateTransformation* toWGS84 = OGRCreateCoordinateTransformation(sr_src, &sr_dest);
196 if (toWGS84 ==
nullptr) {
197 WRITE_WARNING(
"Could not create geocoordinates converter; check whether proj.4 is installed.");
202 layer->ResetReading();
203 while ((feature = layer->GetNextFeature()) !=
nullptr) {
204 OGRGeometry* geom = feature->GetGeometryRef();
207 OGRwkbGeometryType gtype = geom->getGeometryType();
208 if (gtype == wkbPolygon) {
209 assert(std::string(geom->getGeometryName()) == std::string(
"POLYGON"));
211 geom->transform(toWGS84);
212 OGRLinearRing* cgeom = ((OGRPolygon*) geom)->getExteriorRing();
214 assert(cgeom->getNumPoints() == 4);
216 for (
int j = 0; j < 3; j++) {
217 Position pos((
double) cgeom->getX(j), (
double) cgeom->getY(j), (
double) cgeom->getZ(j));
218 corners.push_back(pos);
224 WRITE_WARNINGF(
"Ignored heightmap feature type %", geom->getGeometryName());
257 OGRFeature::DestroyFeature(feature);
259 #if GDAL_VERSION_MAJOR < 2
260 OGRDataSource::DestroyDataSource(ds);
264 OCTDestroyCoordinateTransformation(
reinterpret_cast<OGRCoordinateTransformationH
>(toWGS84));
269 WRITE_ERROR(
"Cannot load shape file since SUMO was compiled without GDAL support.");
279 GDALDataset* poDataset = (GDALDataset*)GDALOpen(file.c_str(), GA_ReadOnly);
280 if (poDataset == 0) {
285 const int xSize = poDataset->GetRasterXSize();
286 const int ySize = poDataset->GetRasterYSize();
287 double adfGeoTransform[6];
288 if (poDataset->GetGeoTransform(adfGeoTransform) == CE_None) {
289 Position topLeft(adfGeoTransform[0], adfGeoTransform[3]);
293 boundary.
add(topLeft);
294 boundary.
add(topLeft.
x() + horizontalSize, topLeft.
y() + verticalSize);
296 WRITE_ERROR(
"Could not parse geo information from " + file +
".");
299 const int picSize = xSize * ySize;
300 int16_t* raster = (int16_t*)CPLMalloc(
sizeof(int16_t) * picSize);
302 for (
int i = 1; i <= poDataset->GetRasterCount(); i++) {
303 GDALRasterBand* poBand = poDataset->GetRasterBand(i);
304 if (poBand->GetColorInterpretation() != GCI_GrayIndex) {
305 WRITE_ERROR(
"Unknown color band in " + file +
".");
310 assert(xSize == poBand->GetXSize() && ySize == poBand->GetYSize());
311 if (poBand->RasterIO(GF_Read, 0, 0, xSize, ySize, raster, xSize, ySize, GDT_Int16, 0, 0) == CE_Failure) {
318 double min = std::numeric_limits<double>::max();
319 double max = -std::numeric_limits<double>::max();
320 for (
int i = 0; i < picSize; i++) {
321 min =
MIN2(min, (
double)raster[i]);
322 max =
MAX2(max, (
double)raster[i]);
324 GDALClose(poDataset);
327 +
" for geo boundary [" +
toString(boundary)
330 rasterData.
raster = raster;
332 rasterData.
xSize = xSize;
333 rasterData.
ySize = ySize;
341 WRITE_ERROR(
"Cannot load GeoTIFF file since SUMO was compiled without GDAL support.");
355 CPLFree(item.raster);
381 return myCorners.around(pos);
399 Position side1 = myCorners[1] - myCorners[0];
400 Position side2 = myCorners[2] - myCorners[0];
#define WRITE_WARNINGF(...)
#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 AbstractPoly 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)
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.
const Boundary & getBoundary()
returns the convex boundary of all known triangles
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.
void set(double x, double y)
set positions x and y
double dotProduct(const Position &pos)
returns the dot product (scalar product) between this point and the second one
void sub(double dx, double dy)
Substracts 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.