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

 

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

初めに

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

処理

タイルデータの取得

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

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

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

画像を付ける

手順1.画像のダウンロード

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

手順2.STLファイルを作成する

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

手順3.STL→OBJ変換する

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

手順4. 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を使いたいと思います。
/img/2020-05-14-001/2020-05-14-001-mini.jpg
ただ、これを使うには、タイルの番号ではなく、座標でしていする必要があります。
そのため、タイル番号→座標の変換をします。

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))

結果

こんな感じになりました。
標高データの問題ですが、プロメテウス火山が周りの建物も飲み込んでいますw
/img/2020-10-04-001/2020-10-04-001-009.jpg

ディズニーランドの方をアップしてみると、結構いい感じにも見えます。
意外に小高い丘みたいなのが沢山あるんですね。
/img/2020-10-04-001/2020-10-04-001-010.jpg

縮尺まで本日手が回りませんでした。
縮尺は結構簡単に計算できますが、とりあえずこのままで…。

おすすめ記事

Pandas徹底解説 - Python
Pandas徹底解説 - Python
Django エラー『TemplateDoesNotExist』対策
Django エラー『TemplateDoesNotExist』対策
音楽ファイルの時間をpythonで取得(ffmpeg使用)
音楽ファイルの時間をpythonで取得(ffmpeg使用)
Blender Bone - Python
Blender Bone - Python
文字列による条件抽出 - Python Pandas
文字列による条件抽出 - Python Pandas
SEO対策として、セキュリティー対策ソフトでチェック
SEO対策として、セキュリティー対策ソフトでチェック
Supponsered

このエントリーをはてなブックマークに追加

Comments

comments powered by Disqus