【Python】mpl_toolkitsのBasemapで等高線だけが表示されない原因 - 加賀百万石ですが何か？

【Python】mpl_toolkitsのBasemapで等高線だけが表示されない原因

IT

Pythonを使って等高線を描こうとしていたときに，X, Y, Zのデータを間違いなく準備できているのになぜか等高線が描けないという問題があったのですが，その原因が分かりにくかったので備忘録として残しておきます。

 

現象（等高線が表示されない例）

こういう風にグリッドデータを用意して等高線を描こうとしても何故か等高線だけが描けないという状況でした。

import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

# 上限／下限の設定
lat_min = 20
lat_max = 50
lon_min = 120
lon_max = 180

# 緯度・経度の範囲／単位を取得
lats = np.linspace(lat_min, lat_max, 7)  # 緯度 : 20〜50, 5度毎
lons = np.linspace(lon_min, lon_max, 13) # 経度 : 120〜180, 5度毎

# グリッドデータに変換
x, y = np.meshgrid(lons, lats)
z = np.sqrt(x**2 + y**2)

# Basemapオブジェクトの作成
bmap = Basemap(
    projection='merc',
    resolution='h',
    llcrnrlon=lon_min,
    llcrnrlat=lat_min,
    urcrnrlon=lon_max,
    urcrnrlat=lat_max)

# 等高線を描く
bmap.contourf(x, y, z)
# 陸地の塗り潰し
bmap.fillcontinents(color='lightgray', lake_color='white')
# 緯線を描く
bmap.drawparallels(lats, labels=[1,0,0,0], fontsize=10, color='gray')
# 経線を描く
bmap.drawmeridians(lons, labels=[0,0,0,1], fontsize=10, color='gray')

# 表示
plt.show()

このソースでの表示結果はこうなります。

何故か等高線だけが描かれません。

ちなみにx, y, zのデータはこんな感じで正常に作成されています。

x
array([[120., 125., 130., 135., 140., 145., 150., 155., 160., 165., 170.,
        175., 180.],
       [120., 125., 130., 135., 140., 145., 150., 155., 160., 165., 170.,
        175., 180.],
       [120., 125., 130., 135., 140., 145., 150., 155., 160., 165., 170.,
        175., 180.],
       [120., 125., 130., 135., 140., 145., 150., 155., 160., 165., 170.,
        175., 180.],
       [120., 125., 130., 135., 140., 145., 150., 155., 160., 165., 170.,
        175., 180.],
       [120., 125., 130., 135., 140., 145., 150., 155., 160., 165., 170.,
        175., 180.],
       [120., 125., 130., 135., 140., 145., 150., 155., 160., 165., 170.,
        175., 180.]])

y
array([[20., 20., 20., 20., 20., 20., 20., 20., 20., 20., 20., 20., 20.],
       [25., 25., 25., 25., 25., 25., 25., 25., 25., 25., 25., 25., 25.],
       [30., 30., 30., 30., 30., 30., 30., 30., 30., 30., 30., 30., 30.],
       [35., 35., 35., 35., 35., 35., 35., 35., 35., 35., 35., 35., 35.],
       [40., 40., 40., 40., 40., 40., 40., 40., 40., 40., 40., 40., 40.],
       [45., 45., 45., 45., 45., 45., 45., 45., 45., 45., 45., 45., 45.],
       [50., 50., 50., 50., 50., 50., 50., 50., 50., 50., 50., 50., 50.]])

z
array([[121.65525061, 126.58988901, 131.52946438, 136.47344064,
        141.42135624, 146.37281168, 151.3274595 , 156.28499608,
        161.24515497, 166.20770139, 171.17242769, 176.13914954,
        181.10770276],
       [122.57650672, 127.47548784, 132.38202295, 137.29530218,
        142.21462653, 147.1393897 , 152.06906326, 157.00318468,
        161.94134741, 166.88319268, 171.82840277, 176.7766953 ,
        181.72781845],
       [123.69316877, 128.54960132, 133.41664064, 138.29316686,
        143.17821063, 148.07092895, 152.97058541, 157.87653404,
        162.78820596, 167.70509831, 172.62676502, 177.55280905,
        182.48287591],
       [125.        , 129.80754986, 134.62912018, 139.46325681,
        144.3086969 , 149.1643389 , 154.02921801, 158.90248582,
        163.78339354, 168.67127793, 173.56554958, 178.46568298,
        183.37120821],
       [126.49110641, 131.24404748, 136.01470509, 140.8012784 ,
        145.60219779, 150.41608956, 155.24174696, 160.07810594,
        164.92422502, 169.77926846, 174.64249197, 179.51323071,
        184.39088915],
       [128.16005618, 132.85330256, 137.56816492, 142.30249471,
        147.0544117 , 151.82226451, 156.60459763, 161.40012392,
        166.20770139, 171.02631376, 175.85505395, 180.69311   ,
        185.53975315],
       [130.        , 134.62912018, 139.28388277, 143.96180049,
        148.66068747, 153.3786165 , 158.11388301, 162.86497475,
        167.63054614, 172.4093965 , 177.20045147, 182.00274723,
        186.81541692]])

x.shape
(7, 13)

y.shape
(7, 13)

z.shape
(7, 13)

この状態で数時間ドツボっていたのですが，いざ解決してみるとものすごく単純な話だったようです。

 

解決方法

正しく等高線が描ける全ソースは以下の通りです。

import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

# 上限／下限の設定
lat_min = 20
lat_max = 50
lon_min = 120
lon_max = 180

# 緯度・経度の範囲／単位を取得
lats = np.linspace(lat_min, lat_max, 7)  # 緯度 : 20〜50, 5度毎
lons = np.linspace(lon_min, lon_max, 13) # 経度 : 120〜180, 5度毎

# グリッドデータに変換
x, y = np.meshgrid(lons, lats)
z = np.sqrt(x**2 + y**2)

bmap = Basemap(
    projection='merc',
    resolution='h',
    llcrnrlon=lon_min,
    llcrnrlat=lat_min,
    urcrnrlon=lon_max,
    urcrnrlat=lat_max)

# Basemap用のデータに変換
x, y = bmap(x, y)

# 等高線を描く
bmap.contourf(x, y, z)
# 陸地の塗り潰し
bmap.fillcontinents(color='lightgray', lake_color='white')
# 緯線を描く
bmap.drawparallels(lats, labels=[1,0,0,0], fontsize=10, color='gray')
# 経線を描く
bmap.drawmeridians(lons, labels=[0,0,0,1], fontsize=10, color='gray')

# 表示
plt.show()

結果はこんな感じで等高線が描かれています。

最初の例で等高線が描かれなかったのは，正しく動くソースの28行目のBasemap用のデータに変換する処理が抜けていたためです。

# Basemap用のデータに変換
x, y = bmap(x, y)

普通の感覚では緯度・経度とそこの値を単純に配列で用意すれば等高線が描けそうですが，Basemapを使う場合は緯度・経度の座標をBasemap用の値に変換しないといけないようです。

今回はメルカトル図法で描いたのでその処理の必要性がイマイチ分かりませんでしたが，正距方位図法といったような図を描く場合があることを考えると，確かに単純な緯度・経度から専用のデータに変換してやる必要があるのかなという気はしますね。

ちなみに，pyplotのcontour（contourf）で等高線を描くときはこのような処理は必要ありません。そのせいで，余計にドツボっていたという背景もありますが…。

スポンサーリンク