diff --git a/src/CaseMap.cpp b/src/CaseMap.cpp index 0de6d8b..cc23b86 100644 --- a/src/CaseMap.cpp +++ b/src/CaseMap.cpp @@ -1,11 +1,14 @@ #include +#include +#include #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) diff --git a/src/CaseMap.hh b/src/CaseMap.hh index 454202d..e0a5a59 100644 --- a/src/CaseMap.hh +++ b/src/CaseMap.hh @@ -1,14 +1,17 @@ +#include +#include +#include /** @brief stores the cases from marching squars for one elevation level */ -#include 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); }; \ No newline at end of file diff --git a/src/HeightMap.cpp b/src/HeightMap.cpp index b82e96e..74eb37c 100644 --- a/src/HeightMap.cpp +++ b/src/HeightMap.cpp @@ -1,3 +1,4 @@ +#include #include #include @@ -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); } diff --git a/src/HeightMap.hh b/src/HeightMap.hh index b10c9d7..94d2112 100644 --- a/src/HeightMap.hh +++ b/src/HeightMap.hh @@ -1,3 +1,4 @@ +#include /** @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); diff --git a/src/main.cpp b/src/main.cpp index 8fd656d..a936309 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,5 +1,7 @@ #include "HeightMap.hh" #include "CaseMap.hh" +#include +#include "ogrsf_frmts.h" #include #include #include @@ -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; iheight*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(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); } \ No newline at end of file