44 #pragma GCC diagnostic push 45 #pragma GCC diagnostic ignored "-Wpedantic" 47 #include <ogrsf_frmts.h> 49 #include <gdal_priv.h> 51 #pragma GCC diagnostic pop 66 myRTree(&
Triangle::addSelf), myRaster(0) {
90 WRITE_WARNING(
"Cannot supply height since no height data was loaded");
100 corners.push_back(
Position(floor(normX) + 0.5, floor(normY) + 0.5,
myRaster[(
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,
myRaster[(
int)normY * xSize + (
int)normX + 1]));
104 corners.push_back(
Position(floor(normX) - 0.5, floor(normY) + 0.5,
myRaster[(
int)normY * xSize + (
int)normX - 1]));
106 if (normY - floor(normY) > 0.5) {
107 corners.push_back(
Position(floor(normX) + 0.5, floor(normY) + 1.5,
myRaster[((
int)normY + 1) * xSize + (
int)normX]));
109 corners.push_back(
Position(floor(normX) + 0.5, floor(normY) - 0.5,
myRaster[((
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, NULL, NULL, NULL);
190 throw ProcessError(
"Could not open shape file '" + file +
"'.");
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);
204 WRITE_WARNING(
"Could not create geocoordinates converter; check whether proj.4 is installed.");
209 layer->ResetReading();
210 while ((feature = layer->GetNextFeature()) != NULL) {
211 OGRGeometry* geom = feature->GetGeometryRef();
215 assert(std::string(geom->getGeometryName()) == std::string(
"POLYGON"));
217 geom->transform(toWGS84);
218 OGRLinearRing* cgeom = ((OGRPolygon*) geom)->getExteriorRing();
220 assert(cgeom->getNumPoints() == 4);
222 for (
int j = 0; j < 3; j++) {
223 Position pos((
double) cgeom->getX(j), (double) cgeom->getY(j), (double) cgeom->getZ(j));
224 corners.push_back(pos);
261 OGRFeature::DestroyFeature(feature);
263 #if GDAL_VERSION_MAJOR < 2 264 OGRDataSource::DestroyDataSource(ds);
268 OCTDestroyCoordinateTransformation(toWGS84);
272 WRITE_ERROR(
"Cannot load shape file since SUMO was compiled without GDAL support.");
282 GDALDataset* poDataset = (GDALDataset*)GDALOpen(file.c_str(), GA_ReadOnly);
283 if (poDataset == 0) {
287 const int xSize = poDataset->GetRasterXSize();
288 const int ySize = poDataset->GetRasterYSize();
289 double adfGeoTransform[6];
290 if (poDataset->GetGeoTransform(adfGeoTransform) == CE_None) {
291 Position topLeft(adfGeoTransform[0], adfGeoTransform[3]);
296 myBoundary.
add(topLeft.
x() + horizontalSize, topLeft.
y() + verticalSize);
298 WRITE_ERROR(
"Could not parse geo information from " + file +
".");
301 const int picSize = xSize * ySize;
302 myRaster = (int16_t*)CPLMalloc(
sizeof(int16_t) * picSize);
303 for (
int i = 1; i <= poDataset->GetRasterCount(); i++) {
304 GDALRasterBand* poBand = poDataset->GetRasterBand(i);
305 if (poBand->GetColorInterpretation() != GCI_GrayIndex) {
306 WRITE_ERROR(
"Unknown color band in " + file +
".");
310 if (poBand->GetRasterDataType() != GDT_Int16) {
315 assert(xSize == poBand->GetXSize() && ySize == poBand->GetYSize());
316 if (poBand->RasterIO(GF_Read, 0, 0, xSize, ySize,
myRaster, xSize, ySize, GDT_Int16, 0, 0) == CE_Failure) {
322 GDALClose(poDataset);
325 WRITE_ERROR(
"Cannot load GeoTIFF file since SUMO was compiled without GDAL support.");
384 Position side2 = myCorners[2] - myCorners[0];
bool around(const Position &p, double offset=0) const
Returns the information whether the position vector describes a polygon lying around the given point...
double ymin() const
Returns minimum y-coordinate.
double xmax() const
Returns maximum x-coordinate.
int loadShapeFile(const std::string &file)
load height data from Arcgis-shape file and returns the number of parsed features ...
Position normalVector() const
returns the normal vector for this triangles plane
double y() const
Returns the y-position.
double x() const
Returns the x-position.
Triangle(const PositionVector &corners)
void addTriangle(PositionVector corners)
adds one triangles worth of height data
void clearData()
clears loaded data
static NBHeightMapper Singleton
the singleton instance
void set(double x, double y)
set positions x and y
NBHeightMapper()
private constructor and destructor (Singleton)
#define UNUSED_PARAMETER(x)
A class that stores a 2D geometrical boundary.
#define WRITE_WARNING(msg)
int loadTiff(const std::string &file)
load height data from GeoTIFF file and returns the number of non void pixels
static const NBHeightMapper & get()
return the singleton instance (maybe 0)
bool isSet(const std::string &name, bool failOnNonExistant=true) const
Returns the information whether the named option is set.
class for cirumventing the const-restriction of RTree::Search-context
static void loadIfSet(OptionsCont &oc)
loads heigh map data if any loading options are set
TRIANGLE_RTREE_QUAL myRTree
The RTree for spatial queries.
Position crossProduct(const Position &pos)
returns the cross product between this point and the second one
void addSelf(const QueryResult &queryResult) const
callback for RTree search
std::string toString(const T &t, std::streamsize accuracy=gPrecision)
A point in 2D or 3D with translation and scaling methods.
bool contains(const Position &pos) const
checks whether pos lies within triangle (only checks x,y)
std::vector< const Triangle * > Triangles
std::vector< std::string > getStringVector(const std::string &name) const
Returns the list of string-vector-value of the named option (only for Option_String) ...
double xmin() const
Returns minimum x-coordinate.
#define PROGRESS_BEGIN_MESSAGE(msg)
static MsgHandler * getMessageInstance()
Returns the instance to add normal messages to.
int16_t * myRaster
raster height information in m
Position mySizeOfPixel
dimensions of one pixel in raster data
Boundary myBoundary
convex boundary of all known triangles;
void reset()
Resets the boundary.
double dotProduct(const Position &pos)
returns the dot product (scalar product) between this point and the second one
bool around(const Position &p, double offset=0) const
Returns whether the boundary contains the given coordinate.
PositionVector myCorners
the corners of the triangle
A storage for options typed value containers)
double getZ(const Position &geo) const
returns the projection of the give geoCoordinate (WGS84) onto triangle plane
Boundary getBoxBoundary() const
Returns a boundary enclosing this list of lines.
double getZ(const Position &geo) const
returns height for the given geo coordinate (WGS84)
void add(double x, double y, double z=0)
Makes the boundary include the given coordinate.
const Boundary & getBoundary()
returns the convex boundary of all known triangles
double ymax() const
Returns maximum y-coordinate.
bool ready() const
returns whether the NBHeightMapper has data
Set z-values for all network positions based on data from a height map.
void sub(double dx, double dy)
Substracts the given position from this one.
void endProcessMsg(std::string msg)
Ends a process information.