個人的ノート集トップ
PyGMT
第6章　二次元プロット

第6章　二次元プロット

ここでは２次元平面状の情報をプロットする方法をまとめる．

GMTにおいて2次元のカラーメッシュ図や等高線図のプロットは， グリッドデータを扱うことで実現する． ここでいうグリッドデータとは，緯度経度単位で区切られたグリッドの各々に，標高などのデータが入ったデータのことである． グリッドの間隔は緯度方向と経度方向それぞれで等間隔でなければならないようだ．

PyGMTにおいてグリッドデータは， netCFDフォーマットのグリッドファイルを読み込むか， プログラム中でxarray.DataArray型の変数を作成することで使用できる．

以下ではまずグリッドデータの準備について記し，その後でグリッドデータの二次元カラーメッシュまたは等値線によるプロット方法を記す．

描画用グリッドデータの準備

グリッドデータの型について

PyGMT上でのグリッドデータは緯度，経度を軸に持つ2次元のxarray.DataArrayである． したがって自分でxarrayを使ってグリッドデータを作成することが可能．

標高グリッドのダウンロード

次のようにしてデフォルトのデータセットを入手できる：

grid_data = pygmt.datasets.load_earth_relief( 
                    resolution = '01m', 
                    region = insrgn
                    )

オンラインでなければならない．

grdファイルの代わりにxarray.DataArrayを利用

import xarray as xr
<data> = xr.Datarray(
data = <2D ndarray>,
dims = ['<index name 1>', '<index name 2>', ……],
cords = {
        '<index name 1>': <ndarray 1>,
        '<index name 2>': <ndarray 2>,
        ……
    }

２次元以上のグリッド状のデータがndarray形式であっても上記のように座標（index）の情報を追加してxarrayに変換すればプロットすることができる．デフォルトではxarrayがインストールされていない可能性があるので，その場合はxarrayのインストールが必要．

ここで作成したxarray形式の変数をgrdimage()のgridに与えればプロットできる．

xyzファイルからグリッドデータへの変換

ここでいうxyzデータとは，名前の通りx,y,z座標が並べられたテキストデータである．各グリッド（点）は改行で区切られている． xyzデータで標高データなどを与えられた場合，それをgrdデータに変換することで，前章と同様に等値線などでプロットできる．次のコマンドで変換できる：

grid_data = pygmt.xyz2grd(
     data =’<~.xyz>', [x=<lat>,y=<lon>, z=<value>でも可]
     region = <region>,
     spacing = <e.g. '15s'>
)

このコマンドでgrdデータとして扱うことができるようになる．このgrid_dataをgrdcontour()に与えればよい．

nanの部分の処理

地理院からダウンロードした5 mメッシュの標高地図などで，値がない(nan)メッシュが含まれることがある． その部分の値をnodata=‘0’のように任意に指定することができる． 指定しなければnanのままである．

なおnanのままにすればgrdimageにおいてnanの部分を透過させることができる．

xyz2grdのパラメータにおける注意

spacingとregionの指定は必須だが，両者の値の組み合わせに条件があることに注意． regionはspacingで指定した長さで割り切れる必要がある． なので適当に指定するとまずエラーが出る． region=[xs,xe,ys,ye]であるが，spacingをregionの幅xe-xs, ye-ysを整数で割った値にすると確実．

次の例では，あらかじめspacing用にdx, dyを計算して用いている． 下記の例で示すdx, dyよりも細かいspacingを指定すると，きれいに充填されたグリッドにならず隙間だらけに値が入ったグリッドになる．

xs = min(dat[:,1])
xe = max(dat[:,1])
ys = min(dat[:,0])
ye = max(dat[:,0])
dx = (xe-xs)/nx
dy = (ye-ys)/ny
dx = str(dx)
dy = str(dy)
（略）
grd = pygmt.xyz2grd(
    x = dat[:,1],
    y = dat[:,0],
    z = dat[:,2],
    region = [xs,xe,ys,ye],
    spacing = dx+'/'+dy,
    nodata = '0',
    )
xyz2grdを用いて変換したデータをプロットした例．適切な設定を行わなければこのように正常にプロットできない．(標高データ: 国土地理院基盤地図情報数値標高モデル)

xyz2grdを用いて変換したデータをプロットした例．適切な設定を行わなければこのように正常にプロットできない．(標高データ: 国土地理院基盤地図情報数値標高モデル)

グリッドデータのプロット

等値線図

(追記v0.12.0)

fig.grdcontour(
    grid = <グリッドファイル名/グリッドデータ変数>,
    annotation = '<アノテーション (計曲線) 間隔>+f<フォント設定>',
    label_placement = '<アノテーションの配置設定>',
    levels = <等値線 (主曲線) 間隔>,
    limit = [<下限>, <上限>]    # 等値線描画間隔
    pen = [
        'a<0.5p,black>',    # 計曲線の線種
        'c0.1p,grey80',    # 主曲線の線種
        ],
)

引数gridにグリッドデータをファイル名または変数で与える．

annotationのオプション（抜粋）:
アノテーションの文字の向きなどの設定を行う．
(参考:GMT公式マニュアル https://docs.generic-mapping-tools.org/dev/grdcontour.html)

  • +f: アノテーションのフォント設定を行う．
    「f10p」のようにして文字サイズを調整するなど．
  • +a: 等値線にたいするアノテーションの角度 (angle) を指定する．
    「+an」は等値線に直交，「+ap」で等値線と平行とする． 「+ap」の場合はさらに「u」または「d」を後ろに加えると，それぞれ勾配の上側，下側にアノテーションの上下を向けることができる．
  • +w: アノテーションの角度を決める際に等値線を構成する点を何点使用するかを定める．
    適切に設定することでアノテーションの角度が自然になる． 少なすぎると等値線の局所的な角度に従ってしまい，感覚に合わない角度になりうる．
  • +r: アノテーションをプロットする最小の曲率を指定する．
    急すぎるカーブにアノテーションを配置しないようにできる．
  • +u: アノテーションの数値の直後に加える単位を指定する．
  • +v: 単に「+v」と加えると，アノテーションの文字を等値線に沿わせることができる．
    しかしよほどなめらかな等値線でないとうまくいかないように感じる．
  • +i: ラインを非表示 (invisible) にし，アノテーションだけを描画する．

label_placementのオプション (抜粋):
等値線に対するアノテーションラベルの位置や数を指定する．
（参考: GMT公式マニュアル https://docs.generic-mapping-tools.org/6.0/cookbook/contour_annotations.html）

  • d: 「d<間隔>」で等値線に沿った間隔を指定することでアノテーションの位置を調整する．
    間隔にはc,i,p（センチメートル，インチ，ポイント）で単位を与えることができる．（例: d5c）
  • D: dオプションとほぼ同じだが，単位に紙面上の単位ではなく地図内の単位を使用できる．
    単位はd,m,s,e,f,k,M,n（度，分，秒，メートル，フィート，キロメートル，法定マイル，海里）.
  • f: 「f<ファイル名>」の形で与えたファイルに従ってアノテーションを配置する．
    ファイルはテキストファイルで，１列目と2列目にアノテーションを配置する座標を与える． 実際の等値線がこのファイル名の座標を通る場合にのみアノテーションが配置される． 厳密に等値線が通る座標を指定することは困難な場合が多いであろうから，「f<ファイル名>/<slop>」とスラッシュで区切ってslopに許容誤差を与えることができる．単位はc,i,pから選択可．
  • l: 「l<始点x>/<始点y>/<終点x>/<終点y>」と記すことで，始点と終点を結ぶ紙面上の線分と等値線が交わる位置にアノテーションを配置する．
    文字列内で「,」で区切ることで複数の線分を与えることもできる．
  • L: lオプションとほぼ同じだが，線分を大円経路で考える．
  • n: 「n<アノテーション数>」各等値線上のアノテーションの数を指定する．
    「n1」とすれば，各等値線 (計曲線) に1つずつアノテーションが追加される． 「n<アノテーション数>/<最少間隔>」とすればアノテーションが詰まりすぎないように最少間隔もc,i,pのいずれかの単位で指定できる．
  • N: nオプションとほぼ同じだが，アノテーションが等分された等値線の各区間の終点に配置される．

levelsやlimitをプロットするデータの値域やregionの領域に合わせないと，

grdcontour [WARNING]: No contours within specified -L

というエラーが出てプロットできないので注意．

（参照）
https://www.pygmt.org/latest/tutorials/advanced/contour_map.html#sphx-glr-tutorials-advanced-contour-map-py

例）等深線プロット．アノテーションラベルのフォントサイズや向きなどを指定している．

import pygmt
prj = 'M5c'
rgn = [140,146,39,43]

# グリッドデータダウンロード
grid_data = pygmt.datasets.load_earth_relief( 
                    resolution = '15s', 
                    region = rgn,
                    )

fig = pygmt.Figure()
# 海岸線のプロットと陸の塗りつぶし
fig.coast(
    shorelines = '0.1p',
    land = 'grey70',
    projection = prj,
    region = rgn,
)

# 等深線のプロット
fig.grdcontour(
    grid = grid_data,
    levels = 500,
    annotation = '2000+f6p+ap+w200+e',
    pen = ['a0.4p,black','c0.2p,grey30,-'],
    limit = [-7000,-1],
    frame = ['neWS','agf'],
    label_placement = 'l140/40/146/40',
)
fig.show()
水深データの等深線プロット．アノテーションの向きや位置，等値線の線のスタイルなどを設定している．

水深データの等深線プロット．アノテーションの向きや位置，等値線の線のスタイルなどを設定している．

カラーメッシュ図

pygmt.grdgradientで傾斜を計算できる．これを用いると後で陰影をつけられる：

gradient_data = pygmt.grdgradient(
    grid = grid_data,
    azimuth = [45,135],
    normalize = 'e0.8')

makecptでカラーパレットを作り，grdimageでプロット．Figure.colorbarでカラーバー：

pygmt.makecpt(
    cmap = 'gebco')

fig.grdimage(
    projection = prj,
    region = rgn,
    grid = grid_data,
    shading = gradient_data
    )

fig.colorbar(
    frame = ['','y+lm'],
    position='JMR+o1c/0c'
    #J: 枠の外，MR:middle right, o1c/0c:右に1cm，上に0cmに描画
    )

カラーパレットはGMTのページで確認できる．

grdimageではnan_transparent=Trueとすると，グリッドのnanになっているところを透過させられる．

例) カラーメッシュ図のプロット

import pygmt
prj = 'M5c'
rgn = [139,146,41,46]

# 標高グリッドデータダウンロード
grid_data = pygmt.datasets.load_earth_relief( 
                    resolution = '30s', 
                    region = rgn,
                    )

fig = pygmt.Figure()
with pygmt.config(
    FONT = '8p',
):
    # 標高プロット
    fig.grdimage(
        grid = 'earth_relief.nc',
        cmap = 'geo',
        projection = prj,
        region = rgn,
        frame = [
            'neWS',
            'agf',
            ],
    )
# 画像高さ取得
with pygmt.clib.Session() as session:
    with pygmt.helpers.GMTTempFile() as tmp:
        session.call_module("mapproject",f"-Wh ->{tmp.name}")
        h = tmp.read()

# 画像の高さと同じ長さのカラーバーを追加
fig.colorbar(
    frame = ['x+lElevation [m]'],
    position = f'JMR+w{h}c',
)
fig.savefig(
    '../fig/grdimage.png',
)
fig.show()
標高グリッドのgrdimageによるプロット．

標高グリッドのgrdimageによるプロット．