matplotlibでつくるJoy Divisionの "Unknown Pleasures"

アンノウン・プレジャーズ

アンノウン・プレジャーズ(Unkwon Pleasures)は1979年にJoy Divisionによって発表されたアルバムである。若者が持つ不安定な衝動を感じさせる本作品はポストパンクの金字塔であり、その影響力は今なお強い。そんな本作品のジャッケットはその洗練されたデザインも相まって、ポストパンクのアイコンになっている。このジャケットはpeter saville氏によってデザインされたものであり、“The Cambridge Encyclopaedia of Astronomy”に掲載されていたパルサーCP1919の電磁波の画像を白黒反転したものである。

このジャケットに自分も惹かれ、波形を変えたパロディー画像をmatplotlibで作りたくなった。とは言っても、すでに誰かが作っている可能性もあるなと考えて検索をかけてみた。そしたらとなんとmatplotlibの公式がコードを公開していた！！

MATPLOTLIB UNCHAINED

これには驚いたのと同時に作成意欲が少し削られた。しかし自分オリジナルの画像を作りたかったので続けることとした。

画像にしたい波形の取得

今回はパロディとして日本の地形をプロットする。標高を取得するにあたっては、日本の国土地理院が提供する経緯度を渡すと標高を返してくれるAPIを使用した。

import pandas as pd
import numpy as np
import requests

#APIで使用するURL
base_url = "http://cyberjapandata2.gsi.go.jp/general/dem/scripts/getelevation.php?"
data_type_url = "outtype=JSON"

#　取得したい地域の指定 （以下の例は伊豆半島)
leftbottom = [138.622787, 34.579943] #プロットしたい領域の左下の経度と緯度
righttop = [139.197175,35.175773]　#プロットしたい領域の右上の経度と緯度
longitude_n = 100 # 経度の解像度
latitude_n = 100 # 緯度の解像度
x = np.linspace(leftbottom[0],righttop[0],longitude_n) # 経度の点の作成
y = np.linspace(leftbottom[1],righttop[1],latitude_n) # 緯度の点の作成
X, Y = np.meshgrid(x, y)

#　各経緯度の標高の取得
Elevation=np.empty((latitude_n,longitude_n))
for i in range(longitude_n):
    for j in range(latitude_n):
        lon = X[j,i]
        lat = Y[j,i]
        request_information_url = "lon="+str(lon)+"&lat="+str(lat)+"&"+data_type_url
        url = base_url + request_information_url
        json = requests.get(url).json()
        elevation = json["elevation"]
        if elevation == '-----': # 多分、海上の標高はデータがないので'-----'となっている
            Elevation[j][i] = 0
        else:
            Elevation[j][i] = elevation

使用上の注意点(国土地理院HPより引用)

• サーバに過度の負担を与えないでください。過度の負担を与えると判断したアクセスについて、国土地理院は予告なく遮断を行う場合があります。
• 国土地理院は、本プログラムを予告なく内容変更したり提供停止したりする可能性があります。

matplotlibで描写する

matplotlib公式がコードを公開しているのでそれをベースに作成してはいるが、細かな点を修正した。

import matplotlib as mpl
import matplotlib.pyplot as plt

#　マージンの指定
mpl.rcParams['axes.xmargin'] = 0.3
mpl.rcParams['axes.ymargin'] = 0.1

fig = plt.figure(figsize=(8, 8), facecolor='black')
ax = plt.subplot(111,frameon=False)

lw = 1.5 #　linewidth
pitch = 150 # plot間のピッチの調整(標高のアスペクト比の調整)
alpha = 1.0 # 見かけの調整
xposi = np.linspace(-1,1,Elevation.shape[-1])

for i in range(Elevation.shape[0]):
    j = Elevation.shape[0]-1-i
    ax.stackplot(xposi, j*pitch+Elevation[j,:]**alpha, edgecolor='white', colors=('black',),lw=lw)

#力技による細かい修正
ax.plot([-1,-1],[0,(latitude_n-1)*pitch+Elevation[latitude_n-1][0]**alpha],color="black", lw=3)
ax.plot([1,1],[0,(latitude_n-1)*pitch+Elevation[latitude_n-1][longitude_n-1]**alpha],color="black", lw=3)
ax.plot([-1,1],[-1.5,-1.5],color="black", lw=3)
ax.set_xlim(-1.8,2)

出力結果

伊豆半島

北アルプス

阿蘇山

以上！！