国土地理院の標高タイルだけで、住所の傾斜・谷/尾根・川までの距離を計算する
住所を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つだけです。
-
シグネチャ確認(
89 50 4E 47 0D 0A 1A 0A) -
チャンクを読む。
length(4) + type(4) + data + crc(4)の繰り返し。
必要なのはIHDR(幅高さ・bit深度・カラータイプ)とIDAT(本体) -
IDATを連結してinflateSync。zlib は標準ライブラリにある -
フィルタ復元。各行の先頭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ピクセルの距離は緯度で変わる。
metersPerPixel(lat, zoom)を必ず通す -
標高タイルは符号付き24bitで読む。
(128,0,0)は無効値。
0m として扱わずnullにして確度に反映する - 表示用と判定用でデコーダを分けてよい。 deck.gl は線形式しか受けない
-
入力の形が固定なら、PNGデコーダは自分で書ける(zlib + フィルタ復元)。
汎用PNG相手には勧めない -
傾斜はスケール付きの量。 基線長を決めずに「◯度」と言ってはいけない。
本物の急傾斜は基線を変えても安定し、ノイズは崩れる -
谷/尾根は「周囲の平均との差」と「4近傍のうち何方向が高いか」で見る。
標高が低いことと周囲より低いことは別 -
ラスタの色で距離を測るのは成立する。 ただし配色変更で壊れ、
暗渠は写らず、池と川を区別できない。だから帯に丸めて確度を下げる - フォールバックはHTTPステータスではなく「欲しい値が読めたか」で決める
一番の学びは5番でした。ライブラリを使わずに自分で計算すると、
出てきた数字を疑う責任も自分に来ます。 ライブラリなら
「基線長」のようなパラメータが引数として表に出てくるので、
決めないまま使うことができません。自分で書くと、暗黙の定数として
埋まったまま本番に出ます。
計算式が合っていても、測り方(スケール)が用途に合っていなければ
出力は嘘になります。 そしてその嘘は、例外も型エラーも出しません。