静かなる名辞

pythonとプログラミングのこと


xyzの点データを内挿してmeshgridにしmatplotlibでプロットする

はじめに

 pythonでmatplotlibを使って作図するとき、三次元のデータでpcolormeshとかcontourでやるような等高線プロットを作りたいんだけど、手持ちのデータはxyzが紐付いた点のバラバラな離散データだけ……ということがままあります。

 散布図ならそれでも良いのですが、等高線などをプロットしようとするとmeshgrid的な形式のデータが要ると思います。困ります。

 なので、内挿してなんとかしてみます。

方針

 matplotlibのドキュメントを読んでいると、それらしいものを見つけました。

matplotlib.mlab — Matplotlib 3.1.1 documentation

matplotlib.mlab.griddata(x, y, z, xi, yi, interp='nn')[source]
Deprecated since version 2.2: The griddata function was deprecated in Matplotlib 2.2 and will be removed in 3.1. Use scipy.interpolate.griddata instead.

 scipyのscipy.interpolate.griddataを使え、ということらしいです。

scipy.interpolate.griddata — SciPy v0.18.1 Reference Guide

 ちょっと使い方がわかりづらいですが、仮に

  • xy:shape=(n_samples, 2)のxy座標のデータ
  • z:shape=(n_samples,)のz座標のデータ
  • X, Y:xy空間のmeshgrid

 という変数を置くとすると

Z = griddata(xy, z, (X, Y))

 でZ(zのmeshgrid)が得られる、ということのようです。詳細はドキュメントのサンプルを見て確認してみてください。

スポンサーリンク



実験

 二次元の正規分布を作り、真値と内挿で得られた値をpcolormesh, contourでプロットしました。

 コードを以下に示します。

import numpy as np
from scipy import stats
from scipy import interpolate
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def main():
    # data (samples)
    xy = np.random.uniform(low=-10, high=10, size=(500, 2))
    norm = stats.multivariate_normal(mean=[2.0, 3.0], cov=[[3, 1],[1,3]])
    z = norm.pdf(xy)
    
    # 3d scatter of samples
    fig = plt.figure()
    ax = Axes3D(fig)
    ax.scatter(xy[:,0], xy[:,1], z)
    plt.savefig("scatter_3d.png")

    # mesh x-y
    x = y = np.linspace(-10, 10, 500)
    X, Y = np.meshgrid(x, y)

    # interpolated Z
    i_Z = interpolate.griddata(xy, z, (X, Y))

    # true Z
    t_Z = norm.pdf(np.vstack([X.ravel(), Y.ravel()]).T).reshape(X.shape)

    # plot
    fig, axes = plt.subplots(ncols=2)
    plt.subplots_adjust(wspace=0.4)

    # plot true Z
    im = axes[0].pcolormesh(X, Y, t_Z, cmap="jet")
    plt.colorbar(im, ax=axes[0])
    axes[0].contour(X, Y, t_Z, colors=["black"])
    axes[0].set_title("true")

    # plot interpolation Z
    im = axes[1].pcolormesh(X, Y, i_Z, cmap="jet")
    plt.colorbar(im, ax=axes[1])
    axes[1].contour(X, Y, i_Z, colors=["black"])
    axes[1].set_title("interpolation")

    plt.savefig("result.png")

if __name__ == "__main__":
    main()

 stats.multivariate_normalについてはこちらの記事もご覧くだい。

scipyで確率分布のサンプルと確率密度関数を生成する - 静かなる名辞

 それほど難しいことはやっていないので、詳細はコードを読んでいただければわかると思います。出力される画像を以下に示します。

 まず、参考のために出力したサンプルの三次元散布図です。

scatter_3d.png
scatter_3d.png

 けっこうスカスカな印象を受けます。

 次に、真値と内挿で得られた値をプロットした結果です。

result.png
result.png

 完璧とは言えませんが、そこそこ良さそうな結果が得られています。なお、内挿のグラフの外側の白い部分は内挿できなくてnanになっている領域です。griddataのfill_valueオプションなどで挙動を変えられるので、検討してみてください。

まとめ

 このように内挿してプロットすることができます。いまいちなデータしか手持ちにないとき威力を発揮します。

 ただし、元のデータの性質を歪めている側面があるので、ある程度注意が必要になります。それさえ気をつければ十分使えると思います。

追記

 内挿方法による比較を行いました。よろしければこちらも御覧ください。

scipy.interpolate.griddataの内挿方法による違いを比較 - 静かなる名辞