分散した地点データを用いて地図上に等値図を描く

matplotlibで地図上に等値図を描画するためには、ばらばらに分散した地点ごとに存在するデータを編集する必要がある。その手続きと、Basemapを用いたマップでの描画の練習、覚書として、今回は2017年の3か月ごと、および年間の降水量について日本列島に等値図を描くことを目標とする。

気象データの準備

データは気象庁のHPページから取得。それを編集し、下記のような表にしてcsvで保存した（'precipitation.csv'とする）。1，2列目は観測サイトの名称、3，4列目が緯度経度である。5列目から8列目が2017年における3か月ごとの積算降水量である（単位はmm）。最後の9行目は4つのデータの合計値として年間雨量を出したものである。

prefecture	site	lat	lon	Jan-Mar	Apr-Jun	Jul-Sep	Oct-Dec	Ann
宗谷	稚内	45.415	141.6783333	128	188	282	434.5	1032.5
宗谷	北見枝幸	44.94	142.585	142.5	258.5	364	353.5	1118.5
:	:	:	:	:	:	:	:	:

データ読込

次に、pythonでこのデータから描画に必要なデータの配列を作成する。

データの読み込み（一行目は列ラベルなので[1:]で飛ばす）

import numpy as np

def loaddata():
	data = open('precipitation.csv','r').readlines()
	data = [l.split(',') for l in data]
	lats = [np.float(l[2]) for l in data[1:]]
	lons = [np.float(l[3]) for l in data[1:]]
	JanMar = [np.float(l[4]) for l in data[1:]]
	AprJun = [np.float(l[5]) for l in data[1:]]
	JulSep = [np.float(l[6]) for l in data[1:]]
	OctDec = [np.float(l[7]) for l in data[1:]]
	Ann = [np.float(l[8]) for l in data[1:]]

	return {'lats':lats,'lons':lons,'JanMar':JanMar,'AprJun':AprJun,'JulSep':JulSep,'OctDec':OctDec,'Ann':Ann}

Basemap上に観測点のプロット

次に、観測地点をBasemapの白地図上に描画する。

地図上に観測所のプロット

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

lonmin = 120
lonmax = 155
latmin = 20
latmax = 50
gspan = 5

lw = 0.5
glw = 0.3
ms = 3
ticsize = 16
labsize = 18
titsize = 20

def station_map():
	data = loaddata()

	xlab = r'longitude ($^\circ$)'　# '$^\circ$'：〇を上付き文字にして角度の単位
	ylab = r'latitude ($^\circ$)'
	tit = 'Meteorological stations'

# 描画画面の設定
	fig = plt.figure() 
	ax = fig.add_subplot(111)

# Basemapの設定（領域指定、海岸線描画、緯度経度グリッド線描画）
	m = Basemap(resolution='l',llcrnrlat=latmin,llcrnrlon=lonmin,
					urcrnrlat=latmax,urcrnrlon=lonmax,ax=ax)
	m.drawcoastlines(linewidth=lw,color='k')
	m.drawparallels(np.arange(latmin,latmax+1,gspan),linewidth=glw)
	m.drawmeridians(np.arange(lonmin,lonmax+1,gspan),linewidth=glw)

# プロットの描画（ks：黒い（blac'k'）四角（'s'quare））
	m.plot(data['lons'],data['lats'],'ks',ms=ms)

# 軸ラベル・表題の設定
	ax.set_xticks(np.arange(lonmin,lonmax+1,gspan))
	ax.set_xticklabels(map(str,np.arange(lonmin,lonmax+1,gspan)),fontsize=ticsize)
	ax.set_yticks(np.arange(latmin,latmax+1,gspan))
	ax.set_yticklabels(map(str,np.arange(latmin,latmax+1,gspan)),fontsize=ticsize)
	ax.set_xlabel(xlab,fontsize=labsize)
	ax.set_ylabel(ylab,fontsize=labsize)
	ax.set_title(tit,fontsize=titsize)

	plt.show()

全国的に観測地点が散らばっており、南鳥島などの小さな島でも観測が行われていることが分かる。

点群データからメッシュデータへの変換

変換には、matplotlib.mlabライブラリの関数griddata()を用いる。近傍の点同士の間で値が線形的に変化すると仮定して面的にデータを拡張する。メッシュデータにするためには事前にメッシュサイズとその領域を用意しておく必要がある。

メッシュデータの作成

from matplotlib.mlab import griddata

cspan = 0.1

def makemesh(season): #seasonは'JanMar'等の文字列
	data = loaddata()

# 作成するメッシュ範囲と解像度を定義し、メッシュ枠組の作成
	lonmesh, latmesh = np.meshgrid(np.arange(lonmin,lonmax+1,cspan), np.arange(latmin,latmax+1,cspan))

# 関数griddata()を用いて、点群データをメッシュデータへ変換（線形内挿）
	return griddata(data['lons'],data['lats'], data[season], lonmesh, latmesh, interp='linear')

メッシュデータを用いて等値図の描画

作成したメッシュデータを用いて、先ほどの観測所のマップに等値図を加える。
年間雨量含め5種類のマップを作成することができる（画像保存の際はplt.show()の部分を書き換える）。

等値図の描画

import matplotlib.cm as cm

vminmax = (0,1000) #[mm/3か月]
annvminmax = (0,2500) #[mm/年]

def makecontour():
	seasons = ('JanMar','AprJun','JulSep','OctDec','Ann')

	xlab = r'longitude ($^\circ$)' # '$^\circ$'：〇を上付き文字にして角度の単位
	ylab = r'latitude ($^\circ$)'

# 季節ごとの積算雨量で描画するためのforループ
	for season in seasons:
		tit = '%s precipitation' %season

# メッシュデータの作成
		lonmesh, latmesh = np.meshgrid(np.arange(lonmin,lonmax+1,cspan), np.arange(latmin,latmax+1,cspan))
		precmesh = makemesh(season)

# 描画画面の設定
		fig = plt.figure()
		ax = fig.add_subplot(111)

# Basemapの設定（領域指定、海岸線描画、緯度経度グリッド線描画）
		m = Basemap(resolution='l',llcrnrlat=latmin,llcrnrlon=lonmin,
						urcrnrlat=latmax,urcrnrlon=lonmax,ax=ax)
		m.drawcoastlines(linewidth=lw,color='k')
		m.drawparallels(np.arange(latmin,latmax+1,gspan),linewidth=glw)
		m.drawmeridians(np.arange(lonmin,lonmax+1,gspan),linewidth=glw)

# 等値図（contourf）の描画
# levelsは、等値図における最大最小値によらず同じカラーバーを作成するために必要（256段階で色が変化）
# 上限を9999でなくvmaxにすると、vmaxよりも高い値が非表示になる
# extend = 'max'でカラーバーの最大値を三角形に
# cm.jet：青から赤のグラデーション
		vmin,vmax = annvminmax if season == 'Ann' else vminmax
		cf = m.contourf(lonmesh, latmesh, precmesh, 
				levels = np.append(np.linspace(vminmax[0],vminmax[1],255),9999),
				extend = 'max', vmin=vmin, vmax=vmax,cmap = cm.jet)

# カラーバーの作成（Basemap内のどの等値図と対応させるか（はじめのcf）を忘れずに）
		cb = m.colorbar(cf, ticks = range(vminmax[0],vminmax[1]+1,250), location='right', size='3%')
		cb.ax.set_ylabel('[mm]', fontsize=labsize)

# 観測所の描画
		m.plot(data['lons'],data['lats'],'ks',ms=ms)

# 軸ラベル・表題の設定
		ax.set_xticks(np.arange(lonmin,lonmax+1,gspan)) 
		ax.set_xticklabels(map(str,np.arange(lonmin,lonmax+1,gspan)),fontsize=ticsize)
		ax.set_yticks(np.arange(latmin,latmax+1,gspan))
		ax.set_yticklabels(map(str,np.arange(latmin,latmax+1,gspan)),fontsize=ticsize)
		ax.set_xlabel(xlab,fontsize=labsize)
		ax.set_ylabel(ylab,fontsize=labsize)
		ax.set_title(tit,fontsize=titsize)

		plt.show()

代表として、年間雨量の等値図を示す。等値図は点群データすべてが含まれる多角形内でのみ描かれる。太平洋側では、島での高い雨量にその周辺の雨量の値が引っ張られていることが分かる。陸域にだけ着目すると、なんとなく現実的な等値図が描かれているように見える。

海面にマスクをかける

mpl_toolkits.basemapライブラリの関数maskoceans()を用いて、等値図の描画範囲を陸域の限定する。上のコードとほとんど同じであるが、一部追加と修正を加える。

等値図に海面マスクをかける

from mpl_toolkits.basemap import maskoceans

def makecontour2():
:  :  :  :  :  
# 追加：陸や海に着色（0.8は灰色の濃さ）
		m.drawlsmask(land_color='white', ocean_color='0.8', lakes=False)

# 修正：等値図に海面マスク（maskoceans）をかける
		cf = m.contourf(lonmesh, latmesh, maskoceans(lonmesh, latmesh, precmesh),
				levels = np.append(np.linspace(vmin,vmax,255),9999),
				extend = 'max', vmin=vmin, vmax=vmax,cmap = cm.jet)
:  :  :  :  :

再び年間雨量での描画結果を示す。日本海側や九州での雨量が高く、内陸部や北海道での雨量が低いことが分かる。

3か月ごとの積算降水等値図の描画

3か月ごとの降水量については、一つのfigureに4つのaxesを作成して描画した。カラーバーの作成は省略したが、これを行うためには別にaxisを作成し、そこに等値図と同条件となるようなカラーバーを再作成する必要がある。

一つの図に集約（各ウィンドウの軸を重ねる）

def makecontour3():
	data = loaddata()

	lonmesh, latmesh = np.meshgrid(np.arange(lonmin,lonmax+1,cspan), np.arange(latmin,latmax+1,cspan))

# 全体の描画画面の設定（forループの手前で）
	fig = plt.figure(figsize=(10,10))

	for (isea,season) in enumerate(seasons): #enumerate()でisea（数字id）を利用
		if season == 'Ann': continue # 年間雨量はパス

# メッシュデータの作成
		lonmesh, latmesh = np.meshgrid(np.arange(lonmin,lonmax+1,cspan), np.arange(latmin,latmax+1,cspan))
		precmesh = makemesh(season)

# 季節ごとの対象描画画面の設定（2×2の4領域をiseaで指定）
# add_axes([左端x,左端y,全長x,全長y])。いずれもfig画面内の比率で指定。
# 今回は恣意的に縦横比が30(度):35(度)となるよう設定。
		ax = fig.add_axes([0.05+0.42*(isea%2),0.05+0.36*(1-isea/2),0.42,0.36])

:  :  :  :  :

# 軸ラベルの非表示
		ax.set_xticks([])
		ax.set_yticks([])

# 表題（配置場所の自由度が小さい）の代わりに文字列の挿入
		ax.text(lonmin+gspan, latmax-gspan, season, fontsize = titsize)

	plt.show()

色は0mm（青）～1000mm（赤）である。夏場の降水量が全国的に多いほか、中部から近畿地方にかけては10月～12月の雨量も多いことが分かる。