Started on function to write output file

src/CaseMap.cpp

@@ -1,11 +1,14 @@
 #include <cstdint>
+#include <gdal/ogr_spatialref.h>
+#include <ogr_spatialref.h>
 #include "CaseMap.hh"
 
-CaseMap::CaseMap(int width, int height, uint8_t* cases)
+CaseMap::CaseMap(int width, int height, uint8_t* cases, OGRSpatialReference reference_system)
 {
     this->width = width;
     this->height = height;
     this->cases = cases;
+    this->reference_system = reference_system;
 
 }
 int CaseMap::get_case(int x, int y)

src/CaseMap.hh

@@ -1,14 +1,17 @@
+#include <gdal/ogr_spatialref.h>
+#include <ogr_spatialref.h>
+#include <cstdint>
 /**
     @brief stores the cases from marching squars for one elevation level
  */
-#include <cstdint>
 class CaseMap
 {
     public:
         uint8_t* cases; //!< pointer to the first case in the array. uint8_t is a 8 bit unsigned integer 
         int width; //!< width of image in cases
         int height; //!< height of image in cases
+        OGRSpatialReference reference_system;
 
-        CaseMap(int width, int height, uint8_t* cases);
+        CaseMap(int width, int height, uint8_t* cases, OGRSpatialReference reference_system);
         int get_case(int x,int y);
 };

src/HeightMap.cpp

@@ -1,3 +1,4 @@
+#include <iostream>
 #include <stdexcept>
 
 #include <gdal/gdal.h>
@@ -26,6 +27,7 @@ HeightMap::HeightMap(const char* filepath)
     this->height = band->GetYSize();
     this->min = band->GetMinimum();
     this->max = band-> GetMaximum();
+    this->reference_system = *(file->GetSpatialRef());
     
     this->data = (float *) CPLMalloc(sizeof(float)*width*height);
     CPLErr error = band->RasterIO( GF_Read, 0, 0, width, height,
@@ -38,5 +40,6 @@ float HeightMap::get_pixel(int x, int y)
 {
     // all the pixels are in an array of floats from left to right, top to bottom
     int offset = ((this->width * y) + x);
+    std::cout << " " << x << ","<< y;
     return *(this->data + offset);
 } 

src/HeightMap.hh

@@ -1,3 +1,4 @@
+#include <gdal/ogr_spatialref.h>
 /**
     @brief stores the contents of a tif file with float32 values
  */
@@ -9,6 +10,7 @@ class HeightMap
         int height; //!< height of image in pixels
         float min; //!< Minimum value in image
         float max; //!< Maximum value in image
+        OGRSpatialReference reference_system;
 
         HeightMap(const char* filepath);
         float get_pixel(int x,int y);

src/main.cpp

@@ -1,5 +1,7 @@
 #include "HeightMap.hh"
 #include "CaseMap.hh"
+#include <gdal/ogr_api.h>
+#include "ogrsf_frmts.h"
 #include <iostream>
 #include <ostream>
 #include <vector>
@@ -12,8 +14,8 @@ CaseMap produce_casemap(HeightMap* heightmap, float z)
 {
     uint8_t *cases = (uint8_t *) CPLMalloc(sizeof(uint8_t)*(heightmap->width-1)*(heightmap->height-1));
     for (int i = 0; i<(heightmap->width-1)*(heightmap->height-1); i++) {
-        int y = i/(heightmap->width-1);
-        int x = i%(heightmap->width-1);
+        int y = i/(heightmap->height-1);
+        int x = i%(heightmap->height-1);
         uint8_t result = (heightmap->get_pixel(x,y+1)>z) + 
         (heightmap->get_pixel(x+1,y+1)>z)*2 +
         (heightmap->get_pixel(x+1,y)>z)*4 +
@@ -21,11 +23,103 @@ CaseMap produce_casemap(HeightMap* heightmap, float z)
         *(cases + i) = result;
     }
     
-    return CaseMap(heightmap->width-1, heightmap->height-1, cases);
+    return CaseMap(heightmap->width-1, heightmap->height-1, cases, heightmap->reference_system);
 
 
 }
 
+void write_output_file(CaseMap* casemap)
+{
+
+    const char *pszDriverName = "ESRI Shapefile";
+    GDALDriver *poDriver;
+
+    GDALAllRegister();
+
+    poDriver = GetGDALDriverManager()->GetDriverByName(pszDriverName );
+    if( poDriver == NULL )
+    {
+        printf( "%s driver not available.\n", pszDriverName );
+        exit( 1 );
+    }
+
+    GDALDataset *poDS;
+
+    poDS = poDriver->Create( "point_out.shp", 0, 0, 0, GDT_Unknown, NULL );
+    if( poDS == NULL )
+    {
+        printf( "Creation of output file failed.\n" );
+        exit( 1 );
+    }
+
+    OGRLayer *poLayer;
+
+    poLayer = poDS->CreateLayer( "point_out", &casemap->reference_system, wkbPoint, NULL );
+    if( poLayer == NULL )
+    {
+        printf( "Layer creation failed.\n" );
+        exit( 1 );
+    }
+
+    OGRFieldDefn oField( "Name", OFTString );
+
+    oField.SetWidth(32);
+
+    if( poLayer->CreateField( &oField ) != OGRERR_NONE )
+    {
+        printf( "Creating Name field failed.\n" );
+        exit( 1 );
+    }
+
+    for (int i = 0; i<casemap->height*casemap->width; i++) {
+        if (*(casemap->cases + i) != 0 && *(casemap->cases + i) != 15)
+        {
+            int x_int = i%casemap->height;
+            int y_int = casemap->height*casemap->width - i/casemap->height;
+            double x = double(x_int);
+            double y = double(y_int);
+
+            OGRFeature *poFeature;
+
+            poFeature = OGRFeature::CreateFeature( poLayer->GetLayerDefn() );
+            poFeature->SetField( "Name", *(casemap->cases + i) );
+
+            //OGRSpatialReference local;
+
+            //auto poCT = OGRCreateCoordinateTransformation( &local, &casemap->reference_system );
+            
+            /*
+            if( poCT == NULL || !poCT->Transform( 1, &x, &y ) )
+                printf( "Transformation failed.\n" );
+            */
+            OGRPoint pt;
+            pt.setX( x );
+            pt.setY( y );
+
+            poFeature->SetGeometry( &pt );
+
+            if( poLayer->CreateFeature( poFeature ) != OGRERR_NONE )
+            {
+                printf( "Failed to create feature in shapefile.\n" );
+                exit( 1 );
+            }
+
+            OGRFeature::DestroyFeature( poFeature );
+        }
+    }
+
+    GDALClose( poDS );
+    /*
+    OGRDataSourceH datasource = OGR_Dr_CreateDataSource(OGRGetDriverByName("GeoJSON"), "contour.geojson", new char*);
+    OGRSpatialReference reference_system = casemap->reference_system;
+    OGRLayerH layer = OGR_DS_CreateLayer(datasource, "contour", &reference_system, wkbLineString25D, new char*);
+    
+    auto feature = OGR_F_Create(OGR_L_GetLayerDefn(static_cast<OGRLayerH>(layer)));
+
+    OGRGeometryH geometry = OGR_G_CreateGeometry(wkbLineString25D);
+    */
+}
+
 int main(int argc, const char* argv[])
 {
     const char* filepath = argv[1];
@@ -35,7 +129,7 @@ int main(int argc, const char* argv[])
     std::cout << "max: " << map.max << "    min: " << map.min << "\n";
     
     auto casemap = produce_casemap(&map, 40.0);
-    
+    /*
     for (int y = 0; y < casemap.height; y++)
     {
         for (int x = 0; x < casemap.width; x++)
@@ -46,4 +140,6 @@ int main(int argc, const char* argv[])
         }
        //std::cout << "\n";
     }
+    */
+    write_output_file(&casemap);
 }