37 #pragma GCC diagnostic push
38 #pragma GCC diagnostic ignored "-Wpedantic"
40 #include <ogrsf_frmts.h>
42 #include <gdal_priv.h>
44 #pragma GCC diagnostic pop
83 WRITE_WARNING(
"Cannot supply height since no height data was loaded");
87 const Boundary& boundary = item.first;
88 int16_t* raster = item.second;
90 if (boundary.
around(geo)) {
95 corners.push_back(
Position(floor(normX) + 0.5, floor(normY) + 0.5, raster[(int)normY * xSize + (
int)normX]));
96 if (normX - floor(normX) > 0.5) {
97 corners.push_back(
Position(floor(normX) + 1.5, floor(normY) + 0.5, raster[(int)normY * xSize + (
int)normX + 1]));
99 corners.push_back(
Position(floor(normX) - 0.5, floor(normY) + 0.5, raster[(int)normY * xSize + (
int)normX - 1]));
101 if (normY - floor(normY) > 0.5) {
102 corners.push_back(
Position(floor(normX) + 0.5, floor(normY) + 1.5, raster[((int)normY + 1) * xSize + (int)normX]));
104 corners.push_back(
Position(floor(normX) + 0.5, floor(normY) - 0.5, raster[((int)normY - 1) * xSize + (int)normX]));
108 if (result > -1e5 && result < 1e5) {
115 minB[0] = (float)geo.
x() - 0.00001f;
116 minB[1] = (float)geo.
y() - 0.00001f;
117 maxB[0] = (float)geo.
x() + 0.00001f;
118 maxB[1] = (float)geo.
y() + 0.00001f;
120 int hits =
myRTree.Search(minB, maxB, queryResult);
122 assert(hits == (
int)result.size());
125 for (Triangles::iterator it = result.begin(); it != result.end(); it++) {
128 return triangle->
getZ(geo);
141 const float cmin[2] = {(float) b.
xmin(), (float) b.
ymin()};
142 const float cmax[2] = {(float) b.
xmax(), (float) b.
ymax()};
143 myRTree.Insert(cmin, cmax, triangle);
149 if (oc.
isSet(
"heightmap.geotiff")) {
151 std::vector<std::string> files = oc.
getStringVector(
"heightmap.geotiff");
152 for (std::vector<std::string>::const_iterator file = files.begin(); file != files.end(); ++file) {
156 " done (parsed " +
toString(numFeatures) +
160 if (oc.
isSet(
"heightmap.shapefiles")) {
162 std::vector<std::string> files = oc.
getStringVector(
"heightmap.shapefiles");
163 for (std::vector<std::string>::const_iterator file = files.begin(); file != files.end(); ++file) {
167 " done (parsed " +
toString(numFeatures) +
177 #if GDAL_VERSION_MAJOR < 2
179 OGRDataSource* ds = OGRSFDriverRegistrar::Open(file.c_str(), FALSE);
182 GDALDataset* ds = (GDALDataset*)GDALOpenEx(file.c_str(), GDAL_OF_VECTOR | GA_ReadOnly, NULL, NULL, NULL);
185 throw ProcessError(
"Could not open shape file '" + file +
"'.");
189 OGRLayer* layer = ds->GetLayer(0);
190 layer->ResetReading();
194 OGRSpatialReference* sr_src = layer->GetSpatialRef();
195 OGRSpatialReference sr_dest;
196 sr_dest.SetWellKnownGeogCS(
"WGS84");
197 OGRCoordinateTransformation* toWGS84 = OGRCreateCoordinateTransformation(sr_src, &sr_dest);
199 WRITE_WARNING(
"Could not create geocoordinates converter; check whether proj.4 is installed.");
204 layer->ResetReading();
205 while ((feature = layer->GetNextFeature()) != NULL) {
206 OGRGeometry* geom = feature->GetGeometryRef();
210 assert(std::string(geom->getGeometryName()) == std::string(
"POLYGON"));
212 geom->transform(toWGS84);
213 OGRLinearRing* cgeom = ((OGRPolygon*) geom)->getExteriorRing();
215 assert(cgeom->getNumPoints() == 4);
217 for (
int j = 0; j < 3; j++) {
218 Position pos((
double) cgeom->getX(j), (
double) cgeom->getY(j), (
double) cgeom->getZ(j));
219 corners.push_back(pos);
256 OGRFeature::DestroyFeature(feature);
258 #if GDAL_VERSION_MAJOR < 2
259 OGRDataSource::DestroyDataSource(ds);
263 OCTDestroyCoordinateTransformation(
reinterpret_cast<OGRCoordinateTransformationH
>(toWGS84));
268 WRITE_ERROR(
"Cannot load shape file since SUMO was compiled without GDAL support.");
278 GDALDataset* poDataset = (GDALDataset*)GDALOpen(file.c_str(), GA_ReadOnly);
279 if (poDataset == 0) {
284 const int xSize = poDataset->GetRasterXSize();
285 const int ySize = poDataset->GetRasterYSize();
286 double adfGeoTransform[6];
287 if (poDataset->GetGeoTransform(adfGeoTransform) == CE_None) {
288 Position topLeft(adfGeoTransform[0], adfGeoTransform[3]);
292 boundary.
add(topLeft);
293 boundary.
add(topLeft.
x() + horizontalSize, topLeft.
y() + verticalSize);
295 WRITE_ERROR(
"Could not parse geo information from " + file +
".");
298 const int picSize = xSize * ySize;
299 int16_t* raster = (int16_t*)CPLMalloc(
sizeof(int16_t) * picSize);
300 for (
int i = 1; i <= poDataset->GetRasterCount(); i++) {
301 GDALRasterBand* poBand = poDataset->GetRasterBand(i);
302 if (poBand->GetColorInterpretation() != GCI_GrayIndex) {
303 WRITE_ERROR(
"Unknown color band in " + file +
".");
307 if (poBand->GetRasterDataType() != GDT_Int16) {
312 assert(xSize == poBand->GetXSize() && ySize == poBand->GetYSize());
313 if (poBand->RasterIO(GF_Read, 0, 0, xSize, ySize, raster, xSize, ySize, GDT_Int16, 0, 0) == CE_Failure) {
319 GDALClose(poDataset);
320 myRasters.push_back(std::make_pair(boundary, raster));
324 WRITE_ERROR(
"Cannot load GeoTIFF file since SUMO was compiled without GDAL support.");
338 CPLFree(item.second);
364 return myCorners.around(pos);
382 Position side1 = myCorners[1] - myCorners[0];
383 Position side2 = myCorners[2] - myCorners[0];