はじめに
高等学校の地学の教科書や図説には、日本付近の地形が分かりやすい立体段彩図が載っている。そんな立体段彩図をPyGMTで描いてみたいと思っていたが、なかなかうまくできなかった。最近になって、やっとある程度満足できるものが描けるようになったので、ここでその手順を紹介する。
立体段彩図を描くことが目的なので、PyGMTのコマンドおよびオプションの説明は、プログラムの中にコメントで記述する程度にとどめている。
海岸線を描く
まずは簡単な日本付近の海岸線を描いてみる。メルカトール図法を用いて、北緯28°～50°、東経128°～150°の範囲を描く。
PyCharmを起動し、メニューの①ファイルから➁新規を選択する。
新規ウィンドウが開いたら、Pythonファイルを選ぶ。
新規Pythonファイルのウィンドウで、ファイル名を入力する。ここではpygmt01とする。拡張子.pyは自動的に付けられるので、入力する必要はない。
次のコードをコピペ。
import pygmt
# 地図を描く範囲の設定
reg = [128, 150, 28, 50]
fig = pygmt.Figure()
# 基図を描く
fig.basemap(
region=reg,
# メルカトール図法により横幅を15cmで描く
projection = "M15c",
# 5度ごとに緯度経度線、10度ごとに数値、白黒パターンの枠を描く
frame = ["g5a10f"]
)
# 海岸線を描く
fig.coast(
# 海岸線のみを0.5ポイントの黒実線で描く
shorelines="1/0.5p,black"
)
# png形式でファイルを出力
fig.savefig("japan01.png")①がpygmt01になっていなければ、押して➁現在のファイルに変更する。
③実行をクリック。
プログラムが終了したら、実行ウィンドウが開いてメッセージが表示される。このとき、”④プロセスは終了コード0で終了しました”、になっていれば成功である。
次のような地図が、PyGMTフォルダの中にjapan01.pngというファイル名で作成される。
次に、海と陸を指定した色で塗分ける。
新規にpythonファイルを作成する。ここではpygmt02とする。
次のコードをコピペして実行する。
import pygmt
# 地図を描く範囲の設定
reg = [128, 150, 28, 50]
fig = pygmt.Figure()
# 基図を設定
fig.basemap(
region=reg,
# メルカトール図法により横幅を15cmで描く
projection="M15c",
)
# 陸と海を塗分け海岸線を描く
fig.coast(
# 陸をlightgrayで塗る
land="lightgray",
# 水域をlightcyanで塗る
water="lightcyan",
# 海岸線のみを0.1ポイントの黒実線で描く
shorelines="1/0.5p,black"
)
fig.basemap(
# 5度ごとに緯度経度線、10度ごとに数値、白黒パターンの枠を描く
frame=["g5a10f"]
)
fig.savefig("japan02.png")次のような地図が、PyGMTフォルダの中にjapan02.pngというファイル名で作成される。
平面段彩図を描く
日本付近の平面段彩図を描く。
必要になる標高のデータは、ネットを通して自動的に得られる。データ間隔を狭くすると、読み込むのに時間がかかるので、適当な値にする。
次のコードをコピペして実行する。
import pygmt
reg = [128, 150, 28, 50]
fig = pygmt.Figure()
# 標高データを取得
grid_data = pygmt.datasets.load_earth_relief(
region=reg,
# 1分間隔でデータを取得
resolution="01m"
)
fig.grdimage(
# メルカトール図法により横幅を15cmで描く
projection = "M15c",
grid = grid_data,
# 影なしの地図
shading = False,
# 10度ごとに数値、白黒パターンの枠を描く
frame=["a10f"],
# カラーマップとしgeoを使用
cmap = "geo"
)
fig.coast(
# 海岸線のみを0.1ポイントの黒実線で描く
shorelines = "1/0.5p,black"
)
fig.savefig("japan03.png")次のような地図がjapan03.pngというファイル名で作成される。
地形の凹凸があまり明瞭ではない。もう少し地形の凹凸が表現されるように陰影をつける。
次のコードをコピペして実行する。
import pygmt
reg = [128, 150, 28, 50]
fig = pygmt.Figure()
# 標高データを取得
grid_data = pygmt.datasets.load_earth_relief(
region=reg,
# 1分間隔でデータを取得
resolution="01m"
)
fig.grdimage(
# メルカトール図法により横幅を15cmで描く
projection = "M15c",
grid = grid_data,
# 影ありの地図
shading = True,
# 10度ごとに数値、白黒パターンの枠を描く
frame=["a10f"],
# カラーマップとしgeoを使用
cmap = "geo"
)
fig.coast(
# 海岸線のみを0.1ポイントの黒実線で描く
shorelines = "1/0.5p,black"
)
fig.savefig("japan03-1.png")次のような地図がjapan03-1.pngというファイル名で作成される。だいぶ見栄えがよくなった。
より凝った段彩図の描き方がTakuto Maedaさんのページに載っていたので、参考にさせて頂いた。
次のコードをコピペして実行する。
import pygmt
reg = [128, 150, 28, 50]
fig = pygmt.Figure()
# 標高データを取得
grid_data = pygmt.datasets.load_earth_relief(
region=reg,
# 1分間隔でデータを取得
resolution="01m"
)
# 斜度データの作成
gradient_data = pygmt.grdgradient(
grid = grid_data,
# 方位
azimuth = [45,145],
normalize = "e0.7"
)
fig.grdimage(
# メルカトール図法により横幅を15cmで描く
projection = "M15c",
grid = grid_data,
# 影ありの地図
shading = gradient_data,
# 10度ごとに数値、白黒パターンの枠を描く
frame=["a10f"],
# カラーマップとしgeoを使用
cmap = "geo"
)
# 海岸線描画と海域を半透明シェーディング
fig.coast(
area_thresh = '100',
water = 'black@80',
)
# 海岸線の再描画
fig.coast(
# 海岸線のみを0.1ポイントの黒実線で描き直す
shorelines = "1/0.5p,black"
)
fig.savefig("japan03-2.png")次のような地図がjapan03-2.pngというファイル名で作成される。
2つの地図を比較すると、japan03-2.pngのほうが細かい地形まで表現されている。
違いがよく分かるのは、富山湾から日本海の海底に延びる海底谷だ。japan03-2.pngには、海底谷がよりはっきりと表現されている。
立体段彩図を描く
いよいよ日本付近の立体段彩図を描く。
次のコードをコピペして実行する。
import pygmt
reg = [128, 150, 28, 50]
reg_3d = [128, 150, 28, 50, -10000, 5000]
fig = pygmt.Figure()
# 標高データを取得
grid_data = pygmt.datasets.load_earth_relief(
region=reg,
# 1分間隔でデータを取得
resolution = "01m"
)
fig.grdview (
region=reg_3d,
# メルカトール図法により横幅を15cmで描く
projection="M15c",
grid = grid_data,
# 見る方位と仰角
perspective = [150, 30],
# 高さは2cm
zsize="2c",
# 地形の表面を塗色
surftype="s",
# 10度ごとに数値、白黒パターンの枠を描く
frame=["a10f"],
# カラーマップとしgeoを使用
cmap = "geo",
# -10000m～0mの範囲の壁面をgrayで描く
plane = "-10000+ggray"
)
# 海岸線の再描画
fig.coast(
perspective = [150, 30, 0],
shorelines = "1/0.5p,black"
)
fig.savefig("japan04.png")次のような地図がjapan04.pngというファイル名で作成される。
地形の凹凸が強調されるように陰影をつける。
次のコードをコピペして実行する。
import pygmt
reg = [128, 150, 28, 50]
reg_3d = [128, 150, 28, 50, -10000, 5000]
fig = pygmt.Figure()
# 標高データを取得
grid_data = pygmt.datasets.load_earth_relief(
region=reg,
# 1分間隔でデータを取得
resolution = "01m"
)
fig.grdview (
region=reg_3d,
# メルカトール図法により横幅を15cmで描く
projection="M15c",
grid = grid_data,
# 見る方位と仰角
perspective = [150, 30],
# 影あり
shading=True,
# 高さは2cm
zsize="2c",
# 地形の表面を塗色
surftype="s",
# 10度ごとに数値、白黒パターンの枠を描く
frame=["a10f"],
# カラーマップとしgeoを使用
cmap = "geo",
# -10000m～0mの範囲の壁面をgrayで描く
plane = "-10000+ggray"
)
# 海岸線の再描画
fig.coast(
perspective = [150, 30, 0],
shorelines = "1/0.5p,black"
)
fig.savefig("japan04-1.png")
次のような地図がjapan04-1.pngというファイル名で作成される。
