地球観測衛星「だいち」の標高モデルで立体富士山をPythonで描いてみる

宇宙航空研究開発機構（JAXA）の陸域観測技術衛星「だいち」（ALOS）による、30mメッシュの無料の標高データセット(DSM) 「ALOS World 3D - 30m (AW3D30)」で立体富士山を描いてみました

ALOS全球数値地表モデル(DSM) "ALOS World 3D - 30m (AW3D30)"

プロダクト: ALOS World 3D - 30m (AW3D30)
第2.1版解像度: 緯度・経度1秒 (30m相当)
配布単位: 緯度1度×経度1度タイル
高さ精度: 目標5m (標準偏差)

今回使用したのは富士山周辺の以下データで、

ラスタサイズ　= 3600 x 3600 x 1
左上経緯度 = 138.0, 36.0
解像度 = 0.000277777777778, -0.000277777777778
回転方向 = 0.0, 0.0
データタイプ = Int16

描画には3600 x 3600の３Dプロットも高速描画のPython Mayaviを使用。

memomemokun.hateblo.jp

これだけの解像度で、無料で自由に使えるのはありがたいですねー！

Mayaviの描画画面で立体富士山をいろんな方向からぐりぐりしてみる（その１）

Mayaviの描画画面で立体富士山をいろんな方向からぐりぐりしてみる（その２）

ソースは以下

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

import numpy as np
from osgeo import gdal

raster_file = 'N035E138_AVE_DSM.tif'
dataset = gdal.Open(raster_file , gdal.GA_ReadOnly)

# データセット情報の取得
xsize = dataset.RasterXSize
ysize = dataset.RasterYSize

print("ラスタサイズ 　= {} x {} x {}".format(xsize,
                                    ysize,
                                    dataset.RasterCount))

geotransform = dataset.GetGeoTransform()
print("左上経緯度     = {}, {}".format(geotransform[0], geotransform[3]))
print("解像度         = {}, {}".format(geotransform[1], geotransform[5]))
print("回転方向       = {}, {}".format(geotransform[2], geotransform[4]))


# ラスタバンドをフェッチする
band = dataset.GetRasterBand(1)
dataType = band.DataType
dataTypeName = gdal.GetDataTypeName(dataType)
print("データタイプ = {}".format(dataTypeName))

# ラスタデータ取得
Raster = band.ReadAsArray(0, 0, xsize, ysize).astype(np.int16)

from mayavi import mlab

mlab.figure(size=(1200, 1000), bgcolor=(0.8, 0.8, 1.0))
mlab.surf(Raster, colormap='gist_earth', warp_scale=0.05,
            vmin=-100, vmax=2500)

mlab.view(azimuth=15.9, elevation=83, distance=9000)
mlab.savefig('1.jpg')
mlab.show()

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

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

地球観測衛星「だいち」の標高モデルで立体富士山をPythonで描いてみる