Convert GeoTiff height map to 16-bit png

parseGeoTiff.cpp

// Converts geotiff height map to 16-bit png
// g++ parseGeoTiff.cpp -lgdal -lpng -o parseGeoTiff
// use: ./parseGeoTiff input.tif output.png
// License: Simplified BSD License
#include <iostream>
#include <string>
#include <vector>
#include <algorithm>
#include <gdal/gdal.h>
#include <gdal/gdal_priv.h>
#include <gdal/cpl_conv.h>
#include <png.h>

void saveElevationAsPNG(const std::string& filename, const std::vector<double>& elevationData, int width, int height) {
    FILE* fp = fopen(filename.c_str(), "wb");
    if (!fp) {
        std::cerr << "Error opening PNG file for writing: " << filename << std::endl;
        return;
    }

    png_structp png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, nullptr, nullptr, nullptr);
    if (!png_ptr) {
        std::cerr << "Error creating PNG write struct" << std::endl;
        fclose(fp);
        return;
    }

    png_infop info_ptr = png_create_info_struct(png_ptr);
    if (!info_ptr) {
        std::cerr << "Error creating PNG info struct" << std::endl;
        png_destroy_write_struct(&png_ptr, nullptr);
        fclose(fp);
        return;
    }

    if (setjmp(png_jmpbuf(png_ptr))) {
        std::cerr << "Error during PNG init_io" << std::endl;
        png_destroy_write_struct(&png_ptr, &info_ptr);
        fclose(fp);
        return;
    }

    png_init_io(png_ptr, fp);

    // Set PNG image information
    png_set_IHDR(
        png_ptr, info_ptr, width, height,
        16, PNG_COLOR_TYPE_GRAY, PNG_INTERLACE_NONE,
        PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT
    );
    png_write_info(png_ptr, info_ptr);

    // Set byte order to little-endian
    png_uint_16 le = 0x0001;
    // Test byte order first
    if ((*(png_const_bytep) & le) != 0)
       png_set_swap(png_ptr);

    // Normalize elevation data to fit within 16-bit range (0-65535)
    std::vector<png_uint_16> normalizedElevationData(elevationData.size());
    double minElevation = *std::min_element(elevationData.begin(), elevationData.end());
    double maxElevation = *std::max_element(elevationData.begin(), elevationData.end());

    // If the minimum elevation is negative, find the next positive elevation manually
    if(minElevation < 0.0) {
        minElevation = 99999.0;
        for(size_t i = 0; i < elevationData.size(); ++i) {
            if(elevationData[i] < minElevation && elevationData[i] > 0.0) {
                minElevation = elevationData[i];                
            }
        }
    }

    std::cout << "Min elevation: " << minElevation << std::endl;
    std::cout << "Max elevation: " << maxElevation << std::endl;
  
    
    for (size_t i = 0; i < elevationData.size(); ++i) {
        double value = (std::max(0.0,std::min(1.0,(elevationData[i]-minElevation) / (maxElevation - minElevation))));
        normalizedElevationData[i] = static_cast<png_uint_16>(value*65535.0);
    }
    
    // Write image data row by row 
    std::vector<png_bytep> row_pointers(height);
    for (int y = 0; y < height; y++)
    {
        row_pointers[y] = (png_bytep)&normalizedElevationData[y * width];
    }
    png_write_image(png_ptr, row_pointers.data());
    
    png_write_end(png_ptr, nullptr);

    png_destroy_write_struct(&png_ptr, &info_ptr);
    fclose(fp);
    std::cout << "Elevation data saved as " << filename << std::endl;
}

int main(int argc, char* argv[]) {

  // Read input and output file paths from command line arguments
  if (argc != 3) {
    std::cout << "Usage: " << argv[0] << " <input_geotiff_file> <output_png_file>" << std::endl;
    return 1;
  }


  // Specify the path to the geotiff file
  std::string filePath(argv[1]);

  // Register GDAL drivers
  GDALAllRegister();

  // Open the geotiff file
  GDALDataset* dataset = (GDALDataset*)GDALOpen(filePath.c_str(), GA_ReadOnly);
  if (dataset == nullptr) {
    std::cout << "Failed to open the geotiff file." << std::endl;
    return 1;
  }

  // Get the number of raster bands in the geotiff file
  int numBands = dataset->GetRasterCount();
  std::cout << "Number of raster bands: " << numBands << std::endl;

  std::vector<double> elevationData;
  int width,height;
  // Access the raster data for each band
  for (int bandIndex = 1; bandIndex <= numBands; ++bandIndex) {
    GDALRasterBand* band = dataset->GetRasterBand(bandIndex);

    // Get the size of the raster band
    width = band->GetXSize();
    height = band->GetYSize();
    std::cout << "Raster band " << bandIndex << " size: " << width << " x " << height << std::endl;

    // Read the raster data
    double* data = new double[width * height];
    band->RasterIO(GF_Read, 0, 0, width, height, data, width, height, GDT_Float64, 0, 0);

    // Process the raster data here...
    for (int i = 0; i < width * height; ++i) {
      elevationData.push_back(data[i]);
    }

    // Clean up
    delete[] data;
  }

  // Save the elevation data as a PNG image
  saveElevationAsPNG(std::string(argv[2]), elevationData, width, height);

  // Close the geotiff file
  GDALClose(dataset);

  return 0;
}