0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

国土地理院の標高タイルだけで、住所の傾斜・谷/尾根・川までの距離を計算する(そして「傾斜」を短い基線で測ると平地が坂になる)

0
Posted at

国土地理院の標高タイルだけで、住所の傾斜・谷/尾根・川までの距離を計算する

住所を1つもらって、その地点の傾斜・谷か尾根か・川までの距離を返す処理を
書いて運用しています。使うのは国土地理院の標高タイルと地図タイルだけです。

  • APIキー不要、無料
  • PostGIS も GDAL も使わない
  • ベクタデータ(河川のライン等)も使わない
  • 依存パッケージを増やさない(PNGデコーダも自分で書いた)

やってみると「思ったより簡単な部分」と「油断すると平気で嘘の値を出す部分」が
はっきり分かれました。特に傾斜は、基線の長さを決めないと意味がありません。
最後にその実測を載せます。自分のコードのバグも見つかりました。


1. タイル座標と「1ピクセルが何メートルか」

標高タイルは https://cyberjapandata.gsi.go.jp/xyz/{layer}/{z}/{x}/{y}.png で取れます。
まず緯度経度をタイル内のピクセルに落とします。Webメルカトルの標準的な式です。

const EARTH_RADIUS_METERS = 6_378_137;
const TILE_SIZE = 256;

export function latLngToWorldPixel(latitude: number, longitude: number, zoom: number) {
  const scale = TILE_SIZE * 2 ** zoom;
  const x = ((longitude + 180) / 360) * scale;
  const latRadians = (latitude * Math.PI) / 180;
  const y =
    ((1 - Math.log(Math.tan(latRadians) + 1 / Math.cos(latRadians)) / Math.PI) / 2) * scale;
  return { x, y };
}

worldPixel を 256 で割れば tileX/tileY、余りがタイル内座標です。
余りは負にならないようにしてください(西経・北半球外でズレます)。

function floorMod(value: number, divisor: number) {
  const remainder = value % divisor;
  return remainder >= 0 ? remainder : remainder + divisor;
}

そしてここが重要です。1ピクセルが表す距離は緯度で変わります。

export function metersPerPixel(latitude: number, zoom: number) {
  return (
    (Math.cos((latitude * Math.PI) / 180) * 2 * Math.PI * EARTH_RADIUS_METERS) /
    (TILE_SIZE * 2 ** zoom)
  );
}

z=15、緯度34.7度(大阪あたり)では 1px ≒ 3.93m。
札幌(43度)なら約 3.5m、那覇(26度)なら約 4.3m です。
傾斜を出すときにこれを定数にすると、北と南で答えが変わります。


2. RGB から標高へ(負の標高と無効ピクセル)

国土地理院の標高タイルは、標高をRGBに詰めたPNGです。公式はこれだけです。

elevation_m = (R * 65536 + G * 256 + B) * 0.01

素直に実装すると海中でおかしくなります。**値は符号付き24bit(2の補数)**です。

function readElevationPixel(image: DecodedPng, pixelX: number, pixelY: number) {
  const rgb = readRgbPixel(image, pixelX, pixelY);
  if (!rgb) return null;
  const [red, green, blue] = rgb;

  // 無効値。データが無い場所を 0m として扱ってはいけない
  if (red === 128 && green === 0 && blue === 0) return null;

  const rawValue = red * 65536 + green * 256 + blue;
  const signedValue = rawValue < 2 ** 23 ? rawValue : rawValue - 2 ** 24;  // ★
  return signedValue * 0.01;
}

2点あります。

  • (128, 0, 0) は無効値。 これを 0m として扱うと、データが無い海や
    範囲外を「標高0mの平地」と誤認します。null にして、
    「有効サンプルが何点取れたか」を確度の判定に使うべきです
  • 符号付きとして読む。 そうしないと標高 -1m の土地が
    約 167,772m になります

余談ですが、同じプロジェクト内に線形デコーダも別に持っています。

/** deck.gl TerrainLayer 用の elevationDecoder 係数。*/
export const GSI_ELEVATION_DECODER = {
  rScaler: 65536 / 100, gScaler: 256 / 100, bScaler: 1 / 100, offset: 0,
};

deck.gl は R*rScaler + G*gScaler + B*bScaler + offset という線形式しか
受け取れない
ので、符号反転を表現できません。そちらは「無効ピクセルと負の標高は
タイル生成の段階で 0 にクランプする」という前処理とセットにしています。

同じデータでも、用途(表示用 / 判定用)でデコーダを分けるのは
悪いことではありません。分けるなら理由をコメントに書いておいてください。


3. PNGデコーダを自分で書いた(203行)

Node には canvas がありません。node-canvas や sharp を入れれば済みますが、
ネイティブ依存が増えてデプロイが重くなります。標高タイルは
8bit RGB/RGBA の小さなPNGしか来ないので、必要な範囲だけ自分で書きました。

import { inflateSync } from "node:zlib";

export function decodePng(bytes: Uint8Array): DecodedPng {
  assertSignature(bytes);
  // IHDR / PLTE / tRNS / IDAT を読む → inflate → フィルタ復元
}

やることは4つだけです。

  1. シグネチャ確認(89 50 4E 47 0D 0A 1A 0A)
  2. チャンクを読む。length(4) + type(4) + data + crc(4) の繰り返し。
    必要なのは IHDR(幅高さ・bit深度・カラータイプ)と IDAT(本体)
  3. IDAT を連結して inflateSync。zlib は標準ライブラリにある
  4. フィルタ復元。各行の先頭1バイトがフィルタ種別(0〜4)で、
    上の行・左のピクセルから復元する

4番だけ少し手間ですが、仕様は短く、テストもしやすいです
(同じPNGを別のデコーダと突き合わせれば済む)。

ネイティブ依存を避けたいときは選択肢になります。
ただし汎用PNGを相手にするなら(16bit、インターレース、
巨大画像)やめたほうがいいです。入力の形が固定だから成立する近道です。


4. 傾斜の計算 — そして基線長の罠

標高が読めれば傾斜は出ます。4近傍の中央差分です。

const dzdx = (east - west) / (4 * metersPerPixelValue);
const dzdy = (south - north) / (4 * metersPerPixelValue);
const slopeDegrees = (Math.atan(Math.sqrt(dzdx ** 2 + dzdy ** 2)) * 180) / Math.PI;

4 * は「±1ピクセルではなく ±2ピクセル離れた点を使っているので、
東西の距離は 2px × 2 = 4px」という意味です。ここは実装に合わせてください。

問題はこの「±2ピクセル」です。

z=15 では 1px ≒ 3.93m なので、±2px は基線わずか15.7m。
この距離で標高差を測ると、地形ではなく construction(建物・護岸・盛土)や
DEMのノイズを拾います。

実際に、基線の長さを変えて同じ3地点を測りました。

地点 ±2px(15.7m) ±4px(31.4m) ±8px(62.8m) ±16px(125.6m)
六甲山上 32.6度 32.6度 35.5度 38.0度
大阪の低地 7.6度 5.6度 (無効値) 0.4度
淀川左岸 9.8度 9.4度 8.7度 1.9度

読み取れることがはっきりしています。

  • 本物の急傾斜(六甲山上)は基線を変えても崩れない。 32.6 → 38.0度
  • 平地は基線を伸ばすと崩壊する。 大阪の低地は 7.6度 → 0.4度。
    つまり短基線の 7.6度は地形ではない
  • 中間のケースもある。 淀川左岸は 62.8m で 8.7度、125.6m で 1.9度。
    これは河川の堤防を拾っていて、どちらも「嘘」ではありません。
    「徒歩の負荷」を知りたいなら堤防込みの短基線が正しく、
    「土地の傾き」を知りたいなら長基線が正しい

これは自分のコードのバグでした

自分の実装は ±2px(15.7m)を使っていて、判定はこうです。

function classifySlope(slopeDegrees: number | null): TerrainSlopeClass {
  if (slopeDegrees === null || slopeDegrees < 2) return "flat";
  if (slopeDegrees < 5) return "gentle";
  if (slopeDegrees < 10) return "moderate";
  return "steep";
}

大阪の低地(真の傾斜 0.4度)が moderate に分類されます。
そしてレポートには「移動負荷:高」の側の文が出ます。
まったいらな土地に「坂の負荷があります」と書いていたわけです。

同じ地点の局所起伏(半径約94mでの標高差)は 4.2m しかないので、
データの中に矛盾は最初からありました。見ていなかっただけです。

教訓

「傾斜」はスカラーではなくスケール付きの量です。
基線長を決めずに「傾斜◯度」と言うのは、単位を書かずに「3」と言うのと同じです。

実装するなら、

  • 基線長を定数として明示し、名前に入れる
    (slopeDegreesOver30m のように)
  • 用途ごとに複数の基線で出す。 徒歩負荷なら短め、地形分類なら長め
  • 矛盾検査を入れる。 「局所起伏が小さいのに傾斜が大きい」は
    ほぼノイズなので、そこで確度を落とすか flat に寄せる
  • 短基線しか使わないなら、閾値を上げる(2/5/10度は長基線向けの感覚値)

5. 谷か尾根かを判定する

傾斜の「向き」だけでは谷と尾根が区別できません(どちらも斜面)。
周囲より低いか高いかと近傍の落ち方を見ます。

// 半径24px・4px刻みでサンプリングした平均との差
const relativeElevationMeters = centerElevation - meanMeters;
// サンプル内の最大 - 最小
const localReliefMeters = maxMeters - minMeters;
// 4近傍のうち、自分より高い/低い点の数
const lowerNeighbors  = [n, s, e, w].filter((v) => centerElevation < v).length;
const higherNeighbors = [n, s, e, w].filter((v) => centerElevation > v).length;

const valleyTendency =
  (relativeElevationMeters <= -4 && lowerNeighbors >= 3) ||
  (localReliefMeters >= 12 && lowerNeighbors >= 3);
const ridgeTendency =
  (relativeElevationMeters >= 4 && higherNeighbors >= 3) ||
  (localReliefMeters >= 12 && higherNeighbors >= 3);
  • 相対標高(周囲の平均との差)が谷/尾根の主要な手がかり
  • lowerNeighbors >= 3、つまり4方向のうち3方向が自分より高ければ
    「窪んでいる」。1方向だけなら単なる斜面
  • 局所起伏を併用して、平坦地での微小な差で谷判定が出ないようにする

低地判定は別条件です。標高が低いことと、周囲より低いことは違います。

const lowlandTendency =
  relativeElevationMeters <= -3 ||
  (centerElevation <= 10 && relativeElevationMeters <= -1.5);

「標高10m以下」だけで低地と決めないのがポイントです。
海沿いの台地は標高が低くても周囲より高いので、浸水リスクの意味が変わります。


6. ベクタデータなしで「川までの距離」を測る

河川のラインデータを取ってきて最近傍距離を計算するのが正攻法です。
やりたくなかったので、地図タイルの色で測りました。

淡色地図(pale)を z=16 で取って、水色っぽいピクセルを同心円状に探します。

for (let offsetY = -WATER_SEARCH_RADIUS; offsetY <= WATER_SEARCH_RADIUS; offsetY += WATER_SEARCH_STEP) {
  for (let offsetX = -WATER_SEARCH_RADIUS; offsetX <= WATER_SEARCH_RADIUS; offsetX += WATER_SEARCH_STEP) {
    const rgb = await readRgbAtWorldPixel(center.worldPixelX + offsetX, center.worldPixelY + offsetY, "pale", MAP_ZOOM, fetchBytes);
    if (!rgb || !isWaterLikePixel(rgb[0], rgb[1], rgb[2])) continue;

    const distance = Math.sqrt(offsetX ** 2 + offsetY ** 2) * metersPerPixelValue;
    if (nearestDistanceMeters === null || distance < nearestDistanceMeters) {
      nearestDistanceMeters = distance;
    }
  }
}

判定は閾値の組み合わせだけです。

function isWaterLikePixel(red: number, green: number, blue: number) {
  return blue >= 150 && blue - red >= 35 && blue - green >= 10 && red <= 180 && green <= 220;
}

blue >= 150 だけでは空色の記号や凡例に反応するので、
**「青が赤より十分強く、緑より少し強い」**という差分条件を足しています。

得られた距離は帯に丸めます。ピクセル由来の値を1m単位で出すと精度を偽装するので、
最初から段階にしてしまうほうが誠実です。

function classifyWaterBand(distanceMeters: number): TerrainWaterBand {
  if (distanceMeters <= 60) return "adjacent";
  if (distanceMeters <= 180) return "near";
  if (distanceMeters <= 320) return "moderate";
  return "distant";
}

この方法の限界(正直に)

  • 地図のデザイン変更で壊れます。 タイルの配色は仕様ではありません
  • 暗渠・地下河川は写りません。 見えないものは測れない
  • プールや池も「水域」になります。 河川と区別できない
  • 探索半径(96px × z=16)の外に水域があれば null。
    「川が遠い」ではなく「分からない」として扱うべき

なので確度を低めに固定しています。

confidence: roundedDistance <= 320 ? "medium" : "low",

近くで見つかったときだけ medium。 遠い結果は「見つからなかった」に近いので
low にしています。ラスタから測った値を high で出してはいけません。


7. DEMレイヤのフォールバックは「無効ピクセル」で決める

標高タイルは1種類ではありません。5mメッシュは市街地中心で、
山間部では404だったり無効値だったりします。

const DEM_LAYERS = ["dem5a_png", "dem5b_png", "dem10b_png"] as const;

async function chooseDemLayer(tileX, tileY, pixelX, pixelY, fetchBytes) {
  for (const layerId of DEM_LAYERS) {
    const tile = await loadPngTile(buildGsiTileUrl(layerId, DEM_ZOOM, tileX, tileY), fetchBytes);
    const elevation = readElevationPixel(tile, pixelX, pixelY);
    if (elevation !== null) return { layerId };   // ★ 値が読めた層を採用
  }
  return { layerId: "dem10b_png" as const };
}

HTTPステータスで判定しないのが要点です。 タイルは 200 で返ってくるのに
中身が無効値 (128,0,0) というケースがあります。
「実際に欲しい1点が読めたか」で採用を決めるのが確実です。

そして採用したレイヤーを確度に反映します(10mメッシュなら確度を下げる)。
どのデータで答えたかを出力に残しておくと、後から検証できます。


8. 実測(3地点)

最後に、実際に動かした結果です。標高タイルは無料なので誰でも再現できます。

六甲山上   1918ms  標高=855.8m  傾斜=32.6度(steep)     局所起伏=146.4m  水域=194m(moderate)  確度=high
大阪の低地   766ms  標高=3.6m    傾斜=7.6度(moderate)   局所起伏=4.2m    水域=24m(adjacent)   確度=high
淀川左岸     656ms  標高=3.7m    傾斜=9.8度(moderate)   局所起伏=7.8m    水域=40m(adjacent)   確度=high

タイルをキャッシュしていない初回で1〜2秒。標高の絶対値と水域距離は妥当です。
傾斜だけが前述のとおり信用できません。 そして確度は high と出ています。
確度の計算に「傾斜の妥当性」が入っていなかったわけで、ここも直す対象です。


まとめ

  1. 1ピクセルの距離は緯度で変わる。 metersPerPixel(lat, zoom) を必ず通す
  2. 標高タイルは符号付き24bitで読む。 (128,0,0) は無効値。
    0m として扱わず null にして確度に反映する
  3. 表示用と判定用でデコーダを分けてよい。 deck.gl は線形式しか受けない
  4. 入力の形が固定なら、PNGデコーダは自分で書ける(zlib + フィルタ復元)。
    汎用PNG相手には勧めない
  5. 傾斜はスケール付きの量。 基線長を決めずに「◯度」と言ってはいけない。
    本物の急傾斜は基線を変えても安定し、ノイズは崩れる
  6. 谷/尾根は「周囲の平均との差」と「4近傍のうち何方向が高いか」で見る。
    標高が低いことと周囲より低いことは別
  7. ラスタの色で距離を測るのは成立する。 ただし配色変更で壊れ、
    暗渠は写らず、池と川を区別できない。だから帯に丸めて確度を下げる
  8. フォールバックはHTTPステータスではなく「欲しい値が読めたか」で決める

一番の学びは5番でした。ライブラリを使わずに自分で計算すると、
出てきた数字を疑う責任も自分に来ます。
ライブラリなら
「基線長」のようなパラメータが引数として表に出てくるので、
決めないまま使うことができません。自分で書くと、暗黙の定数として
埋まったまま本番に出ます。

計算式が合っていても、測り方(スケール)が用途に合っていなければ
出力は嘘になります。
そしてその嘘は、例外も型エラーも出しません。

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?