バーチャル長崎構築RTAをする
長崎の3D点群データの一部が3/14から公開されたようなので、これを使ってバーチャル長崎を作ってみます
Termsを見る限りCC-BYのようなので表記を行なって利用していきます
ただ、
①データセットの作成者欄に記載する組織名（長崎県等）
②リソースの名称
③リソースのURL
果たしてこの３つに何が入るのか...その謎を探るべくアマゾンの奥地へと...とは行かないので一旦以下の表記を利用します。
開始時刻 3/14 18:10
終了時刻 3/16 20:00
点群を一旦全部ダウンロードする
全データで260GBぐらいあるため、自宅PCにダウンロードするとやばいのでさくらのクラウド(東京リージョン)を利用します
最初DigitalOceanを使ってましたが、帯域が細くて辛いのでやめました。(12:30から始めて17時時点で半分も行ってなかったため、再走しました。)
スクリプトを書く
詳しくは書きませんが、DevToolsを開いたままMapに行くといいかんじにデータ一覧を取得してるリクエストがあるのでこのjsonを保存しておきます。
CDNとかが前段にないApacheサーバのようなので、流石に一気にダウンロードのリクエストを送るのはやばいと言うことでPythonでスクリプトを書いてある程度自重した感じでダウンロードします。
urlretrieve使おうかなと思ったらこんな記事を見つけたので、urlopenを使って書いていきます
import json
import urllib.request
import os
import time
temp_folder = "./tmp"
las_folder = "./data"
os.makedirs(temp_folder, exist_ok=True)
os.makedirs(las_folder, exist_ok=True)
i = 0
with open("./response.json", "r") as f:
json_data = json.load(f)
value_length = len(json_data.get("result").values())
for v in json_data.get("result").values():
file_name = v.get("path").split("/")[-1]
remote_path = v.get("path").replace("data/", "https://opennagasaki.nerc.or.jp/storage/") + ".zip"
temp_file_name = temp_folder + "/" + v.get("path").split("/")[-1] + ".zip"
data = urllib.request.urlopen(remote_path).read()
with open(temp_file_name, mode="wb") as fa:
fa.write(data)
i += 1
print("{}/{} : {}".format(i, value_length, file_name))
if i % 10 == 0:
time.sleep(10)
気休め程度かもしれませんが、10回に1回、10秒sleepするようにしておきます。
流石に並列でリクエストしなければ大丈夫だとは思いますが、リリース初日でアクセスが集中している可能性も考えてこうしていました。
開始時刻 3/14 18:20
終了時刻 3/14 18:30
ダウンロードする
上で作ったスクリプトを回すだけですね。
めっちゃ時間かかるのでscreenとかで別画面にしておく方が良いと思います。
python downloader.py
開始時刻 3/14 18:30
終了時刻 3/15 01:15
解凍する
tmpというフォルダにzipがダウンロードされているので解凍していきましょう
6コアなので単純に6並列で処理します。(diskIO的に多分あんまり宜しくはない気がする)
find ./tmp -name "*.zip" | xargs -P 6 -IX -n 1 unzip -d ./data X
開始時刻 3/15 01:15
終了時刻 3/15 02:03
点群を3DTilesに変換する
点群を一括でrepareする
何故か毎回opendataのlasデータはぶっ壊れてるので簡易的にlasinfoのrepairにだけかけます
git clone https://github.com/LAStools/LAStools.git
cd LAStools
make
cd ./bin
sudo cp lasinfo /usr/local/bin
find ./data -name '*.las' | xargs -P 6 -IX -n 1 lasinfo -i X -repair
find ./data -name '*.las' | xargs -IX -n 1 pdal info --summary X >> test.log
一応エラーが出ないことだけはチェック
開始時刻 3/15 02:03
終了時刻 3/15 03:48
pntsファイルを作成する
py3dtilesをつかってlasからpntsに変換します
今後差分が出る可能性が大いにあるため、それぞれのpntsを作ってからmergeしたかったんですが、mergeコマンドの使い方がわからなすぎたので諦めて、一旦全部一気に作ろうかと思います。
睡眠を取るために、2core 4gbのインスタンスで以下のコマンドを実行してから寝た(3/15 03:48)んですが、
py3dtiles convert --srs_in 6669 --srs_out 4978 --out ./3dtiles ./data/*.las
起きたらメモリ不足でKilledになっていました。(3/15 06:30)
そこでさくらのクラウドの20Core メモリ32GBのインスタンスに変更し、それらをフル活用することにします。
一部Memory不足に陥る可能性があるので念の為cache_sizeを30720MB (= 30GB)に規制します
py3dtiles convert --srs_in 6669 --srs_out 4978 --cache_size 30720 --out ./3dtiles ./data/*.las
使おうとしていたコマンド
mkdir 3dtiles
find ./data -name '*.las' | awk '{print gensub(/\.\/data\/(.+)\.las/, "\\1", "g")}' | xargs -P 6 -IX -n 1 py3dtiles convert --srs_in 6669 --srs_out 4978 --cache_size 30720 --out 3dtiles/X --jobs 6 data/X.las
開始時刻 3/15 10:15
終了時刻 3/15 16:36
表示してみる
出来上がったので表示してみましょう
海の上やんけ！！！！！
X,Y変換がうまくいってないですね..........
(上に書いてないですがX,Y変換してたんですがうまくいっていない...)
再走しよう
ダウンロードとか諸々済んでしまっているのでレギュレーション違反ですが、ダウンロード〜repairの時間を先に足すことで時間を出しましょう
ubuntu@ubuntu:~$ df
Filesystem 1K-blocks Used Available Use% Mounted on
tmpfs 3287232 968 3286264 1% /run
/dev/vda2 1056688032 851911240 151073544 85% /
tmpfs 16436144 0 16436144 0% /dev/shm
tmpfs 5120 0 5120 0% /run/lock
容量逼迫しちゃってますね...
物理で殴るをコンセプトにやっているので、ディスクサイズをあげようとしたんですが時間がかかりそうなので、諦めてミスったやつを削除しちゃいましょう
ubuntu@ubuntu:~$ df
Filesystem 1K-blocks Used Available Use% Mounted on
tmpfs 3287232 968 3286264 1% /run
/dev/vda2 1056688032 645687984 357296800 65% /
tmpfs 16436144 0 16436144 0% /dev/shm
tmpfs 5120 0 5120 0% /run/lock
んーいけるやろ！
開始時刻 3/15 18:15
終了時刻 3/15 18:30
諸々検証する
- XY変換のみ ng
- SRS指定 & XY変換 & py3dtilesのSRS設定 ng
- SRS指定 & XY変換 ng
- SRS指定 ok
ということでpdalでSRSだけ指定します。
開始時刻 3/15 18:30
終了時刻 3/15 19:15
pdalでSRSを指定する
超ありがたいことにわかりやすい記事があるのでこれを参考に、SRSだけ修正します
世界測地系(JGD2011)／平面直角座標系第1系 なので EPSG:6669 です。
find ./data -name '*.las' | awk '{print gensub(/\.\/data\/(.+)\.las/, "\\1", "g")}' | xargs -P 20 -IX -n 1 pdal translate -i ./data/X.las -o /nagasaki-data/translate/X.las --writers.las.a_srs="EPSG:6669"
:::失敗例
以下のjsonを使用してXY変換して試したんですが、SRSの変換のみでないと何故かさらにずれてしまうようなので一旦やらないように
[
{
"type": "readers.las",
"filename": "rescale.las",
"spatialreference": "EPSG:6669"
},
{
"type": "filters.reprojection",
"in_srs": "EPSG:6669",
"out_srs": "EPSG:6669",
"in_axis_ordering": "2, 1"
},
{
"type": "writers.las",
"filename": "translated.las",
"forward": "header,scale,vlr",
"offset_x": "auto",
"offset_y": "auto",
"offset_z": "auto"
}
]
:::
何度か試してdisk不足で落ちたのでディスクを2TB増設しました
開始時刻 3/15 19:15
終了時刻 3/16 02:18
py3dtilesで再生成する
cd /nagasaki-data
py3dtiles convert --srs_out 4978 --cache_size 30720 --out ./3dtiles ./translate/*.las
100.0 % in 13678 sec [est. time left: 0 sec]
開始時刻 3/16 02:18
終了時刻 3/16 09:22
生成完了時の最終的なサイズは
205G ./3dtiles
205GBでした。
生成を確認する
cp ./3dtiles /var/www/html/3dtiles
一旦とりあえず表示してみるようのコード
<!DOCTYPE html>
<html>
<head>
<meta charset="UTF-8">
<title>Cesium</title>
<script src="https://cesium.com/downloads/cesiumjs/releases/1.72/Build/Cesium/Cesium.js"></script>
<link href="https://cesium.com/downloads/cesiumjs/releases/1.72/Build/Cesium/Widgets/widgets.css" rel="stylesheet">
<style>
#cesiumContainer {
position: absolute;
top: 0;
left: 0;
height: 100%;
width: 100%;
margin: 0;
overflow: hidden;
padding: 0;
font-family: sans-serif;
}
html {
height: 100%;
}
body {
padding: 0;
margin: 0;
overflow: hidden;
heght: 100%;
}
</style>
</head>
<body>
<div id="cesiumContainer"></div>
<script>
const raster_tile = new Cesium.OpenStreetMapImageryProvider({
url: 'https://cyberjapandata.gsi.go.jp/xyz/std/',
credit: new Cesium.Credit('地理院タイル', '', 'https://maps.gsi.go.jp/development/ichiran.html'),
maximumLevel: 18,
})
const viewer = new Cesium.Viewer('cesiumContainer', {
imageryProvider: raster_tile,
baseLayerPicker: false,
geocoder: false,
homeButton: false
})
viewer.camera.setView({
destination : Cesium.Cartesian3.fromDegrees(140.00, 36.14, 20000000.0)
})
const tileset = viewer.scene.primitives.add(
new Cesium.Cesium3DTileset({
url: './3dtiles/tileset.json'
})
)
tileset.style = new Cesium.Cesium3DTileStyle({
pointSize: 5
})
viewer.zoomTo(tileset)
</script>
</body>
</html>
いい感じに表示されてますね！！！
Zoominした感じも良さそう
ということでこれでやっと生成が完了です。長かった....
開始時刻 3/16 09:22
終了時刻 3/16 12:16
R2にアップロードする
R2にsyncします
cd ./3dtiles
aws s3 sync . s3://bucketname/3dtiles/ --profile r2 --endpoint-url https://userid.r2.cloudflarestorage.com
基本はこれでいいと思いますが、色々とファイル数が多い(2444495ファイル)などあるので、include/excludeでよしなに分散してもらうといいと思います。
aws s3 sync . s3://bucketname/3dtiles/ --profile r2 --endpoint-url https://userid.r2.cloudflarestorage.com --exclude "*" --include "04*"
こんな感じですかね
開始時刻 3/16 12:16
終了時刻 3/16 18:30
提供用のCloudflare Workerを作成する
WranglerでInitします
> wrangler init nagasaki-3dtiles Thu Mar 16 17:54:00 2023
⛅️ wrangler 2.12.3
--------------------
Using npm as package manager.
✨ Created nagasaki-3dtiles/wrangler.toml
✔ Would you like to use git to manage this Worker? … yes
✨ Initialized git repository at nagasaki-3dtiles
✔ No package.json found. Would you like to create one? … yes
✨ Created nagasaki-3dtiles/package.json
✔ Would you like to use TypeScript? … yes
✨ Created nagasaki-3dtiles/tsconfig.json
✔ Would you like to create a Worker at nagasaki-3dtiles/src/index.ts? › Fetch handler
✨ Created nagasaki-3dtiles/src/index.ts
✔ Would you like us to write your first test with Vitest? … no
npm WARN deprecated rollup-plugin-inject@3.0.2: This package has been deprecated and is no longer maintained. Please use @rollup/plugin-inject.
npm WARN deprecated sourcemap-codec@1.4.8: Please use @jridgewell/sourcemap-codec instead
added 104 packages, and audited 105 packages in 14s
11 packages are looking for funding
run `npm fund` for details
found 0 vulnerabilities
✨ Installed @cloudflare/workers-types and typescript into devDependencies
To start developing your Worker, run `cd nagasaki-3dtiles && npm start`
To publish your Worker to the Internet, run `npm run deploy`
wrangler.toml を編集してR2をbindingします。
name = "nagasaki-3dtiles"
main = "src/index.ts"
compatibility_date = "2023-03-16"
account_id = "0977733090346c25bbe7b1d1c522e83f"
workers_dev = true
[[r2_buckets]]
binding = "BUCKET"
bucket_name = "virtual-nagasaki-data"
index.tsもCacheAPIとR2を使うように変更します
export interface Env {
BUCKET: R2Bucket;
}
export default {
async fetch(
request: Request,
env: Env,
ctx: ExecutionContext
): Promise<Response> {
switch (request.method) {
case 'GET':
const url = new URL(request.url);
const cacheKey = new Request(url.toString(), request);
const cache = caches.default;
let response = await cache.match(cacheKey);
if (response) {
return response;
}
const objectKey = url.pathname.slice(1);
const object = await env.BUCKET.get(`3dtiles/${objectKey}`);
if (object === null) {
return new Response('Object Not Found', { status: 404 });
}
const headers = new Headers();
object.writeHttpMetadata(headers);
headers.set('etag', object.httpEtag);
headers.set('Access-Control-Allow-Origin', '*')
headers.append('Cache-Control', 's-maxage=7200');
response = new Response(object.body, {
headers,
});
ctx.waitUntil(cache.put(cacheKey, response.clone()));
return response;
default:
return new Response('Method Not Allowed', {
status: 405,
headers: {
Allow: 'GET',
},
});
}
},
};
これを wrangler publish します。
合わせてコントロールパネルでCustomDomainを設定したので
でアクセスできるようになりました
開始時刻 3/16 18:30
終了時刻 3/16 18:55
実際に表示する
DeckGLでいい感じに表示します
Vite + React + DeckGLでやってみます
とりあえずプロジェクトを作ります
ここと同じことをやるので省略します
Mapライブラリを入れます
pnpm install react-map-gl maplibre-gl
pnpm install mapbox-gl@npm:empty-npm-package@1.0.0
DeckGLを入れます
pnpm install @deck.gl/react @deck.gl/core @deck.gl/geo-layers @deck.gl/extensions @loaders.gl/core @luma.gl/core @luma.gl/shadertools @deck.gl/layers @loaders.gl/gltf @luma.gl/engine @luma.gl/webgl @deck.gl/mesh-layers @loaders.gl/images @luma.gl/gltools @luma.gl/constants @math.gl/core @math.gl/web-mercator @loaders.gl/3d-tiles
App.tsxを描きます
import 'maplibre-gl/dist/maplibre-gl.css';
import { Tile3DLayer } from '@deck.gl/geo-layers/typed';
import { default as DeckGL } from '@deck.gl/react/typed';
import { Tiles3DLoader } from '@loaders.gl/3d-tiles';
import { default as maplibregl } from 'maplibre-gl';
import type { FC } from 'react';
import { default as Map } from 'react-map-gl';
const INITIAL_VIEW_STATE = {
longitude: 129.87361,
latitude: 32.74472,
zoom: 11,
pitch: 30,
bearing: 0,
};
const App: FC = () => {
const layers = [
new Tile3DLayer({
id: 'tile-3d-layer',
pointSize: 2,
data: 'https://nagasaki-3dtiles.nekoyasan.dev/tileset.json',
loaders: [Tiles3DLoader],
}),
];
return (
<DeckGL initialViewState={INITIAL_VIEW_STATE} controller={true} layers={layers}>
<Map
mapLib={maplibregl}
mapStyle={{
version: 8,
sources: {
photo: {
type: 'raster',
tiles: ['https://cyberjapandata.gsi.go.jp/xyz/seamlessphoto/{z}/{x}/{y}.jpg'],
tileSize: 256,
attribution:
'<a href="http://www.gsi.go.jp/kikakuchousei/kikakuchousei40182.html" target="_blank" rel="noopener noreferrer">地理院タイル</a>',
},
},
layers: [
{
id: 'photo',
type: 'raster',
source: 'photo',
minzoom: 2,
maxzoom: 18,
},
],
}}
/>
</DeckGL>
);
};
export default App;
これをCloudflare Pagesにデプロイしたものがこちらです
開始時刻 3/16 18:55
終了時刻 3/16 20:00
以上で完走とします！
完走した感想
２回途中再走になってしまいましたが無事完走できてよかったです。
実際にかかった時間は丸二日ほどになりましたが、実際他の走者の方は
- ダウンロード 6時間55分
- 解凍 48分
- LAS修正 by lastools 1時間45分
- LAS修正 by pdal 7時間3分
- タイル生成 7時間4分
- 生成確認 2時間54分
- タイル公開 6時間39分
- デモ作成 1時間5分
合計 34時間13分
で完走できるかと思います。
las -> 3dtilesの変換ははまりどころがかなり多く、辛い部分も多くありましたが先人の知恵のおかげでなんとか正常に表示できるところまで漕ぎ着けることができました。
長崎県全域の点群データに関してはこちらの記事をぜひご確認ください。
今回、さくらのクラウドをフル活用して行った関係もあり結果的にはかなりの額になってしまっていますが、まぁ、勉強代ということで....
地理的に近いのもありダウンロード速度はDigitalOceanと比べ物にならないなという感じです。
このようにかなり綺麗な点群データを表示することもできるのでぜひ皆さんもお試しください！
（そして、もっといいチャート・もっといい方法を見つけていただけると嬉しいです）
また、3dtilesはしばらく、クラウド破産しなさそうであればCC-BYを継承して公開させていただきますのでぜひご利用ください(長崎県の方の表記だけでも構いません)
Discussion