はじめに
USGSのSRTMデータから、RGB標高タイルを作った記録です。
動機
MapLibreでも3D Terrain 表示がサポートされ、今やMapbox GL JS でもMapLibre GL JSでも、3D表現やHillshadeの表現が可能になりました。そして、これらの機能の利用のためにはウェブ地図向けの標高データ、つまりRGB標高タイルが必要です。
世の中には自由に使える標高データ、例えばUSGSからリリースされているパブリックドメインの標高データなどがありますが、自分で自由に使うことができるRGB 標高タイル形式のデータはあまりないので、それなら自分で作ろうと思いました。これまでにZL8までのタイルは地球地図から作ったところですが、30メートル解像度（1 arc Second）のSRTM DEMを使ってZL12までのRGB標高タイルを作りたいと思います。
また、EsriさんもArcGISシリーズ用に等高線のデータやHillshadeのラスタデータを提供していますが、こちらもEsriさんの製品のユーザーの方向けですので、自分でRGB標高タイルを作っておくことには一定の価値があると思います。
私の環境
- データ処理
- nodejs: v16.16.0
- npm: 8.11.0
- GDAL: 3.0.4
- Platform: Ubuntu 20.04.4 LTS (built on Docker for windows)
- データダウンロード
- Windows 10 のPCでやりました。（Google Chrome）
【参考】RHELの環境だと、gdalをインストールするところでつまずいたので、Dockerfileを作ってWindowsPCでやった方が楽だと思います。
https://github.com/unvt/rgbify-node
SRTMについて
Shuttle Rader Topography Missionの略で、スペースシャトルのレーダーで作ったグローバルな標高データです。ただし、グローバルといっても高緯度のデータはありません(だいたい北緯60度～南緯56度)。
グローバルな30メートル（１秒）解像度のデータが公開されています。データにはいくつかバージョンがありますが、バージョン３で、データの欠損地域をASTER GDEMやUSGSのGMTEDで補完しているデータを使いました。データはUSGSが公開していて、EarthExplorerなどからダウンロード出来ます。
方針の検討
解像度とズームレベル
ソースデータが1秒（30メートル）解像度とういことを考えれば（下の表）、基本的にはZL12までのものを作る方針とし、データサイズに余裕があればZL13を目指すということでよいと思います。
（参考） ズームレベルと全世界のピクセル数（256ピクセルタイルの場合）
| ズームレベル | タイル数 | 全世界のピクセル数 | 東西方向・1ピクセルの長さ（赤道） | 東西方向・1ピクセルの長さ（北緯60度） | 備考 |
|---|---|---|---|---|---|
| 0 | 1 by 1 | 256 by 256 | 156,259 m | 78,125 m | |
| 1 | 2 by 2 | 512 by 512 | 78,125 m | 39,063 m | |
| ... | ... | ... | ... | ... | |
| 10 | 1,024 by 1,024 | 262,144 by 262,144 | 153 m | 76 m | |
| 11 | 2,048 by 2,048 | 524,288 by 524,288 | 76 m | 38 m | 30メートル解像度のDEMではここまでは十分 |
| 12 | 4,096 by 4,096 | 1,048,576 by 1,048,576 | 38 m | 19 m | 30 メートル解像度だと高緯度で厳しい |
| 13 | 8,192 by 8,192 | 2,097,152 by 2,097,152 | 19 m | 10 m | ピクセルのサイズが解像度より小さいので結構厳しいが、ZL12より元データを再現できるのでここまで作るというのもあり |
データサイズのあらあらの見積もり
以前、RGB標高タイルは平均100kbと仮定して、世界の30％が陸という様な感じで全体のサイズをざっくり見積ってみたところではZL12までで640GB くらいかなと思っていました。これ以上大きいと一つのサーバーには入らないですし、やはりデータサイズの関係からもまずはZL12を目指します。
※実際にやってみたら、この予想よりデータサイズは小さくなることがわかりました。
作業方針を考える
下の図の緑のボックスが、SRTMのデータの範囲です（全部で14,277個）。上にズームレベル6のタイル区画を重ねながら考えてみます。（投影が緯度経度なので、ウェブメルカトルタイルは高緯度で小さく見えます。）
ズームレベル6のタイル区画ごとに、その領域をカバーするSRTMファイルをマージして、それをタイル化していこうと思います。 この作戦は、これまでのベクトルタイル作成の経験でも効果がありました。UNVTの作業ではよく使う基本戦術ですね。
- SRTMは1度×1度ごとのデータになっていて（3601×3601ピクセル）、全部で14,277ファイル、300GB以上のサイズがあります。これを全部マージしてタイル化するのは至難の業です（gdal_mergeも動かないだろうと思います。）。どういう区画で分けるかを考える必要があります。
- ZL6のタイル区画で着れば、1つの区画をカバーする1度×1度のタイルは最大で49（7×7）。SRTMの1ファイルは約25MBなので、最大の区画でも1.2GB程度の画像ファイルを処理すればいい。1.2GBならgdal_mergeでも扱った経験がある。→ZL5だとこの4倍なのでちょっと難しい。
- なお、ズームレベル6のタイル区画は4096あるが、SRTMのデータがあるのは932区画（図中赤枠） で、残りの3162にはデータが含まれていない。
（将来へのメモ）
高緯度のデータの存在範囲の境界では、タイルの境目とデータの存在範囲が一致しない。ZL2－8を地球地図データからのタイルにして、ZL9-12をSRTMタイルにするとすると、一部のところではZL9以降でデータの欠損が目立つかもしれない。一番上と一番下は必要に応じて落とすのかなと思います。（まったくのタイル範囲外であれば地球地図データのオーバーズーミングがされるが、SRTMが半分あって半分ないタイルだとデータ欠損がそのまま出る。）
実際の作業
Step 1: SRTM データダウンロード
USGSのページから、SRTMの1秒グローバルデータ（TIFF形式）をダウンロードします。14,277ファイル、合計315GBでした。
Earth Explorerのバルクダウンロードを使うと便利でした。土日はすいすいダウンロードできましたが、平日はゆっくりだったので、ダウンロードに数日を要しました。(土日を含めて約5日 )
やりかたはこちら↓です。
Step 2： nodejs スクリプトを書く方針を考える
データがダウンロードできたので、以下の様なことを考えながらスクリプトを準備したいと思います。
- Z-X-Yの範囲を指定する（Z=6）。
- その範囲をカバーするのに必要な1度×1度タイルをリストする。
- リストした1度×1度のタイルにソースデータ（SRTMデータ）があるかチェックする。
- （その後、マージして、RGBタイル化）
あるZ-X-Yの中のSRTMタイルのリスト化
あるz-x-yタイルをカバーするために必要な1度×1度のタイルの範囲については、以下の絵の様な考え方になると思います。
z-x-yのリストをkeysというアレーに入れるとして、特定のz-x-yにあるべき1度×1度のタイルをリストして、それがソースファイルにあるかを確認するという方法をとれば、特定のz-x-yで使うべきファイルがわかるはずです。
for (const key of keys){
//for (const key of ['6-31-31','6-32-32']){
let [tilez, tilex, tiley] = key.split('-')
tilex = Number(tilex)
tiley = Number(tiley)
tilez = Number(tilez)
const bbox = tilebelt.tileToBBOX([tilex, tiley, tilez])
modulesObj[key] = []
for (x=Math.floor(bbox[0]); x < bbox[2]; x++ ){
m = x.toString(10) // 10 means decimal
if(x < 0) {
m = m.replace("-","")
if(m.length == 1){
m = `00${m}`
} else if (m.length == 2) {
m = `0${m}`
}
m = `W${m}`
} else {
if(m.length == 1){
m = `00${m}`
} else if (m.length == 2) {
m = `0${m}`
}
m = `E${m}`
} // Then, m has proper string
for (y = Math.floor(bbox[1]); y < bbox[3]; y++){
n = y.toString(10)
if(y<0){
n = n.replace("-","")
if(n.length == 1) {
n = `0${n}`
}
n = `S${n}`
} else {
if(n.length == 1) {
n = `0${n}`
}
n = `N${n}`
}
nm = `${n.toLowerCase()}_${m.toLowerCase()}_1arc_v3.tif`
if(srtmFiles.includes(nm)){
//console.log (`${nm}---> yes(${key})`)
modulesObj[key].push(`${srcDir}/${nm}`)
}
}
}
if (Object.keys(modulesObj[key]).length == 0) {
emptyModules.push(key)
}
if (modulesObj[key].length == 0){
delete modulesObj[key]
}
}
Step 2-1: マージ
作業を進めるため、nodejsを作って、gdal_merge.pyをspawnで動かすスクリプトを書いてみました。全世界を1つのファイルにマージするのではなく、932の区画（ZL6のタイル区画でデータが存在するところ）ごとに、区画をカバーできるTIFファイルをマージして作ります
https://github.com/ubukawa/srtm-rgbify/blob/main/index4merge-re.js
途中で止めて再開してもいいように、出力先フォルダ(merge)にファイルがあるとマージをスキップするようにしています。（ただし、マージを途中で止めてしまったファイルは一度消しておかないといけないので注意しました。）
本当は、spawnで出力したものをそのままrasterioコマンドに流し入れて1つのスクリプトで最後まで処理したかったのですが（spawnのexitのあとに、またspawnを入れる・・・）、better-queueのなかで並列処理を上手く回すことができなかったのでマージの作業とRgbifyの作業を分けて行うことにしました。
上記のスクリプトを実行することで、並列で（concurrent）gdal_mergeをすることが出来ます。1つで処理できる完璧なスクリプトの追求に時間をかけるよりは、成果品を早く得るために妥協することにします。（nodejsを上手に使える人がうらやましいなぁと思います・・・。）
2日くらい回して、1/3強は終わっているのですが、どうも時間がかかったりかからなかったり、何がコンバージョンを律速しているのかよくわかりません。外付けHDDからソースの読み込んで、外付けHDDにマージファイルを書き出しているので、ディスクのリードライトの時間が強く影響しているののかもしれません。（conccurent 3くらいでやったときに、3つのプロセスで処理するタイルが多いと時間が非常にかかる印象。）
とはいいながら、約6日間（5日18時間程度） で作業ができました。基本的に並列処理数conccurentは3で（少しだけ2で進めたときもあります）、出来たファイル数は932、合計サイズは540GB でした。作業中の画面は以下の様な感じでした。
...
[6-30-29,6-30-30,6-31-18] in progress
--- 6-31-18 (357/932) starts: 9 src file/files
--- 6-31-18: 9 src file/files merge ends (2022-07-30T22:26:04.412Z --> 2022-07-30T22:27:55.110Z )
[6-30-29,6-30-30,6-31-19] in progress
--- 6-31-19 (358/932) starts: 18 src file/files
--- 6-30-29: 42 src file/files merge ends (2022-07-30T21:25:43.100Z --> 2022-07-30T22:29:53.708Z )
[6-30-30,6-31-19,6-31-20] in progress
--- 6-31-20 (359/932) starts: 23 src file/files
--- 6-31-19: 18 src file/files merge ends (2022-07-30T22:27:55.117Z --> 2022-07-30T22:37:41.421Z )
[6-30-30,6-31-20,6-31-21] in progress
--- 6-31-21 (360/932) starts: 28 src file/files
--- 6-31-20: 23 src file/files merge ends (2022-07-30T22:29:53.782Z --> 2022-07-30T22:43:34.722Z )
[6-30-30,6-31-21,6-31-22] in progress
--- 6-31-22 (361/932) starts: 16 src file/files
--- 6-30-30: 48 src file/files merge ends (2022-07-30T22:18:11.577Z --> 2022-07-30T22:52:03.729Z )
[6-31-21,6-31-22,6-31-23] in progress
--- 6-31-23 (362/932) starts: 28 src file/files
--- 6-31-22: 16 src file/files merge ends (2022-07-30T22:43:34.724Z --> 2022-07-30T22:55:11.647Z )
[6-31-21,6-31-23,6-31-24] in progress
--- 6-31-24 (363/932) starts: 29 src file/files
--- 6-31-21: 28 src file/files merge ends (2022-07-30T22:37:41.422Z --> 2022-07-30T23:17:49.526Z )
[6-31-23,6-31-24,6-31-25] in progress
--- 6-31-25 (364/932) starts: 35 src file/files
--- 6-31-23: 28 src file/files merge ends (2022-07-30T22:52:03.739Z --> 2022-07-31T02:02:10.165Z )
[6-31-24,6-31-25,6-31-26] in progress
--- 6-31-26 (365/932) starts: 30 src file/files
--- 6-31-24: 29 src file/files merge ends (2022-07-30T22:55:11.650Z --> 2022-07-31T03:53:12.044Z )
[6-31-25,6-31-26,6-31-27] in progress
--- 6-31-27 (366/932) starts: 42 src file/files
--- 6-31-25: 35 src file/files merge ends (2022-07-30T23:17:49.549Z --> 2022-07-31T05:04:23.999Z )
[6-31-26,6-31-27,6-31-28] in progress
--- 6-31-28 (367/932) starts: 36 src file/files
--- 6-31-26: 30 src file/files merge ends (2022-07-31T02:02:10.169Z --> 2022-07-31T05:26:30.602Z )
[6-31-27,6-31-28,6-31-29] in progress
...（途中省略）
--- 6-62-41 (919/932) starts: 8 src file/files
--- 6-62-39: 9 src file/files merge ends (2022-08-03T14:38:11.613Z --> 2022-08-03T14:40:09.559Z )
[6-62-40,6-62-41,6-62-42] in progress
--- 6-62-42 (920/932) starts: 2 src file/files
--- 6-62-42: 2 src file/files merge ends (2022-08-03T14:40:09.563Z --> 2022-08-03T14:40:32.324Z )
[6-62-40,6-62-41,6-62-43] in progress
--- 6-62-43 (921/932) starts: 2 src file/files
--- 6-62-43: 2 src file/files merge ends (2022-08-03T14:40:32.324Z --> 2022-08-03T14:40:49.967Z )
[6-62-40,6-62-41,6-63-20] in progress
--- 6-63-20 (922/932) starts: 5 src file/files
--- 6-63-20: 5 src file/files merge ends (2022-08-03T14:40:50.044Z --> 2022-08-03T14:41:27.226Z )
[6-62-40,6-62-41,6-63-21] in progress
--- 6-63-21 (923/932) starts: 8 src file/files
--- 6-63-21: 8 src file/files merge ends (2022-08-03T14:41:27.227Z --> 2022-08-03T14:42:46.266Z )
[6-62-40,6-62-41,6-63-32] in progress
--- 6-63-32 (924/932) starts: 7 src file/files
--- 6-62-41: 8 src file/files merge ends (2022-08-03T14:39:51.163Z --> 2022-08-03T14:43:05.234Z )
[6-62-40,6-63-32,6-63-33] in progress
--- 6-63-33 (925/932) starts: 9 src file/files
--- 6-63-32: 7 src file/files merge ends (2022-08-03T14:42:46.270Z --> 2022-08-03T14:49:31.054Z )
[6-62-40,6-63-33,6-63-34] in progress
--- 6-63-34 (926/932) starts: 5 src file/files
--- 6-63-33: 9 src file/files merge ends (2022-08-03T14:43:05.234Z --> 2022-08-03T14:49:59.254Z )
[6-62-40,6-63-34,6-63-35] in progress
--- 6-63-35 (927/932) starts: 13 src file/files
--- 6-63-34: 5 src file/files merge ends (2022-08-03T14:49:31.056Z --> 2022-08-03T14:50:35.986Z )
[6-62-40,6-63-35,6-63-38] in progress
--- 6-63-38 (928/932) starts: 5 src file/files
--- 6-63-38: 5 src file/files merge ends (2022-08-03T14:50:35.990Z --> 2022-08-03T14:52:31.686Z )
[6-62-40,6-63-35,6-63-39] in progress
--- 6-63-39 (929/932) starts: 21 src file/files
--- 6-63-35: 13 src file/files merge ends (2022-08-03T14:49:59.255Z --> 2022-08-03T14:53:20.943Z )
[6-62-40,6-63-39,6-63-40] in progress
--- 6-63-40 (930/932) starts: 7 src file/files
--- 6-62-40: 26 src file/files merge ends (2022-08-03T14:38:50.224Z --> 2022-08-03T14:54:48.579Z )
[6-63-39,6-63-40,6-63-41] in progress
--- 6-63-41 (931/932) starts: 1 src file/files
--- 6-63-41: 1 src file/files merge ends (2022-08-03T14:54:48.584Z --> 2022-08-03T14:54:56.985Z )
[6-63-39,6-63-40,6-63-42] in progress
--- 6-63-42 (932/932) starts: 1 src file/files
--- 6-63-42: 1 src file/files merge ends (2022-08-03T14:54:56.986Z --> 2022-08-03T14:55:05.152Z )
--- 6-63-40: 7 src file/files merge ends (2022-08-03T14:53:20.944Z --> 2022-08-03T14:55:32.198Z )
--- 6-63-39: 21 src file/files merge ends (2022-08-03T14:52:31.686Z --> 2022-08-03T15:02:54.228Z )
Production ends: 2022-07-30T11:19:12.241Z --> 2022-08-03T15:02:54.312Z
Step 2-2: RGBify（RGB標高タイル化）
932の各範囲をカバーするマージされたTIFFファイルが出来ているので、これをrasterio rgbifyコマンドで変換していきます。また nodejs のspawnとbetter-queue を使うことで並列で処理をすることができるようにしています。pngで出力するとファイル数が多すぎて扱いにくい mbtiles で出力します。（サーバーに置くときにmapbox/mbtilesモジュールを使ってpngを返すことにしたいと思います。）
https://github.com/ubukawa/srtm-rgbify/blob/main/index4rgb.js
この中で、結局spawnのところがポイントなのですが、rasterio rgbifyをつかうときにbounding tileも指定して、指定範囲でクリップされるようにしています。（小ネタ：x,y,zの間にスペースがあるとコマンドが動かないので注意。）
const rgbify = spawn(rasterioPath, [
'rgbify', '-b','-10000','-i','0.1', '--max-z', maxZ, '--min-z', minZ,
'--format', 'webp', '--bounding-tile', `[${x.toString()},${y.toString()},${z.toString()}]`,
mergedPath, tmpPath
])
これも、concurrentは3 で処理を進めました。今回は最小ズームレベルは6 ，最大ズームレベルは12にしました。
8月7日AM1時でまだ932ファイルのうち600ファイルしか終わっていませんが、全部終わればかかった時間と合計サイズを報告します。
追記： 932ファイルで、179GBありました。そして変換にはだいたい6日間 かかりました。（途中少しスタックしていた時間を含む。conccurentは3）
変換が進み、所定のフォルダに6-x-y.mbtiles という名前でファイルが出てきます（x, yはそれぞれ）。
Step 2-3: mbtilesからpngタイルの配信（必要に応じ）
mbtilesフォーマットのRGB標高タイルなので、インターネット上でpngファイルとしてアクセスできるようにnodejs/expressを使って簡単なサーバーを作ります。
（範囲が狭ければ、mbutilを使ってmbtilesからzxyフォルダ構造のpngファイルに展開してもいいと思います。）
coesite3サーバーを改造するようにしました。ベクトルタイルで同じことをやっていたので、pbfをpngに変えるようなイメージでよいと思います。
...
var rgbElevRouter = require('./routes/rgbElev')
...
app.use('/unvt/rgb-elev', rgbElevRouter)
...
var express = require('express')
var router = express.Router()
const config = require('config')
const fs = require('fs')
const cors = require('cors')
const MBTiles = require('@mapbox/mbtiles')
// config constants
const defaultZ = 6
const rgbElevDir = 'rgbElev'
// global variables
let mbtilesPool = {}
let busy = false
var app = express()
app.use(cors())
//Get tile functions
const getMBTiles = async (z, x, y) => {
let mbtilesPath = ''
if (z < 6) {
mbtilesPath = `${rgbElevDir}/0-0-0.mbtiles`
} else {
mbtilesPath =
`${rgbElevDir}/6-${x >> (z - 6)}-${y >> (z - 6)}.mbtiles`
}
return new Promise((resolve, reject) => {
if (mbtilesPool[mbtilesPath]) {
resolve(mbtilesPool[mbtilesPath].mbtiles)
} else {
if (fs.existsSync(mbtilesPath)) {
new MBTiles(`${mbtilesPath}?mode=ro`, (err, mbtiles) => {
if (err) {
reject(new Error(`${mbtilesPath} could not open.`))
} else {
mbtilesPool[mbtilesPath] = {
mbtiles: mbtiles, openTime: new Date()
}
resolve(mbtilesPool[mbtilesPath].mbtiles)
}
})
} else {
reject(new Error(`${mbtilesPath} was not found.`))
}
}
})
}
const getTile = async (mbtiles, z, x, y) => {
return new Promise((resolve, reject) => {
mbtiles.getTile(z, x, y, (err, tile, headers) => {
if (err) {
reject()
} else {
resolve({tile: tile, headers: headers})
}
})
})
}
/* GET Tile. */
router.get(`/zxy/:z/:x/:y.png`,
async function(req, res) {
busy = true
const z = parseInt(req.params.z)
const x = parseInt(req.params.x)
const y = parseInt(req.params.y)
getMBTiles(z, x, y).then(mbtiles => {
getTile(mbtiles, z, x, y).then(r => {
if (r.tile) {
res.set('content-type', 'image/png')
//res.set('content-encoding', 'gzip')
res.set('last-modified', r.headers['Last-Modified'])
res.set('etag', r.headers['ETag'])
res.send(r.tile)
busy = false
} else {
res.status(404).send(`tile not found: /zxy/${z}/${x}/${y}.png`)
busy = false
}
}).catch(e => {
res.status(404).send(`tile not found (getTile error): /zxy/${z}/${x}/${y}.png`)
busy = false
})
}).catch(e => {
res.status(404).send(`mbtiles not found for /zxy/${z}/${x}/${y}.png`)
})
}
);
module.exports = router;
とすることでRGB標高タイルにアクセス出来るようになったようです。タイルサイズが512×512ピクセルになっているのはちょっと調査が必要です。rasterio rgbifyはベースが256だったのだと思うのですが、mapbox/mbtilesの働きでそうなったのか、この先ちょっと調べてみようと思います。
まとめ
作業日数と、データの保存容量を必要としますが、SRTMデータから自分でRGB標高タイル（ZL12まで）を作ることができました。3次元表現のTerrain model や Hillshadeのために使うことが出来ました。
課題
ただし、もともとのSRTMデータに穴があることがわかりました（例えばイタリアのアルプス（アオスタ地方など））。どうしたものか考えてみます。SRTM V3の記載に欠損領域は埋めてあると書いてあったんですが・・・。とりあえずUSGSのカスタマーセンターに聞いてみています。
謝辞
SRTMデータを公開しているUSGSに感謝します。