QGISを参考にgeotiff.jsである座標のピクセル値を取得する機能を実装

ブラウザでgeotiffを読み込んである地点の値を抽出する必要があり
geotiff.jsを使ってラスターサンプリング機能を実装しました。
実装するにあたりQGISのソースコードを参考にしました。
QGISからWebへの機能の移植の例を共有したいと思います。

この関数がラスターサンプリング機能です。

import type { GeoTIFF, TypedArray } from "geotiff";

/**
 * 指定された経緯度の値をGeoTIFFデータから取得
 * @param latitude
 * @param longitude
 * @param bandIndex
 * @param nodataValue
 * @returns {number} - 値
 */
const getValue = async (
  geoTIFF: GeoTIFF,
  latitude: number,
  longitude: number,
  bandIndex: number,
  nodataValue: number,
): Promise<number> => {
  const image = await geoTIFF.getImage();
  const data = await image.readRasters();
  const [originX, originY] = image.getOrigin();
  const [resolutionX, resolutionY] = image.getResolution();
  const x = Math.round((longitude - originX) / resolutionX);
  const y = Math.round((latitude - originY) / resolutionY);
  const index = y * data.width + x;
  const band = data[bandIndex] as TypedArray;
  const value = band[index];
  if (value === nodataValue) {
    return 0;
  }
  return value;
};

実装方法

実際にQGISのソースコードを読みどこで目的の処理が実装されているか
特定した上でTypeScriptに実装しました。
処理を特定した過程を下記にまとめます。

以下は2023/04/18時点のmasterブランチの話です

これ以下のコードはQGISの引用となっております。

機能の名前を探す

プロセシングツールボックス->ベクタレイヤにラスタ値を付加
マウスオーバーするとアルゴリズムID(native:rastersampling)がわかります。

QGISソースコードを読む

QGISのリポジトリにアクセスして
rastersamplingで検索すると
qgsalgorithmrastersampling.cppというC++ファイルがヒットします。

QgsRasterSamplingAlgorithm::processAlgorithmを読み進めます。

const double value = mDataProvider->sample( point, band, &ok );

上記のコードが見られますこれが前後の処理や引数からこの処理でラスターをサンプリングしていることがわかります。
sampleの実装がまだ分からないので、mDataProviderの宣言をヘッダーファイルで読みます。

ヘッダーファイルを探す

キーボードのtキーを押下するとファイル検索ができるので
qgsalgorithmrastersampling.hを開きます。

std::unique_ptr< QgsRasterDataProvider > mDataProvider;

mDataProviderはQgsRasterDataProviderクラスだとわかりました。
従ってQgsRasterDataProviderのsampleメソッドを読みます。

再度tキーを押下してQgsRasterDataProviderがファイル名に含まれるファイルを探します。
qgsrasterdataprovider.cppを読みます。

QgsRasterDataProviderクラスを読む

QgsRasterDataProvider::sample

const auto res = identify( point, QgsRaster::IdentifyFormatValue, boundingBox, width, height, dpi );
const QVariant value = res.results().value( band );

このコードとsampleメソッドがvalueを浮動小数点数に変換した値を返すことから
identifyを読むと遂に座標から画像の行と列を取得する処理を見つけることができました。

// Calculate the row / column where the point falls
const double xres = ( finalExtent.width() ) / width;
const double yres = ( finalExtent.height() ) / height;

const int col = static_cast< int >( std::floor( ( point.x() - finalExtent.xMinimum() ) / xres ) );
const int row = static_cast< int >( std::floor( ( finalExtent.yMaximum() - point.y() ) / yres ) );

まとめ

今回はQGISの処理の特定は比較的容易であるということを例示してみました。
特定するポイントをまとめました。

• QGISのアルゴリズムはプロセシングバーのアルゴリズムIDから探せる
• メソッドや関数の引数に注目する
• クラスや型を足がかりに探す
• クラスや型を知るために宣言を探す