netconvert: read --heightmap.geotiff spatial reference and transform coordinates for lookup

This adds support for geotiff input data in UTM projection, among
others, and avoids an additional step to reproject the geotiff data to
WGS84 (by using `gdalwarp`, for example), which adds a resampling step
that introduces inaccuracy.

Author: opatut

src/netbuild/NBHeightMapper.cpp

@@ -92,10 +92,28 @@ NBHeightMapper::getZ(const Position& geo) const {
         const Boundary& boundary = item.boundary;
         float* raster = item.raster;
         double result = -1e6;
-        if (boundary.around2D(geo)) {
+
+        double x = geo.x();
+        double y = geo.y(); 
+
+#ifdef HAVE_GDAL
+        // Transform geo coordinates to the coordinate system of this
+        // raster image for lookup in its raster, if applicable.
+        if (item.transform != nullptr) {
+            // Since the input coordinates are always WGS84 (they may be
+            // transformed to it in NBNetBuilder::transformCoordinate), and
+            // WGS84 uses latitude-longitude order (y-x), we have to swap the
+            // input coordinates here.
+            std::swap(x, y);
+
+            item.transform->Transform(1, &x, &y);
+        }
+#endif
+
+        if (boundary.around2D(x, y)) {
             const int xSize = item.xSize;
-            const double normX = (geo.x() - boundary.xmin()) / mySizeOfPixel.x();
-            const double normY = (geo.y() - boundary.ymax()) / mySizeOfPixel.y();
+            const double normX = (x - boundary.xmin()) / mySizeOfPixel.x();
+            const double normY = (y - boundary.ymax()) / mySizeOfPixel.y();
             PositionVector corners;
             corners.push_back(Position(floor(normX) + 0.5, floor(normY) + 0.5, raster[(int)normY * xSize + (int)normX]));
             if (normX - floor(normX) > 0.5) {
@@ -332,17 +350,28 @@ NBHeightMapper::loadTiff(const std::string& file) {
         min = MIN2(min, (double)raster[i]);
         max = MAX2(max, (double)raster[i]);
     }
+
+    // Make a copy, GDALClose will destroy the original
+    OGRSpatialReference spatialRef(*poDataset->GetSpatialRef());
+
     GDALClose(poDataset);
+
     if (ok) {
         WRITE_MESSAGE("Read geotiff heightmap with size " + toString(xSize) + "," + toString(ySize)
                       + " for geo boundary [" + toString(boundary)
                       + "] with elevation range [" + toString(min) + "," + toString(max) + "].");
-        RasterData rasterData;
-        rasterData.raster = raster;
-        rasterData.boundary = boundary;
-        rasterData.xSize = xSize;
-        rasterData.ySize = ySize;
-        myRasters.push_back(rasterData);
+
+        OGRSpatialReference wgs;
+        wgs.SetWellKnownGeogCS("WGS84");
+
+        myRasters.push_back(RasterData{
+            raster,
+            boundary,
+            xSize,
+            ySize,
+            OGRCreateCoordinateTransformation(&wgs, &spatialRef)
+        });
+
         return picSize;
     } else {
         return 0;
@@ -364,6 +393,9 @@ NBHeightMapper::clearData() {
 #ifdef HAVE_GDAL
     for (auto& item : myRasters) {
         CPLFree(item.raster);
+        if (item.transform != nullptr) {
+            delete item.transform;
+        }
     }
     myRasters.clear();
 #endif

src/netbuild/NBHeightMapper.h

@@ -34,6 +34,10 @@ typedef __int16 int16_t;
 #include <utils/geom/Boundary.h>
 #include <utils/common/UtilExceptions.h>
 
+#ifdef HAVE_GDAL
+#include <ogr_spatialref.h>
+#endif
+
 #define TRIANGLE_RTREE_QUAL RTree<NBHeightMapper::Triangle*, NBHeightMapper::Triangle, float, 2, NBHeightMapper::QueryResult>
 
 // ===========================================================================
@@ -128,6 +132,9 @@ class NBHeightMapper {
         Boundary boundary;
         int xSize;
         int ySize;
+#ifdef HAVE_GDAL
+        OGRCoordinateTransformation* transform;
+#endif
     };
 
     /// @brief the singleton instance