Added georeferencing

src/CellMap.cpp

@@ -4,13 +4,13 @@
 #include "contour_creator.hh"
 
 
-CellMap::CellMap(int width, int height, uint8_t* cells, OGRSpatialReference reference_system)
+CellMap::CellMap(int width, int height, uint8_t* cells, OGRSpatialReference reference_system, double* geotransform)
 {
     this->width = width;
     this->height = height;
     this->cells = cells;
     this->reference_system = reference_system;
-
+    this->geotransform = geotransform;
 }
 int CellMap::get_cell(int x, int y)
 {

src/HeightMap.cpp

@@ -1,3 +1,4 @@
+#include <gdal/cpl_conv.h>
 #include <iostream>
 #include <stdexcept>
 
@@ -28,8 +29,12 @@ HeightMap::HeightMap(const char* filepath)
     this->height = band->GetYSize();
     this->min = band->GetMinimum();
     this->max = band-> GetMaximum();
+    //https://gdal.org/api/gdaldataset_cpp.html#_CPPv4N11GDALDataset15GetGeoTransformEPd
+    this->geotransform = (double *) CPLMalloc(sizeof(double)*6);
     this->reference_system = *(file->GetSpatialRef());
 
+    file->GetGeoTransform(this->geotransform);
+    
     this->data = (float *) CPLMalloc(sizeof(float)*width*height);
     CPLErr error = band->RasterIO( GF_Read, 0, 0, width, height,
                     this->data, width, height, GDT_Float32,

src/contour_creator.hh

@@ -9,6 +9,8 @@ class HeightMap
 {
     public:
         float* data;
+        double* geotransform; //!< Six double buffer for storing the affine transformations
+        // https://gdal.org/api/gdaldataset_cpp.html#_CPPv4N11GDALDataset15GetGeoTransformEPd
         int width; //!< width of image in pixels
         int height; //!< height of image in pixels
         float min; //!< Minimum value in image
@@ -26,10 +28,10 @@ class CellMap
 {
     public:
         uint8_t* cells; //!< pointer to the first cell in the array. uint8_t is a 8 bit unsigned integer 
+        double* geotransform; //!< Six double buffer for storing the affine transformations
         int width; //!< width of image in cells
         int height; //!< height of image in cells
         OGRSpatialReference reference_system;
-
+        CellMap(int width, int height, uint8_t* cells, OGRSpatialReference reference_system, double* geotransform);
-        CellMap(int width, int height, uint8_t* cells, OGRSpatialReference reference_system);
         int get_cell(int x,int y);
 };

src/main.cpp

@@ -24,7 +24,7 @@ CellMap produce_cellmap(HeightMap* heightmap, float z)
         *(cells + i) = result;
     }
     
-    return CellMap(heightmap->width-1, heightmap->height-1, cells, heightmap->reference_system);
+    return CellMap(heightmap->width-1, heightmap->height-1, cells, heightmap->reference_system, heightmap->geotransform);
 
 
 }
@@ -50,6 +50,16 @@ std::vector<CellMap> vector_cellmap(HeightMap* heightmap, int interval)
         return vector_contours;
 }
 
+std::tuple<double, double> local_to_projected(double* geotransform, int x, int y)
+{
+    return {
+            geotransform[1] * x + geotransform[2] * y +
+            geotransform[1] * 0.5 + geotransform[2] * 0.5 + geotransform[0],
+            geotransform[4] * x + geotransform[5] * y +
+            geotransform[4] * 0.5 + geotransform[5] * 0.5 + geotransform[3]
+            };
+}
+
 void write_output_file(std::vector<CellMap> cellmaps, const char *filepath)
 {
 
@@ -98,13 +108,7 @@ void write_output_file(std::vector<CellMap> cellmaps, const char *filepath)
         for (int i = 0; i<cellmap->height*cellmap->width; i++) {
             if (*(cellmap->cells + i) != 0 && *(cellmap->cells + i) != 15)
             {
-                int x_int = i%cellmap->width;
-                int y_int = cellmap->height*cellmap->width - i/cellmap->width;
-                double x = double(x_int);
-                double y = double(y_int);
-
                 OGRFeature *poFeature;
-
                 poFeature = OGRFeature::CreateFeature( poLayer->GetLayerDefn() );
                 poFeature->SetField( "Name", *(cellmap->cells + i) );
 
@@ -117,6 +121,11 @@ void write_output_file(std::vector<CellMap> cellmaps, const char *filepath)
                     printf( "Transformation failed.\n" );
                 */
                 OGRPoint pt;
+
+                int x_raw = i%cellmap->width;
+                int y_raw = i/cellmap->width;
+                auto [x, y] = local_to_projected(cellmap->geotransform, x_raw, y_raw);
+
                 pt.setX( x );
                 pt.setY( y );