open a geotiff dtm with gdal

DTM_GEOTIFF_GDAL.cpp

// Compile with  g++ geotiff2.cpp -std=c++11 -lgdal && ./a.out

#include "cpl_conv.h" // for CPLMalloc()
#include "gdal_priv.h"
#include <iostream>

typedef struct {
  double lon;
  double lat;
} coord;

class TerrainTile {
public:
  TerrainTile();
  ~TerrainTile();

  TerrainTile(const char *);
  // TerrainTile(QJsonDocument document);
  // TerrainTile(QByteArray byteArray);

  bool isIn(coord coordinate);

  bool isValid(void) const { return _isValid; }

  double elevation(coord *coordinate);

  double minElevation(void) const { return _minElevation; }

  double maxElevation(void) const { return _maxElevation; }

  double avgElevation(void) const { return _avgElevation; }

  coord centerCoordinate(void);

  // static QByteArray serialize(QByteArray input);

  /// Approximate spacing of the elevation data measurement points
  static constexpr double terrainAltitudeSpacing = 30.0;

private:
  double _isValid;
  double _minElevation;
  double _maxElevation;
  double _avgElevation;

  int lonlatToxy(coord *c);
  int xyTolonlat(coord *c);
  int lonlatToPixel(coord c, int *xy);

  GDALDataset *poDataset;
  GDALRasterBand *poBand;
  double adfGeoTransform[6];

  OGRCoordinateTransformation *lonlatToxyTransformation;
  OGRCoordinateTransformation *xyTolonlatTransformation;
};

bool TerrainTile::isIn(coord coordinate) {
  int xy[2];
  lonlatToPixel(coordinate, xy);

  if (xy[0] < poBand->GetXSize() && xy[1] < poBand->GetYSize()) {
    int ret = poBand->GetDataCoverageStatus(xy[0], xy[1], 1, 1, 0, nullptr);
    if (ret == GDAL_DATA_COVERAGE_STATUS_DATA){
        // If data is NoDataValue, return false
        return (elevation(&coordinate)!=poBand->GetNoDataValue());
    }
  }
  return false;
}

coord TerrainTile::centerCoordinate(void){
    coord result;
    double xSize = poDataset->GetRasterXSize();
    double ySize = poDataset->GetRasterYSize();

    // Fill coord with XY data before transformation
   result.lon =  adfGeoTransform[0]+ xSize*adfGeoTransform[1]/2;
   result.lat =  adfGeoTransform[3]+ xSize*adfGeoTransform[5]/2;

   xyTolonlat(&result);
   return result;

}

int TerrainTile::lonlatToxy(coord *c) {
  int ret = lonlatToxyTransformation->Transform(1, &(c->lon), &(c->lat), NULL);
  return ret;
}

int TerrainTile::xyTolonlat(coord *c) {
  int ret = xyTolonlatTransformation->Transform(1, &(c->lon), &(c->lat), NULL);
  return ret;
}

int TerrainTile::lonlatToPixel(coord c, int *xy) {
  double lon = c.lon;
  double lat = c.lat;
  int ret = lonlatToxyTransformation->Transform(1, &lon, &lat, NULL);
  int x, y;

  xy[0] = (int)(((lon - adfGeoTransform[0]) / adfGeoTransform[1]));
  xy[1] = (int)(((lat - adfGeoTransform[3]) / adfGeoTransform[5]));

  lon = x;
  lat = y;

  return ret;
}

TerrainTile::~TerrainTile() { GDALClose(poDataset); }

TerrainTile::TerrainTile(const char *fname) {
  GDALAllRegister();
  poDataset = (GDALDataset *)GDALOpen("dsm.tif", GA_ReadOnly);
  if (poDataset == NULL) {
    _isValid = false;
  }

  int x;
  int y;

  if (poDataset->GetProjectionRef() != NULL) {
    _isValid = false;
  }
  if (poDataset->GetGeoTransform(adfGeoTransform) != CE_None) {
    _isValid = false;
  }

  OGRSpatialReference src;
  src.SetWellKnownGeogCS("WGS84");
  OGRSpatialReference dest(poDataset->GetProjectionRef());
  lonlatToxyTransformation = OGRCreateCoordinateTransformation(&src, &dest);
  xyTolonlatTransformation = OGRCreateCoordinateTransformation(&dest, &src);

  // Fetch a raster band

  int nBlockXSize, nBlockYSize;
  int bGotMin, bGotMax;
  double adfMinMax[2];
  poBand = poDataset->GetRasterBand(1);

  poBand->ComputeStatistics(true, &_minElevation, &_maxElevation, &_avgElevation, NULL, NULL, NULL);

}

double TerrainTile::elevation(coord *coordinate) {
  int xy[2];
  lonlatToPixel(*coordinate, xy);

  coord curpoint;
  curpoint.lon = coordinate->lon;
  curpoint.lat = coordinate->lat;

  if(!(xy[0] < poBand->GetXSize() && xy[1] < poBand->GetYSize())) {
      std::cout << "elevation error" << std::endl;
      return poBand->GetNoDataValue(nullptr);
  }
  // Reading raster data
  float *pafScanline;
  int nXSize = poBand->GetXSize();
  pafScanline = (float *)CPLMalloc(sizeof(float) * nXSize);
  poBand->RasterIO(GF_Read, 0, xy[1], nXSize, 1, pafScanline, nXSize, 1,
                   GDT_Float32, 0, 0);

  return pafScanline[xy[0]];
  CPLFree(pafScanline);
}

int main() {

  TerrainTile test("dtm.tif");

  coord c;

  // This point is in the tile

  c.lon = 174.941589;
  c.lat = -37.815401;

  std::cout << test.isIn(c) << std::endl;

  //  This point is nowhere
  c.lon = 175.940404;
  c.lat = -37.815851;

  std::cout << test.isIn(c) << std::endl;

  //  This point is in the window, not available
  c.lon = 174.940404;
  c.lat = -37.815851;

  std::cout << test.isIn(c) << std::endl;

  coord c2 = test.centerCoordinate();

  std::cout << test.isIn(c2) << std::endl;

  std::cout << " " << c2.lon << " " << c2.lat << std::endl;

  std::cout << "Center :\t"<<test.elevation(&c2) << " m" << std::endl;
  std::cout << "Outside Window \t" << test.elevation(&c) << " m" << std::endl;
  std::cout << "Min Elev " << test.minElevation() << " m" << std::endl;
  std::cout << "Max Elev " << test.maxElevation() << " m" << std::endl;
  std::cout << "Avg Elev " << test.avgElevation() << " m" << std::endl;
  std::cout << "Center " << test.centerCoordinate().lon << ","
            << test.centerCoordinate().lat << std::endl;

}