国土地理院標高タイルテキスト版を使って、ディズニーシーのプロメテウス火山を可視化

Sun, Oct 4, 2020

${{DZ_TITLE}}$
国土地理院のテキスト版標高タイルを少し見てみました。
ディズニーランドの花火が平日も復旧したので、ディズニーシーにあるプロメテウス火山を3D化しました。

初めに

以前、国土地理院の画像タイルの記事の中で、標高タイルのPNG版の仕様について言及しました。
ですが、もっと簡単な、テキストファイル版の標高タイルもあります。
なので、今回はそのテキスト版の標高タイルを取得してみたいと思います。
今回の記事ではがっつりプログラムをしないで、ラフにこんなことが出来るよ程度の紹介にしようかと思います。
後日、できたら、しっかりプログラム組んで楽ができる様にできたらと思っています。

処理

タイルデータの取得

google mapから座標、 (35.627134, 139.883577) ～ (35.625119, 139.885798)辺りにディズニーシーにあるプロメテウス火山あることを確認。
国土地理院のテキスト版標高タイルは縮尺14より細かいファイルが無いことから、一番細かい14を選択します。
前回の記事のタイルを求める計算式（下に抜粋）より、標高タイルのURLを求められます。
タイル取得部分を抜粋

def long2x_float( long , zoom ):
    """
    Logtitude(経度)より タイル画像の座標を求める。（小数点も活用の必要あり）

    Parameters
    ----------
    long float
        Logtitude(経度)
    zoom int
        拡大値
    Returns
    -------
    x
        タイル画像中の座標(1=タイル1枚分)
    """
    center_point = 2 ** ( zoom + 7 )
    total_pixels = center_point * 2
    pixels_per_long_degree = total_pixels / 360;

    return (center_point + long * pixels_per_long_degree + 0.5) / 256

def lat2y_float( lat , zoom ):
    """
    Latitude(緯度)より タイル画像の座標を求める。（小数点も活用の必要あり）

    Parameters
    ----------
    lat float
        Latitude(緯度)
    zoom int
        拡大値
    Returns
    -------
    y
        タイル画像中の座標(1=タイル1枚分)
    """
    center_point          = 2 ** ( zoom + 7 )
    total_pixels          = center_point * 2
    pixels_per_lng_radian = total_pixels / (2 * math.pi)
    siny = min(max(math.sin(lat * (math.pi / 180)), -0.9999), 0.9999);

    return       (center_point - 0.5 * math.log((1 + siny) / (1 - siny)) * pixels_per_lng_radian + 0.5) / 256

結果、下記の様なURLを求められました。
https://cyberjapandata.gsi.go.jp/xyz/dem/14/14558/6454.txt

ファイル内の確認

256x256のCSVデータで、海がeとなっているようです。
/img/2020-10-04-001/2020-10-04-001-001.jpg

pythonで可視化

海は判りやすいように0ではなく-10としておきました。単位はメートルです。

%matplotlib inline
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

df = pd.read_csv(r'14558-6454.txt',encoding='utf8',header=None)
df = df.replace('e',-10)
df = df.astype(float)
sns.heatmap(df)

/img/2020-10-04-001/2020-10-04-001-002.jpg
なんかそれっぽいのが映っていますね。

Meshlabで確認

pythonでxyzファイルを出力して、Meshlabで確認してみます。

filepath = r'out.xyz'
with open(filepath, 'w+') as fp:
    for idxy,row in df.iterrows():
        lst  = list(row)
        lst.reverse()
        for idxx,val in enumerate(lst):
            fp.write('%d %d %f\n '%(idxx,idxy,val))

正直判りづらい…。
/img/2020-10-04-001/2020-10-04-001-003.jpg
Meshlabでの面張りを調べるのもおっくうなので、別の方法で…。

OpenSCADで確認

filepath = r'out.sur'
with open(filepath, 'w+') as fp:
    for idx1,row in df.iterrows():
        lst  = list(row)
        lst.reverse()
        for val in lst:
            fp.write('%f '%(val))
        fp.write('\n')

OpenSCADでは簡単にテキストデータを読み込めるので、Pythonデータを修正します。

OpenSCADはプロラムで動かすCADなので、上記のデータを読み込むようなプログラムを組みます。
Windowsでフォルダ区切りの\ を使う場合 \\のように2回記載します。
\=￥記号です。

surface(file = "out.sur", center = true, convexity = 5);

見えました！
/img/2020-10-04-001/2020-10-04-001-004.jpg

小さな山なので画質が低いのは仕方ないですね。
ここまでしっかり表示されただけで良しとしましょう。
/img/2020-10-04-001/2020-10-04-001-005.jpg

右の山はディズニーシーにあるプロメテウス火山で右の小高いものは…？
インディージョーンズですかね？

画像を付ける

手順１．画像のダウンロード

国土地理院の画像ファイルをとりあえず貼り付けてみます。
画像のURLは置き換えるだけで単純なので、簡単ですね。
demをortに、txtをjpgにします。
http://cyberjapandata.gsi.go.jp/xyz/ort/14/14558/6454.jpg

手順２．STLファイルを作成する

OpenSCADからSTLを出力
（この程度のファイルなのに、凄く時間がかかります…）

手順３．STL→OBJ変換する

MeshLabでSTLファイルを読み込んで、OBJファイルに変換します。

手順４． OBJにテクスチャを貼り付ける

以前書いた記事よりOBJファイルにテクスチャを貼り付けます。
${{DZ_TITLE}}$

どう見ても180度ずれていますね…。
/img/2020-10-04-001/2020-10-04-001-006.jpg
OpenSCADでZ軸に180度回転させます。
※ 画像回転させる方が楽なのですが、標高を回転させる方が筋だと思うので、このようにします。
rotate(180,[0,0,1]) データ読込のsurfaceの前に定義するだけです。
めちゃくちゃ簡単ですが、またSTL出力に結構時間がかかるのが地味に嫌です。

rotate(180,[0,0,1])
surface(file = "out.sur", center = true, convexity = 5);

で、上記の手順を再度実施して、今度は綺麗にできました。
/img/2020-10-04-001/2020-10-04-001-007.jpg

アップ画像
よく見てみると、小さいですがディズニーランドのビックサンダーマウンテンと思われるものも映っています。
/img/2020-10-04-001/2020-10-04-001-008.jpg

画質の改善

標高メッシュタイルと同じ縮尺14で画像も用意しましたが、この写真のデータは縮尺18まであります。
縮尺4違うと16倍の長さになるので、画像ファイル16x16、画素4096x4096のデータになります。
大きすぎるサイズではないので、この縮尺18のデータを用意します。

イメージデータの作成方法

以前記事にした画像を作るClassを使いたいと思います。

ただ、これを使うには、タイルの番号ではなく、座標でしていする必要があります。
そのため、タイル番号→座標の変換をします。

import math
def x2long_float( x , zoom ):
    """
    タイル画像の座標より、Logtitude(経度)を求める。（小数点も活用の必要あり）

    Parameters
    ----------
    int x
        タイル画像中の座標(1=タイル1枚分)
    zoom int
        拡大値
    Returns
    -------
    long int
        Logtitude(経度)

    """
    center_point = 2 ** ( zoom + 7 )
    total_pixels = center_point * 2
    pixels_per_long_degree = total_pixels / 360;

    return (256 * x - center_point-0.5) /pixels_per_long_degree

def y2lat_float( y , zoom ):
    """
    タイル画像の座標よりLatitude(緯度)を求める。

    Parameters
    ----------
    y int
        タイル画像中の座標(1=タイル1枚分)
    zoom int
        拡大値
    Returns
    -------
    lat float
        Latitude(緯度)
    """
    y_temp = (y / 2.0**zoom) * 2 * math.pi - math.pi
    lat1 = 2 * math.atan(math.e ** (- y_temp)) * 180 / math.pi - 90

    return lat1

x2long_float(14558,14),y2lat_float(6454,14),x2long_float(14559,14),y2lat_float(6455,14)

で結果以下のようにも止まりました。
(139.87788677215576, 35.63944106897392, 139.89985942840576, 35.62158189955967)

このデータを以前の記事のClassで活用します。

pos1         = (139.87788677215576, 35.63944106897392)
pos2         = (139.89985942840576, 35.62158189955967)
zoom         = 18
out_filepath = r'D:\temp\out2.jpg'
GetImage.make_file( pos1, pos2, zoom, out_filepath )

# Jupyter Notebook環境下で画像を表示
from IPython.display import Image, display_jpeg
display_jpeg(Image(out_filepath))