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.first;
86 int16_t* raster = item.second;
88 if (boundary.
around(geo)) {
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) {
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();
208 assert(std::string(geom->getGeometryName()) == std::string(
"POLYGON"));
210 geom->transform(toWGS84);
211 OGRLinearRing* cgeom = ((OGRPolygon*) geom)->getExteriorRing();
213 assert(cgeom->getNumPoints() == 4);
215 for (
int j = 0; j < 3; j++) {
216 Position pos((
double) cgeom->getX(j), (
double) cgeom->getY(j), (
double) cgeom->getZ(j));
217 corners.push_back(pos);
254 OGRFeature::DestroyFeature(feature);
256 #if GDAL_VERSION_MAJOR < 2
257 OGRDataSource::DestroyDataSource(ds);
261 OCTDestroyCoordinateTransformation(
reinterpret_cast<OGRCoordinateTransformationH
>(toWGS84));
266 WRITE_ERROR(
"Cannot load shape file since SUMO was compiled without GDAL support.");
276 GDALDataset* poDataset = (GDALDataset*)GDALOpen(file.c_str(), GA_ReadOnly);
277 if (poDataset == 0) {
282 const int xSize = poDataset->GetRasterXSize();
283 const int ySize = poDataset->GetRasterYSize();
284 double adfGeoTransform[6];
285 if (poDataset->GetGeoTransform(adfGeoTransform) == CE_None) {
286 Position topLeft(adfGeoTransform[0], adfGeoTransform[3]);
290 boundary.
add(topLeft);
291 boundary.
add(topLeft.
x() + horizontalSize, topLeft.
y() + verticalSize);
293 WRITE_ERROR(
"Could not parse geo information from " + file +
".");
296 const int picSize = xSize * ySize;
297 int16_t* raster = (int16_t*)CPLMalloc(
sizeof(int16_t) * picSize);
298 for (
int i = 1; i <= poDataset->GetRasterCount(); i++) {
299 GDALRasterBand* poBand = poDataset->GetRasterBand(i);
300 if (poBand->GetColorInterpretation() != GCI_GrayIndex) {
301 WRITE_ERROR(
"Unknown color band in " + file +
".");
305 if (poBand->GetRasterDataType() != GDT_Int16) {
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) {
317 GDALClose(poDataset);
318 myRasters.push_back(std::make_pair(boundary, raster));
322 WRITE_ERROR(
"Cannot load GeoTIFF file since SUMO was compiled without GDAL support.");
336 CPLFree(item.second);
362 return myCorners.around(pos);
380 Position side1 = myCorners[1] - myCorners[0];
381 Position side2 = myCorners[2] - myCorners[0];
#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< std::pair< Boundary, int16_t * > > 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.