pythonで富士山の斜面勾配図ぽいものを描いてみる（numpy.gradientで遊んでみた）

先日から、国土地理院の数値標高モデルからpythonで地形図を描いてみるのにハマってますが^^::

memomemokun.hateblo.jp

地形図というのは調べてみると凄く奥行きが深いもののようで、ちょっとだけ調べてみたら、地形図には、傾斜量図、 斜面勾配図、 斜面方位図、 地上開度図、地下開度図、尾根谷度図、接峰面図、接谷面図、起伏量図、傾斜ラプラシアン図、高度分散量図、高度分散異常図etc.などなど、地形の特徴をを読み解くためのいろんな解析図があるようで、

それらは、おいおい調べていくことにして、上に出てきた傾斜量図、 斜面勾配図、 斜面方位図という単語から、

勾配と言ったら、NumPyにはスカラー場の勾配（gradient：グラディエント）を求めるためのお誂え向けの配列要素間の勾配を求める関数があるので

numpy.gradient — NumPy v1.15 Manual

標高数値メッシュの隣接する要素間の勾配をgradientで求めて、図にして遊んでみました。

今回も図にしてみるのは、富士山周辺の標高データ。（ズームレベル12、x座標3624〜3628、y座標1615〜1619の標高タイル）

gradientで求めた東西方向の勾配の値から、西側に勾配が大きく下がっている斜面ほど濃いグレーで図示して見たのが以下。なんか、それっぽい、凹凸図になってる！

同様に、gradientで求めた南北方向の勾配の値から、北側が下がっている斜面ほど濃いグレーで図示して見たのが以下。

東西方向の勾配と南北方向の勾配を足して２で割って、上の、２つをあんちょこ合成してみたのが、以下で、南東方向から光が当たって凹凸で陰影ができてる感じでしょうか。

ちょっと、地形解析についてちゃんと調べてみよーっと！

上を描いた時のソースは以下。

import requests
from io import StringIO
import string
import pandas as pd
import numpy as np

import matplotlib.pyplot as plt
import matplotlib.colors as colors


def load_gis(z, x1, x2, y1, y2):

    a=np.array([])
    for x in range(x1, x2+1):
        b=np.array([])
        for y in range(y1, y2+1):
            print  x, y

            #地形データを読み込む
            url = 'http://cyberjapandata.gsi.go.jp/xyz/dem/{0}/{1}/{2}.txt'.format(z,x,y)
            print url
            response = requests.get(url)
            if response.status_code == 404:
                 Z = np.zeros((256, 256))
            else:
                 #標高値がない区画はとりあえず0mに置換する
                 maptxt = string.replace(response.text, u'e', u'0.0')
                 Z = pd.read_csv(StringIO(maptxt), header=None)
                 Z = Z.values
                 #Z = Z.transpose() 
            print Z
            
            if len(b)==0:
                b = Z
            else:
                b = np.append(b,Z,0)
            
    
        if len(a)==0:
            a = b
        else:
            a = np.append(a,b,1)
            
    return a
    



z=12
x1=3624
x2=3628
y1=1615
y2=1619

        
        
Z = load_gis(z, x1, x2, y1, y2)        
print Z


# 勾配を求める
 (Zy,Zx) = np.gradient(Z)
    
# 勾配を図にしてみる
fig = plt.figure(figsize=(12, 12), dpi=72, facecolor='w', edgecolor='k')
plt.imshow(-Zx, cmap="binary_r")
#plt.savefig('1.jpg', dpi=72)
plt.show()

plt.cla()
fig = plt.figure(figsize=(12, 12), dpi=72, facecolor='w', edgecolor='k')
plt.imshow(-Zy, cmap="binary_r")
#plt.savefig('2.jpg', dpi=72)
plt.show()

Z = Zx + Zy
plt.cla()
fig = plt.figure(figsize=(12, 12), dpi=72, facecolor='w', edgecolor='k')
plt.imshow(-Z/2., cmap="binary_r")
#plt.savefig('3.jpg', dpi=72)
plt.show()