メッシュデータの扱い

import pygmt
import numpy as np
import matplotlib.pyplot as plt

等間隔メッシュデータ

\(z(x,y)\) のような形式で，\(x\), \(y\)が等間隔のメッシュで与えられているデータを可視化してみましょう．

まず，サンプルとして NumPy で適当な等間隔データを作成してみます．

nx = 500
ny = 500
x = np.linspace(0, 2, nx+1)
y = np.linspace(0, 2, ny+1)

z = np.zeros([nx+1, ny+1])

for i in range(nx):
    for j in range(ny):
        a = np.pi * ( x[i]    - np.sqrt(3.0) * y[j] )
        b = np.pi * ( 3.*x[i] + np.sqrt(3.0) * y[j] )
        z[i,j] = np.cos(2*a) * np.cos(4*b) + np.cos(5*a) * np.cos(3*b) + np.cos(7*a) * np.cos(b)

試しにmatplotlibでプロットしてみると，以下のような感じです．

g = plt.subplot()
g.pcolormesh(x, y, z, shading='auto')
g.set_aspect('equal')

このデータを pygmt の xyz2grd を用いて grdデータにしてプロットしてみます．

Note

2023年6月現在，xyz2grd はx, y, zそれぞれ1次元配列にしないと渡せません．そのため，np.meshgrid()と ravel() メソッドを通じて nx * ny サイズの1次元配列を作っています．

もうすこしスッキリした書き方としては，x, y, zの代わりにそれらがパッキングされたdataリストを内包表記を用いて以下のように与える方法もあります．

pygmt.xyz2grd(
    data = [[x[i], y[j], z[i,j]] for i in range(len(x)) for j in range(len(y))]
    # 他の要素は省略
    )

fig = pygmt.Figure()

X, Y = np.meshgrid(x, y)

grddata = pygmt.xyz2grd(
    region = [0, 2, 0, 2],     
    spacing = '0.004/0.004', 
    x = X.ravel(),
    y = Y.ravel(), 
    z = z.ravel()
)

pygmt.makecpt(
    cmap       = 'viridis',     # 比較のためカラーパレットはほとんど同じものを使う
    series     = [-3, 3, 1], 
    continuous = True
)

fig.grdimage(
    grid       = grddata, 
    projection = 'X10c/10c', 
    frame      = ['WSen', 'xaf+l"x"', 'yaf+l"y"'],    
)

fig.grdcontour(
    grid     = grddata, 
    interval = 1, 
    pen      = 'default,black'
)

fig.colorbar(
    position = '+e'    
)

fig.show()

ここでは，grdimageのほかにgrdcontourも用いて等値線も描画してみました．使い方は上記の例を見ればほとんど明らかだと思います．

2次元データの補間と可視化

続いて，粗いデータを補間しつつプロットしてみましょう．

サンプルデータとして，Wikipediaの双3次補間の解説にある例を採用します．ここからコードをお借りして（Creditはリンク先参照；CC-BY），以下のようなデータを使います．まずはデータ生成とmatplolibでの標準的可視化の例を示します．

METHODS = [ 'nearest', 'bilinear' ]
COLORS  = 'viridis'

N = 5
np.random.seed(1)
grid = np.arange(0, N, 1)
data = np.round(np.random.rand(N, N), 1)
mesh = np.meshgrid(grid, grid)

for interp in METHODS:

    fig = plt.figure(figsize=(5,5))

    ax = fig.add_axes([0.125, 0.175, 0.75, 0.75])
    plt.imshow(data, interpolation=interp, cmap=COLORS, vmin=0, vmax=1)
    plt.plot(mesh[0], mesh[1], marker='.', ms=8, color='k', lw=0)
    plt.title(interp, weight='bold')
    plt.xlim(grid.min()-0.5, grid.max()+0.5)
    plt.ylim(grid.min()-0.5, grid.max()+0.5)
    plt.xticks(grid)
    plt.yticks(grid)

    cax = fig.add_axes([0.125, 0.075, 0.75, 0.03])
    cb = plt.colorbar(cax=cax, orientation='horizontal',
                      ticks=np.linspace(0, 1, 6))
    cb.solids.set_edgecolor('face')

ここで作られたデータを ravel()メソッド1次元化し，pygmt.surfaceで補間をしてみます．比較のため pygmt.xyz2grd で粗い間隔のままのプロットも作成します．

fig = pygmt.Figure()

griddata_s = pygmt.surface(
    x       = mesh[0].ravel(), 
    y       = mesh[1].ravel(), 
    z       = data.ravel(), 
    region  = [-0.5, 4.5, -0.5, 4.5], 
    spacing = '0.02/0.02',
    tension = 0.5 # 無指定は0
)

griddata_x = pygmt.xyz2grd(
    x       = mesh[0].ravel(), 
    y       = mesh[1].ravel(), 
    z       = data.ravel(), 
    region  = [-0.5, 4.5, -0.5, 4.5], 
    spacing = '1/1',
    registration = 'p'
)


pygmt.makecpt(
    cmap = 'viridis', 
    series = [0, 1, 0.1], 
    continuous = True
)

fig.grdimage(
    grid = griddata_x, 
    projection = 'X10c/10c',
    frame = ['WSen+t"GMT xyz2grd"', 'xaf', 'yaf'], 
)

fig.plot(
    x = mesh[0].ravel(), 
    y = mesh[1].ravel(), 
    style = 'c0.2c', 
    fill  = 'black', 
)

# 右にずれる
fig.shift_origin( xshift = 11 )

fig.grdimage(
    grid = griddata_s, 
    projection = 'X10c/10c',
    frame = ['WSen+t"GMT surface"', 'xaf', 'yaf']    
)

fig.grdcontour( 
    grid = griddata_s, 
    pen = 'default,black', 
    annotation = '-'
)

fig.plot(
    x = mesh[0].ravel(), 
    y = mesh[1].ravel(), 
    style = 'c0.2c', 
    fill  = 'black', 
)

fig.colorbar(
    position = '+e'    
)

fig.show()

Note

pygmt.surface は tension factor \(t\) という量で補間の状況をコントロールします．その値を指定するオリジナルのGMTの -T オプションは，v0.8.0までが正式には実装されていませんでしたが， pygmt.surface{T = } とすると動作する状態でした．この T はオリジナルのGMTのsurfaceモジュールのオプション名です．

v0.9.0から公式に pygmt.surface に GMTのTオプションに相当するtensionオプションが実装されたようです．