気になるメモメモφ(..)

日常気になる事を取り留めもなくメモメモしてます。

地球観測衛星「だいち」の標高モデルで立体富士山を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の3Dプロットも高速描画のPython Mayaviを使用。

 

memomemokun.hateblo.jp

 

 

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

ALOS全球数値地表モデルで富士山を描いてみる(1)

ALOS全球数値地表モデルで富士山を描いてみる(1)

 

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

ALOS全球数値地表モデルで富士山を描いてみる(2)

ALOS全球数値地表モデルで富士山を描いてみる(2)

 

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

ALOS全球数値地表モデルで富士山を描いてみる(3)

ALOS全球数値地表モデルで富士山を描いてみる(3)

 

ソースは以下

 

#!/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()