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

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

OpenStreetMapの地図画像で立体地図をつくってみました。

Googleマップにしろ、国土地理院の地図にしろ、その画像をネット上で自由に加工して配信したいとなると、その利用には何かと著作権が絡むので、みんなで使えるフリーマップOpenStreetMapを加工し、立体地図をPythonで作成してみました。

 

OpenStreetMap Japan | 自由な地図をみんなの手に/The Free Wiki World Map

 

 

もともとのOpenStreetMapの画像に

http://a.tile.openstreetmap.org/10/906/404.png

 

 

 

『pythonで富士山の斜面勾配図ぽいものを描いてみる』の記事で作成した方法で射影陰影図を作成して合成すると

 

 

 

こんな感じになります。

富士山近辺の立体地図

 

 

また、国土地理院の起伏を示した地図に先のOpenStreetMapの地図画像を合成してみると、こんな感じになりました。

https://cyberjapandata.gsi.go.jp/xyz/relief/10/906/404.png

 

 

富士山近辺の立体地図

 

ソースはお待ちください。

 

スペースシャトルの標高データ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のものが公開される。

 

 観測期間中、各地でコーナーキューブで電波を反射させ地形図に文字を書き込む実験が行われた。

 Wikipedia内でもスペースシャトル標高データとランドサット衛星写真を合成した画像が公開されている

 

 

前回まで、国土地理院の標高データを使ってましたが、国土地理院の標高データは著作権上、色々制約が多く、ネット上で利用するとなるとなかなかに使いづらいので、パブリック・ドメインで誰でも利用可能な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のデータを使ってみることに。

http://dds.cr.usgs.gov/srtm/

 

 

今回使用したデータは以下

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の古いデータなので、データがところどころ欠落捨てるところがあるみたいですね。 

SRTP標高データによる富士山付近の立体地形図

 

ソースは以下

 

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()