*[[SRTMデータの可視化]] [#a59dbce3]

SRTM(Shuttle Radar Topography Mission)は、シャトルに搭載される合成開口レーダーを用いたリモートセンシング技術により取得された地表のレーダー画像です。

可視化ソフトウェア[[ParaView:http://www.paraview.org]]を使用して、このSRTMデータを可視化してみましょう。

(注:簡単な変換で可視化するため、変換後の座標などは不正確です)

(注:データとプログラムを最下段のhgtToCSV.zipに置きました。修正してご利用ください)
(注:データとプログラム(Windows用)を最下段の[[hgtToCSV.zip:https://www.rccm.co.jp/icem/pukiwiki/index.php?SRTM%E3%83%87%E3%83%BC%E3%82%BF%E3%81%AE%E5%8F%AF%E8%A6%96%E5%8C%96#k551ac14]]に置きました。run.batを修正して、実行すると動作します)

**SRTMデータの取得 [#pbcc3c07]

SRTMデータは、[[Shuttle Radar Toplography Mission:http://www2.jpl.nasa.gov/srtm/]]からダウンロードできます。

今回は、この中の
https://dds.cr.usgs.gov/srtm/version2_1/SRTM3/Africa/
から、適当に選んだN32E020.hgt.zipを使用します。

このファイルをダウロード、解凍すると
このファイルをダウンロード、解凍すると

 N32E020.hgt

という、HGT形式のファイルが格納されています。

**SRTM(HGT形式)からCSVファイルへの変換 [#f9d4e9c7]

SRTMのHGT形式は、ビッグエンディアン(Big endian)で記入された1201x1201個の座標点からなるバイナリファイルです。

以下のPythonの簡単なプログラムでCSVファイルに変換します。
以下の簡単なPythonプログラムでCSVファイルに変換します。

***hgtToCSV.py [#t46193d5]

#code(python){{

import sys

nrows = ncols = 1201
xdim = ydim = 0.000833

args = sys.argv
filename = args[1]
name, ext = os.path.splitext(filename)
outname = name + '.csv'

longitude = float(name[1:3])
latitude = float(name[4:7])

print('longitude: ' + name[0] + name[1:3])
print('latitude: ' + name[3] + name[4:7])

ulxmap = longitude + (xdim/2.0)
ulymap = (latitude+1.0) - (ydim/2.0)

a = array.array('h')
f = open(filename, 'rb')

a.fromfile(f, nrows*ncols)
f.close()

a.byteswap()
alist = a.tolist()

out = open(outname, 'w')

count = 0
for row in range(nrows):
  for col in range(ncols):
    x = ulxmap + xdim * col
    y = ulymap - ydim * row
    z = max(0.0, alist[count])
    out.write(str(x) + ',' + str(y) + ',' + str(z/10000.0) + '\n')
    count += 1

    if count%100000 == 0:
        print(str(int(100.0 * count/(nrows*ncols))) + '%')

print('100%')
print('node number: ' + str(count))

out.close()

}}

以下のコマンド
  > python hgtToCSV.py N32E020.hgt
を実行すると、N32E020.csvに変換されます。

注意:座標は不正確です。実際に使用するときは、正しい座標に変換してください。


**CSVファイルからVTSファイルへの変換 [#f399ec2c]

次に、CSVファイルをParaViewで読み込めるVTS(VTKの構造格子用データ形式)に変換します。
次に、CSVファイルをParaViewで読み込めるVTS形式(VTKの構造格子用データ形式)に変換します。

こちらも簡単なPythonスクリプトを使用します。なお、標高データは変数zに格納しています。

***csvToVts.py [#d82f8293]

#code(python){{

import os.path
import sys
#### import the simple module from the paraview
from paraview.simple import *
#### disable automatic camera reset on 'Show'
paraview.simple._DisableFirstRenderCameraReset()

args = sys.argv
filename = args[1]

workingDir = './'
fileName = args[1]

# create a new 'CSV Reader'
csvReader = CSVReader(FileName=[workingDir + fileName + '.csv'])
csvReader.DetectNumericColumns = 1
csvReader.UseStringDelimiter = 1
csvReader.HaveHeaders = 0
csvReader.FieldDelimiterCharacters = ','
csvReader.MergeConsecutiveDelimiters = 0

# create a new 'Table To Structured Grid'
tableToStructuredGrid = TableToStructuredGrid(Input=csvReader)
tableToStructuredGrid.WholeExtent = [0, 1200, 0, 1200, 0, 0]
tableToStructuredGrid.XColumn = 'Field 0'
tableToStructuredGrid.YColumn = 'Field 1'
tableToStructuredGrid.ZColumn = 'Field 2'

# create a new 'Calculator'
calculator = Calculator(Input=tableToStructuredGrid)
calculator.AttributeMode = 'Point Data'
calculator.CoordinateResults = 0
calculator.ResultNormals = 0
calculator.ResultTCoords = 0
calculator.ResultArrayName = 'z'
calculator.Function = 'coordsZ'
calculator.ReplaceInvalidResults = 1
calculator.ReplacementValue = 0.0

# save data
SaveData(workingDir + fileName + '.vts', proxy=calculator, 
           Writealltimestepsasfileseries=0, 
           DataMode='Appended', 
           EncodeAppendedData=1, 
           CompressorType='ZLib')

}}

実行は、ParaViewのインストールフォルダ
-C:/Program Files/ParaView 5.2.0-Qt4-OpenGL2-Windows-64bit/bin/pvpython.exe

にあるpvpythonを使用して実行します。
  > pvpython csvToVts.py N32E020

変換が終了すると、N32E020.vtsが出力されます。

**ParaViewによる可視化 [#l13c3534]

上記で作成されたN32E020.vtsファイルをParaViewに読み込むと、以下のように可視化できます。

#ref(N32E020.png,center,nowrap,60%,N32E020);

**ダウンロード [#k551ac14]
hgtToCSVプログラム一式:
#ref(hgtToCSV.zip)


トップ   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS