Постановка задачи⌗

Подготовим простенькую C++ модель, возвращающую глубину Мирового океана или высоту над уровнем моря для заданных географических координат (широта, долгота), а также строгую принадлежность данной точки суше либо океану.

Для глубин и высот будем использовать растровые датасеты (полученные CleanTOPO2, GEBCO, SRTM и прочими проектами), а для строгой классификации суша-вода - векторные карты геометрии суши и океанов (например, проект Natural Earth содержит много полезных физических карт Земли). Наличие двух принципиально разных датасетов обусловлено прежде всего факторами точности измерений: разрешения растра, топографической съёмки, грубостью векторных многоугольников и т.д. Для чтения датасетов используем библиотеку абстракции геопространственных данных GDAL, а для тестирования модели - графические возможности Qt.

Датасеты⌗

Для глубин океанов и высот над уровнем моря я воспользуюсь одним комбинированным растром Blue Earth Bathymetry. Это редактированная версия GEBCO (General Bathymetric Chart of the Oceans), дополненная топографией для всего Земного шара.

Для классификации вода-земля я воспользуюсь векторной картой 1:10м массы суши Земли от Natural Earth Data.

Наша модель, однако, сможет поддерживать раздельные растры для суши и океана, а также незаисимую работу двух векторных слоёв.

Интерфейс к модели⌗ 


#include <filesystem>

class hypsobathymetry_model {
 public:
  virtual ~hypsobathymetry_model() = default;

  /*
   * Query interface
   */
  virtual bool is_land(double lat, double lon) = 0;
  virtual bool is_ocean(double lat, double lon) = 0;
  virtual int value(double lat, double lon) = 0;

  /*
   * Attachments of data files
   */
  virtual void set_ocean_geometry(const std::filesystem::path& path) = 0;
  virtual void set_land_geometry(const std::filesystem::path& path) = 0;
  virtual void set_ocean_values_raster(const std::filesystem::path& path) = 0;
  virtual void set_land_values_raster(const std::filesystem::path& path) = 0;
};

Реализация модели⌗

Конкреция использует С++ биндинги к GDAL. Мы владеем несколькими объектами GDALDataset (тип полностью абстрагирует конкретные растровые/векторные форматы данных).


#include "hypsobathymetry_model.h"

#include <gdal/gdal.h>
#include <gdal/ogrsf_frmts.h>

#include <memory>

class gdal_hypsobathymetry_model : public hypsobathymetry_model {
 private:
  std::unique_ptr<GDALDataset> m_ocean_geometry_dataset;
  std::unique_ptr<GDALDataset> m_land_geometry_dataset;
  std::shared_ptr<GDALDataset> m_ocean_values_dataset;
  std::shared_ptr<GDALDataset> m_land_values_dataset;

 public:
  gdal_hypsobathymetry_model();
  virtual bool is_land(double lat, double lon) override;
  virtual bool is_ocean(double lat, double lon) override;
  virtual int value(double lat, double lon) override;
  virtual void set_ocean_geometry(const std::filesystem::path& path) override;
  virtual void set_land_geometry(const std::filesystem::path& path) override;
  virtual void set_ocean_values_raster(const std::filesystem::path& path) override;
  virtual void set_land_values_raster(const std::filesystem::path& path) override;
};

Инициализация GDAL⌗

В конструкторе gdal_hypsobathymetry_model() инициализируем библиотеку GDAL:

GDALAllRegister();

Открытие датасетов⌗

Окрытие векторных и растровых датасетов происходит соответствующими вызовами GDALOpenEx() и GDALOpen():

GDALDataset* dataset;

dataset = static_cast<GDALDataset*>(
  GDALOpenEx(path.c_str(), GDAL_OF_VECTOR, nullptr, nullptr, nullptr));
  
dataset = static_cast<GDALDataset*>(GDALOpen(path.c_str(), GA_ReadOnly));

Извлечение информации из датасетов⌗

В GDAL применяется концепция слоёв (layer) для векторных и растровых областей (raster band) для растровых датасетов.

Вектор: геометрия суши/океанов⌗

Внутри датасета может быть один или более объектов OGRLayer. Внутри каждого такого слоя содержится коллекция обектов OGRFeature, каждый со своей геометрией OGRGeometry (фундаментальный тип геометрии един для всего слоя). Для тестирования “попадания” или нет заданных широты/долготы в пределы геометрии одного или более feature мы используем конструкцию следующего вида:

OGRPoint point(lon, lat);
for (auto&& layer : m_land_geometry_dataset->GetLayers()) {
  for (auto&& feature : layer) {
    auto geometry = feature->GetGeometryRef();
    if (geometry->Contains(&point) == true) {
      return true;
    }
  }
}

Как видно, мы создаём точку и тестируем её принадлежность каждой из извлекаемых геометрий методом OGRGeometry::Contains().

Растр: численные значения глубин океана и рельефа возвышенностей⌗

Растр читать немного сложнее, так как тут уже не уйти от математической трансформации географических координат: если вершины наших векторных геометрий хранятся в файле прямо в географических координатах (первая проекция, т.е. X как широта, Y как долгота), то прямоугольный растровый файл имеет чёткое количество пикселей по горизонтале и по вертикале. Необходимо сопоставлять точку, т.е. строку и столбец растра реальным географическим координатам и обратно.

Геотрансформация⌗

В терминологии GDAL смещение в растре по горизонтали называется “pixel offset”, по вертикали “line offset”. Сначала из загруженного датасета методом GetGeoTransform() считывается матрица прямого преобразования (маппинг смещений от верхнего левого угла растра в геокоординаты), затем функцией GDALInvGeoTransform() из этой матрицы строится обратная ей.

double geotransform[6], inverse_geotransform[6];
if (dataset->GetGeoTransform(geotransform) != CE_None) {
  throw std::runtime_error("Geo transform not present in the raster");
}

if (!GDALInvGeoTransform(geotransform, inverse_geotransform)) {
  throw std::domain_error("Can't invert the geo transform");
}

int raster_pixel_coord = static_cast<int>(
    floor(inverse_geotransform[0] + inverse_geotransform[1] * lon +
          inverse_geotransform[2] * lat));
int raster_line_coord = static_cast<int>(
    floor(inverse_geotransform[3] + inverse_geotransform[4] * lon +
          inverse_geotransform[5] * lat));
Считывание блока растра⌗

Внутри одного датасета может быть один или более объектов GDALRasterBand (аналогия с OGRLayer выше). Их порядковый счёт идёт с единицы. Извлекаем первый доступный raster band и читаем блок собственно данных методом GetRasterBand(). GDAL оптимизирует чтение растра прямоугольными блоками, однако нас интересует лишь одна точка растра, поэтому в нашем случае размер блока составит 1 x 1 пиксель:

auto band = dataset->GetRasterBand(1);
double raster_pixel_value[2];
if (band->RasterIO(GF_Read, raster_pixel_coord, raster_line_coord, 1, 1,
                   raster_pixel_value, 1, 1, GDT_CFloat64, 0, 0) != CE_None) {
  throw std::runtime_error("Raster read error");
}

int value = (int)raster_pixel_value[0];

Закрытие датасетов⌗

Объекты GDALDataset закрывают все открытые файлы и полностью деинициализируются в деструкторе данного класса.

Полный код имплементации модели⌗

Ниже приводится полный исходный код. Модель пытается использовать всю доступную ей информацию: например, при отсутствии датасета геометрии океанов она применяет карту суши (и наоборот).


#include "gdal_hypsobathymetry_model.h"

#include <iostream>
#include <stdexcept>

static void check_lat_lon(double lat, double lon) {
  if (lat > 90.0 || lat < -90.0 || lon > 180.0 || lon < -180.0) {
    throw std::out_of_range("Invalid lat/lon");
  }
}

bool gdal_hypsobathymetry_model::is_land(double lat, double lon) {
  check_lat_lon(lat, lon);
  if (m_land_geometry_dataset == nullptr &&
      m_ocean_geometry_dataset == nullptr) {
    throw std::domain_error("No geometry data");
  }

  if (m_land_geometry_dataset == nullptr) {
    return !is_ocean(lat, lon);
  }

  OGRPoint point(lon, lat);
  for (auto&& layer : m_land_geometry_dataset->GetLayers()) {
    for (auto&& feature : layer) {
      auto geometry = feature->GetGeometryRef();
      if (geometry->Contains(&point) == true) {
        return true;
      }
    }
  }
  return false;
}
bool gdal_hypsobathymetry_model::is_ocean(double lat, double lon) {
  check_lat_lon(lat, lon);
  if (m_land_geometry_dataset == nullptr &&
      m_ocean_geometry_dataset == nullptr) {
    throw std::domain_error("No geometry data");
  }

  if (m_ocean_geometry_dataset == nullptr) {
    return !is_land(lat, lon);
  }

  OGRPoint point(lon, lat);
  for (auto&& layer : m_ocean_geometry_dataset->GetLayers()) {
    for (auto&& feature : layer) {
      auto geometry = feature->GetGeometryRef();
      if (geometry->Contains(&point) == true) {
        return true;
      }
    }
  }
  return false;
}
int gdal_hypsobathymetry_model::value(double lat, double lon) {
  check_lat_lon(lat, lon);
  if (m_ocean_values_dataset == nullptr && m_land_values_dataset == nullptr) {
    throw std::domain_error("No raster values data");
  }

  // select the appropriate raster dataset
  GDALDataset* dataset = nullptr;
  if (m_ocean_values_dataset == m_land_values_dataset) {
    // combined dataset
    dataset = m_ocean_values_dataset.get();
  } else {
    // both datasets are different; maybe one of them is null
    bool use_land = is_land(lat, lon);
    dataset =
        use_land ? m_land_values_dataset.get() : m_ocean_values_dataset.get();
    if (dataset == nullptr) {
      std::string msg = use_land ? "Land " : "Ocean ";
      msg.append("raster data is null");
      throw std::domain_error(msg);
    }
  }

  double geotransform[6], inverse_geotransform[6];
  if (dataset->GetGeoTransform(geotransform) != CE_None) {
    throw std::runtime_error("Geo transform not present in the raster");
  }

  if (!GDALInvGeoTransform(geotransform, inverse_geotransform)) {
    throw std::domain_error("Can't invert the geo transform");
  }

  int raster_pixel_coord = static_cast<int>(
      floor(inverse_geotransform[0] + inverse_geotransform[1] * lon +
            inverse_geotransform[2] * lat));
  int raster_line_coord = static_cast<int>(
      floor(inverse_geotransform[3] + inverse_geotransform[4] * lon +
            inverse_geotransform[5] * lat));

  // use the first raster band
  auto band = dataset->GetRasterBand(1);
  double raster_pixel_value[2];
  if (band->RasterIO(GF_Read, raster_pixel_coord, raster_line_coord, 1, 1,
                     raster_pixel_value, 1, 1, GDT_CFloat64, 0, 0) != CE_None) {
    throw std::runtime_error("Raster read error");
  }

  int value = (int)raster_pixel_value[0];

  return value;
}

void gdal_hypsobathymetry_model::set_ocean_geometry(
    const std::filesystem::path& path) {
  if (path.empty() == true) {
    m_ocean_geometry_dataset.reset(nullptr);
    return;
  }
  auto dataset = static_cast<GDALDataset*>(
      GDALOpenEx(path.c_str(), GDAL_OF_VECTOR, nullptr, nullptr, nullptr));
  if (dataset == nullptr) {
    throw std::runtime_error("GDALOpenEx() returned null dataset");
  }
  std::cout << "Set ocean geometry from " << path.string() << " ("
            << dataset->GetDriverName() << ")\n";
  m_ocean_geometry_dataset.reset(dataset);
}
void gdal_hypsobathymetry_model::set_land_geometry(
    const std::filesystem::path& path) {
  if (path.empty() == true) {
    m_land_geometry_dataset.reset(nullptr);
    return;
  }
  auto dataset = static_cast<GDALDataset*>(
      GDALOpenEx(path.c_str(), GDAL_OF_VECTOR, nullptr, nullptr, nullptr));
  if (dataset == nullptr) {
    throw std::runtime_error("GDALOpenEx() returned null dataset");
  }
  std::cout << "Set land geometry from " << path.string() << " ("
            << dataset->GetDriverName() << ")\n";
  m_land_geometry_dataset.reset(dataset);
}
void gdal_hypsobathymetry_model::set_ocean_values_raster(
    const std::filesystem::path& path) {
  if (path.empty() == true) {
    m_ocean_values_dataset.reset();
    return;
  }
  auto dataset = static_cast<GDALDataset*>(GDALOpen(path.c_str(), GA_ReadOnly));
  if (dataset == nullptr) {
    throw std::runtime_error("GDALOpen() returned null dataset");
  }
  if (m_land_values_dataset != nullptr) {
    auto other_dataset = m_land_values_dataset;
    if (other_dataset->GetFileList()[0] == path.string()) {
      std::cout << "Reusing land values raster for ocean: " << path.string()
                << "\n";
      m_ocean_values_dataset = other_dataset;
      return;
    }
  }
  std::cout << "Set ocean values from " << path.string() << " ("
            << dataset->GetDriverName() << ")\n";
  m_ocean_values_dataset.reset(dataset);
}

void gdal_hypsobathymetry_model::set_land_values_raster(
    const std::filesystem::path& path) {
  if (path.empty() == true) {
    m_ocean_values_dataset.reset();
    return;
  }
  auto dataset = static_cast<GDALDataset*>(GDALOpen(path.c_str(), GA_ReadOnly));
  if (dataset == nullptr) {
    throw std::runtime_error("GDALOpen() returned null dataset");
  }
  if (m_ocean_values_dataset != nullptr) {
    auto other_dataset = m_ocean_values_dataset;
    if (other_dataset->GetFileList()[0] == path.string()) {
      std::cout << "Reusing ocean values raster for land: " << path.string()
                << "\n";
      m_land_values_dataset = other_dataset;
      return;
    }
  }
  std::cout << "Set land values raster from " << path.string() << " ("
            << dataset->GetDriverName() << ")\n";
  m_land_values_dataset.reset(dataset);
}

gdal_hypsobathymetry_model::gdal_hypsobathymetry_model() { GDALAllRegister(); }

Визуализация контуров суши/береговых линий⌗

Воспользуемся графическим фреймворком QGraphicsScene.

Унаследуем от QGraphicsItemGroup класс, который при инстанциации конвертирует заданный векторный датасет участков суши в одну большую группу QGraphicsPolygonItem, а также экспозит простой метод для установки толщины линий, цвета контура и штриховки всех многоугольников коллекции.


#include <QBrush>
#include <QGraphicsItem>
#include <QGraphicsItemGroup>
#include <optional>

class GdalCoastlineGraphicsItem : public QGraphicsItemGroup {
 private:
 public:
  GdalCoastlineGraphicsItem(const char* coastline_shp);

  void setPenBrushForAllPolygons(std::optional<qreal> penwidth,
                                 std::optional<QColor> pencolor = std::nullopt,
                                 std::optional<QBrush> brush = std::nullopt);
};

Особый интерес представляют статические функции трансформации вершин геометрий датасета (зеркального отражения точки по вертикали, т.к. в Qt графическая ось Y растёт вниз, а в датасете использованы географические координаты -90…+90) при создании объектов QPolygonF. В системе GDAL/GEOS геометрия типа OGRMultiPolygon содержит несколько OGRPolygon, причём объект OGRPolygon может описывать не один, а несколько многоугольников сразу: обязательный exterior ring (многоугольник) и опциональные interior rings (многоугольники) для внутренних “дырок”. Геометрии можно итерировать по точкам.

static void transform_qpoint(QPointF& point) {
  point.setY(90.f - point.y());
  point.setX(point.x() * 60);  // 1 minute -> 1 pixel scale
  point.setY(point.y() * 60);
}

static QList<QPolygonF> convert_polygon(const OGRPolygon& ogrpoly) {
  QPolygonF qextring;
  auto extring = ogrpoly.getExteriorRing();
  for (auto& pt : extring) {
    QPointF qpt{pt.getX(), pt.getY()};
    transform_qpoint(qpt);
    qextring.push_back(qpt);
  }

  QList<QPolygonF> qintrings;
  if (ogrpoly.getNumInteriorRings() > 0) {
    qDebug("Poly has %d interior rings:", ogrpoly.getNumInteriorRings());
    for (int i = 0; i < ogrpoly.getNumInteriorRings(); ++i) {
      auto intring = ogrpoly.getInteriorRing(i);
      qDebug("\tRing %d (%d points)", i, intring->getNumPoints());
      QPolygonF qintring;
      for (auto& pt : intring) {
        QPointF qpt{pt.getX(), pt.getY()};
        transform_qpoint(qpt);
        qintring.push_back(qpt);
      }
      qintrings.push_back(qintring);
    }
  }
  qintrings.push_front(qextring);
  return qintrings;
}

static QList<QPolygonF> convert_multipolygon(
    const OGRMultiPolygon& ogrmultipoly) {
  QList<QPolygonF> qpolys;
  for (const OGRPolygon* poly : ogrmultipoly) {
    qpolys.append(convert_polygon(*poly));
  }
  return qpolys;
}

Вот полная имплементация:


#include "GdalCoastlineGraphicsItem.h"

#include <gdal/ogrsf_frmts.h>

#include <QGraphicsLineItem>
#include <QGraphicsPathItem>
#include <QPen>

static void transform_qpoint(QPointF& point) {
  point.setY(90.f - point.y());
  point.setX(point.x() * 60);  // 1 minute -> 1 pixel scale
  point.setY(point.y() * 60);
}

static QList<QPolygonF> convert_polygon(const OGRPolygon& ogrpoly) {
  QPolygonF qextring;
  auto extring = ogrpoly.getExteriorRing();
  for (auto& pt : extring) {
    QPointF qpt{pt.getX(), pt.getY()};
    transform_qpoint(qpt);
    qextring.push_back(qpt);
  }

  QList<QPolygonF> qintrings;
  if (ogrpoly.getNumInteriorRings() > 0) {
    qDebug("Poly has %d interior rings:", ogrpoly.getNumInteriorRings());
    for (int i = 0; i < ogrpoly.getNumInteriorRings(); ++i) {
      auto intring = ogrpoly.getInteriorRing(i);
      qDebug("\tRing %d (%d points)", i, intring->getNumPoints());
      QPolygonF qintring;
      for (auto& pt : intring) {
        QPointF qpt{pt.getX(), pt.getY()};
        transform_qpoint(qpt);
        qintring.push_back(qpt);
      }
      qintrings.push_back(qintring);
    }
  }

  qintrings.push_front(qextring);
  return qintrings;
}

static QList<QPolygonF> convert_multipolygon(
    const OGRMultiPolygon& ogrmultipoly) {
  QList<QPolygonF> qpolys;
  for (const OGRPolygon* poly : ogrmultipoly) {
    qpolys.append(convert_polygon(*poly));
  }
  return qpolys;
}

GdalCoastlineGraphicsItem::GdalCoastlineGraphicsItem(
    const char* coastline_shp) {
  GDALAllRegister();

  auto dataset = static_cast<GDALDataset*>(
      GDALOpenEx(coastline_shp, GDAL_OF_VECTOR, nullptr, nullptr, nullptr));
  if (dataset == nullptr) {
    throw std::runtime_error("GDALOpenEx() returned null dataset");
  }

  QList<QPolygonF> total_qpolys;
  auto layer = dataset->GetLayer(0);
  for (auto& feature : layer) {
    auto multipolygon = feature->GetGeometryRef()->toMultiPolygon();
    if (multipolygon == nullptr) {
      qDebug("Geometry for FID %lld is not MultiPolygon", feature->GetFID());
      break;
    }
    auto qpolys = convert_multipolygon(*multipolygon);
    qDebug(
        "Constructed list of total %d polys from a multypolygon for FID %lld",
        qpolys.size(), feature->GetFID());
    total_qpolys.append(qpolys);
  }

  // create the items
  for (auto& qpoly : total_qpolys) {
    auto qpolyitem = new QGraphicsPolygonItem{qpoly};
    this->addToGroup(qpolyitem);
  }

  delete dataset;
}

void GdalCoastlineGraphicsItem::setPenBrushForAllPolygons(
    std::optional<qreal> penwidth, std::optional<QColor> pencolor,
    std::optional<QBrush> brush) {
  for (auto item : childItems()) {
    auto shapeitem = dynamic_cast<QAbstractGraphicsShapeItem*>(item);
    if (shapeitem) {
      QPen pen = shapeitem->pen();
      if (penwidth.has_value()) {
        pen.setWidthF(penwidth.value());
      }
      if (pencolor.has_value()) {
        pen.setColor(pencolor.value());
      }
      shapeitem->setPen(pen);
      if (brush.has_value()) {
        shapeitem->setBrush(brush.value());
      }
    }
  }
}

Тестирование модели⌗

Для удобства просмотра возьмём QGraphicsView и добавим в него испускание сигнала при клике по холсту (трансформируем мышиные координаты в широту/долготу), а также поддержку zoom-in/zoom-out колёсиком мышки с автокорректировкой толщины контура GdalCoastlineGraphicsItem. Весь остальной функционал вроде скроллинга получается автоматически при подключении объекта-сцены благодаря фреймворку Qt Graphics.


#include <QGraphicsView>

class ZoomableGraphicsView : public QGraphicsView {
  Q_OBJECT
 public:
  ZoomableGraphicsView(QWidget* parent);

 signals:
  void lat_lon_clicked(double lat, double lon);

 protected:
  virtual void wheelEvent(QWheelEvent* event) override;
  virtual void mousePressEvent(QMouseEvent* event) override;
};

Полная реализация просмотрщика сцены:


#include "ZoomableGraphicsView.h"

#include <QScrollBar>
#include <QWheelEvent>
#include <QtMath>

#include "GdalCoastlineGraphicsItem.h"

void ZoomableGraphicsView::wheelEvent(QWheelEvent* event) {
  static float total_zoom = 1.f;
  static const qreal ZOOM_STEP = .5;
  if (event->modifiers() & Qt::ControlModifier) {
    // zoom
    const ViewportAnchor anchor = transformationAnchor();
    setTransformationAnchor(QGraphicsView::AnchorUnderMouse);
    int angle = event->angleDelta().y();
    qreal factor;
    if (angle > 0) {
      factor = 1. + ZOOM_STEP;
    } else {
      factor = 1. - ZOOM_STEP;
    }
    scale(factor, factor);
    total_zoom *= factor;

    // rescale outlines of certain graphics objects in the scene
    auto scn = scene();
    if (scn != nullptr) {
      for (auto item : scn->items()) {
        if (auto coastlineitem =
                dynamic_cast<GdalCoastlineGraphicsItem*>(item)) {
          coastlineitem->setPenBrushForAllPolygons({1.f / total_zoom});
        }
      }
    }
    qDebug("Zoom: total zoom %f", total_zoom);

    setTransformationAnchor(anchor);
  } else {
    QGraphicsView::wheelEvent(event);
  }
}

ZoomableGraphicsView::ZoomableGraphicsView(QWidget* parent)
    : QGraphicsView(parent) {}

void ZoomableGraphicsView::mousePressEvent(QMouseEvent* event) {
  QPointF mousepos = event->localPos();
  auto scenepos = mapToScene(event->x(), event->y());
  double lat = 90. - (scenepos.y() / 60);
  double lon = scenepos.x() / 60.;

  emit lat_lon_clicked(lat, lon);

  QGraphicsView::mousePressEvent(event);
}

В главном окне создаём графику береговых линий, устанавливаем цвет штриховки.

auto coastline = new GdalCoastlineGraphicsItem{land_shp};
coastline->setPenBrushForAllPolygons(
    std::nullopt, QColor{0, 128, 0},
    QBrush{QColor{15, 15, 15}});  // landmass bg color

Добавляем её в наш просмотрщик, подписываемся на его сигнал о клике по холсту (т.е. широте/долготе).

auto scene = new QGraphicsScene;
scene->addItem(coastline);
// oceans background color
ui->graphicsView->setBackgroundBrush(QBrush{QColor{27, 38, 49}});
ui->graphicsView->setScene(scene);
QObject::connect(ui->graphicsView, &ZoomableGraphicsView::lat_lon_clicked,
                 this, &MainWindow::lat_lon_requested);

В конструкторе окна мы инициализировали нашу гипсобатиметрическую модель – пришло время её опросить по кликам мышки:

void MainWindow::lat_lon_requested(double lat, double lon) {
  auto on_land = model->is_land(lat, lon);
  ui->label->setText(QString::asprintf(
      "LatLon: %f, %f, %s, %s %d meters", lat, lon, on_land ? "land" : "ocean",
      on_land ? "altitude" : "depth", model->value(lat, lon)));
}

Берега Южной Америки:

Где-то в Тихом океане:

Апеннинский полуостров: