プログラム日記φ(..)

おもにPython関係のプログラムメモ

X3DOMでグリグリ動く3D世界地図を作ってみる。

こういうのが欲しかった!

 

X3DOMのBVHRefinerタグに標高データと地図画像テクスチャと法線マップをタイル形式で与えてあげると、簡単にグリグリ動く3Dマップを作れることを発見!

https://doc.x3dom.org/tutorials/largeModels/BVHRefiner/index.html

 

標高データと法線マップは只今作成中で、以下、地形の凹凸はあってませんが、先日、X3Dで地球儀を作った際のNASAの地球画像から、地図画像タイルを生成して

 

memomemokun.hateblo.jp

 

3D世界地図を表示してみるとこんな感じ。

 

http://chamu.org/map/

ご覧の環境では、object要素がサポートされていないようです。外部文書を別ウィンドウで開いてください

 

 

自由自在に地形の拡大や

f:id:memomemokun:20190118232024j:plain

 

ワイヤーフレームの表示もできちゃいます。

f:id:memomemokun:20190118232144j:plain

 

 

作り方は、少々お待ちください。

 

地球観測衛星「だいち」の標高モデルやETOPO1グローバル地形データセットの水深データを使い標高データと法線マップを作成し、

 

 

memomemokun.hateblo.jp

 

memomemokun.hateblo.jp

 

 

OpenStreetMapの地図画像とかランドサットの衛星画像などを組み合わせると、完璧、ネット上で自由に公開できる自分の目的に応じた3D世界地図を作れそうですね〜〜!

 

memomemokun.hateblo.jp

 

 

なを、上の3D世界地図の地図画像タイルは、NASAのBlue Marbleの画像を使い、ズームレベル5までのタイル画像を生成し表示しています。

NASA Visible Earth: December, Blue Marble Next Generation w/ Topography and Bathymetry

 

#!/usr/bin/env python
# -*- coding: utf-8 -*- 

import os

import numpy as np
from scipy import misc
from PIL import Image

from math import pi
from math import tanh
from math import sin
from math import asin
from numpy import arctanh


# ピクセル座標を緯度経度に変換する
def fromPointToLatLng(pixelLat,pixelLon, z):
    L = 85.05112878
    lon = 180 * ((pixelLon / 2.0**(z+7) ) - 1)
    lat = 180/pi * (asin(tanh(-pi/2**(z+7)*pixelLat + arctanh(sin(pi/180*L)))))
    return lat, lon
    

# イメージ処理可能サイズを21600*10800に指定
Image.MAX_IMAGE_PIXELS=21600*10800

texturefile = 'world.topo.bathy.200412.3x21600x10800.jpg'
im = Image.open(texturefile)


for z in range(6):

    # ズームレベル用ヂレクトリーを作成する
    zoom_dir = 'map/texture/%d' % z
    if not os.path.isdir(zoom_dir): os.mkdir(zoom_dir)
    
    # ズームレベル描画緯度経度範囲にblue marble画像を切り出し
    lat, lon = fromPointToLatLng(0, 0, z)
    upper = im.size[1] * (90 - lat) / 180
    lower = im.size[1] * (90 + lat) / 180
    im_crop = im.crop((0, upper, im.size[0], lower))
    
    # 縦横のタイル分割数
    partition = 2**z

    # タイル画像の分割サイズ
    xsize, ysize = im_crop.size
    xsize = float(xsize / partition)
    ysize = float(ysize / partition)
    print im_crop.size, xsize, ysize

    
    # 地図画像を分割する
    for column in range(partition):
        column_dir = 'map/texture/%d/%d' % (z, column)
        if not os.path.isdir(column_dir): os.mkdir(column_dir)

        for row in range(partition):
            
            left  = xsize * column
            upper = ysize * row
            right = xsize * (column + 1)
            lower = ysize * (row + 1)
            
            print left, upper, right, lower
            
            im_cropRow = im_crop.crop((left, upper, right, lower))
            img_resize = im_cropRow.resize((256,256), Image.LANCZOS)

            # 出力タイル名
            tile = 'map/texture/%d/%d/%d.jpg' % (z, column, row)
            img_resize.save(tile, quality=100)