スペースシャトルの標高データSRTMで富士山近辺の立体地形図を描いてみた
毛利さんも参加した、スペースシャトルのCバンド・Xバンド干渉合成開口レーダーで測量された標高データSRTMで富士山近辺の立体地形図を描いてみました。
SRTMの説明はWikipediaさんから引用するとして
Shuttle Radar Topography Mission - Wikipedia
Shuttle Radar Topography Mission (シャトル レーダー トポグラフィー ミッション、SRTM)はスペースシャトルに搭載したレーダーで、地球の詳細な数値標高モデルを作製することを目的としたミッションである。
2000年、毛利衛も参加したエンデバーのSTS-99ミッションで行なわれ、ペイロードベイのメインアンテナと長さ60mのマストの先の外部アンテナで構成されるCバンド及びXバンド干渉合成開口レーダーで、緯度60度内を測量し陸地の80%の標高データを得た。Cバンドの標高データは無償でダウンロードでき、3秒角(約90m)メッシュのSRTM-3、30秒角(約900m)メッシュのSRTM-30、及びアメリカ国内の1秒角(約30m)メッシュのSRTM-1が公開されている。
高解像度の地形標高データは、これまで米国内でしか公開されていなかったが、2015年から世界中で公開するとホワイトハウスが2014年9月に発表した。SRTMで得られた解像度90mの低解像度データは2003年に公開されていたが、オリジナルデータの解像度30mのものが公開される。
観測期間中、各地でコーナーキューブで電波を反射させ地形図に文字を書き込む実験が行われた。
前回まで、国土地理院の標高データを使ってましたが、国土地理院の標高データは著作権上、色々制約が多く、ネット上で利用するとなるとなかなかに使いづらいので、パブリック・ドメインで誰でも利用可能なSRTMを使ってみることに。
最新データは以下で取得できるのですが、利用登録が必要なので
SRTM1 Version 3
https://e4ftl01.cr.usgs.gov/MEASURES/SRTMGL1.003/2000.02.11/
SRTM3 Version 3
https://e4ftl01.cr.usgs.gov/MEASURES/SRTMGL3.003/2000.02.11/
地形データがところどころ欠落しているようなのですが、今回は、古いVer2.1のデータを使ってみることに。
今回使用したデータは以下
https://dds.cr.usgs.gov/srtm/version2_1/SRTM3/Eurasia/N35E138.hgt.zip N35E138.hgt Driver: SRTMHGT/SRTMHGT File Format Size is 1201 x 1201 x 1 Projection is GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AUTHORITY["EPSG","4326"]] Origin = (137.999583333, 36.0004166667) Pixel Size = (0.000833333333333, -0.000833333333333) Band Type=Int16 Min=-12.000, Max=3130.000
結果はこんな感じ
Ver2.1の古いデータなので、データがところどころ欠落捨てるところがあるみたいですね。
ソースは以下
import requests import zipfile from StringIO import StringIO import struct from osgeo import gdal import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D import matplotlib.colors as colors from matplotlib import cm # SRTM3データをネットからダウンロード url = 'https://dds.cr.usgs.gov/srtm/version2_1/SRTM3/Eurasia/N35E138.hgt.zip' print url response = requests.get(url) zipfile = zipfile.ZipFile(StringIO(response.content), 'r') # ZIPから必要ファイルを一時ファイルに一旦書き出し fname = zipfile.infolist()[0].filename zipdata = zipfile.read(fname) print fname tmpf = 'tmp/'+fname f = open(tmpf, mode='w') f.write(zipdata) f.close() # SRTMファイルを読み込む dataset = gdal.Open(tmpf , gdal.GA_ReadOnly) # データセット情報の取得 print("Driver: {}/{}".format(dataset.GetDriver().ShortName, dataset.GetDriver().LongName)) print("Size is {} x {} x {}".format(dataset.RasterXSize, dataset.RasterYSize, dataset.RasterCount)) print("Projection is {}".format(dataset.GetProjection())) geotransform = dataset.GetGeoTransform() if geotransform: print("Origin = ({}, {})".format(geotransform[0], geotransform[3])) print("Pixel Size = ({}, {})".format(geotransform[1], geotransform[5])) # ラスタバンドをフェッチする band = dataset.GetRasterBand(1) print("Band Type={}".format(gdal.GetDataTypeName(band.DataType))) min = band.GetMinimum() max = band.GetMaximum() if not min or not max: (min,max) = band.ComputeRasterMinMax(True) print("Min={:.3f}, Max={:.3f}".format(min,max)) if band.GetOverviewCount() > 0: print("Band has {} overviews".format(band.GetOverviewCount())) if band.GetRasterColorTable(): print("Band has a color table with {} entries".format(band.GetRasterColorTable().GetCount())) # 標高データの読み込み Raster = np.array([]) for Y in range(dataset.RasterYSize): scanline = band.ReadRaster(xoff=0, yoff=Y, xsize=band.XSize, ysize=1, buf_xsize=band.XSize, buf_ysize=1, buf_type=gdal.GDT_Float32) #tuple_of_floats = struct.unpack('f' * band.XSize, scanline) Raster = np.append(Raster, struct.unpack('f' * band.XSize, scanline)) Raster = np.reshape(Raster, (band.XSize, band.YSize)) # 標高0メートル以下はデータ欠損によるものかもしれないので、取り敢えず0メートルに置換して標高データとする Z = np.where(Raster < 0, 0, Raster) #X,Y軸のグリッドを生成 X, Y = np.meshgrid(np.linspace(0, band.XSize-1, band.XSize), np.linspace(band.YSize-1, 0, band.YSize)) # メッシュを間引く X=X[::3,::3] Y=Y[::3,::3] Z=Z[::3,::3] # 立体地図を作成する fig = plt.figure(figsize=(16, 12), dpi=72, facecolor='w', edgecolor='k') plt.subplots_adjust(left=0, right=1, bottom=0, top=1) ax = fig.add_subplot(111, projection='3d') ax.view_init(elev=65, azim=-70) # 横を1/0.5=2倍長く設定 ax.set_aspect(0.5, adjustable='box') # 水深・標高用カラーマップ作成 cmap = colors.ListedColormap(['blue', 'royalblue', 'cornflowerblue', 'dodgerblue', 'limegreen', 'g', 'darkkhaki', 'papayawhip', 'moccasin', 'wheat', 'lightsalmon', 'tomato', 'orangered', 'red', 'brown', 'snow']) cmap.set_over('white') cmap.set_under('darkblue') bounds = [-3000, -2000, -1000, -500, 0, 10, 20, 70, 100, 200, 400, 800, 1200, 1600, 2000, 3000, 4000] norm = colors.BoundaryNorm(bounds, cmap.N) # 標高に応じて塗りつぶす surf = ax.plot_surface(X, Y, Z, norm=norm, rstride=1, cstride=1, linewidth=0, alpha=.7, cmap=cmap) # 標高100間隔で等高線を描く elevation = range(-10000,10000,100) cont = ax.contour(X, Y, Z, levels=elevation, edgecolor='none', cmap='binary', alpha=1.0) # 水深・標高バーをつける cb = fig.colorbar(surf, shrink=0.4, aspect=20, extend='both', ticks=bounds, spacing='proportional') cb.set_label('elevation') plt.savefig('1.jpg', dpi=72) plt.show()