研究室提供のチュートリアル

Python総合講座

本講座は2020年まで公開していた「Anti-Pythonian教授によるPython作図講座」をリニューアルし、ある程度、系統的にPythonコードの書き方を解説するものです。

Ⅰ.気象データの出入力

(1) 所謂、pure binaryデータの読み込みと書き込みです。データはここからdownloadしてください。気温の気候値です。 このデータはlittle endianですが、big endianの場合はmode=">f4"とします。

def read_binary_data():
    import numpy as np
    nx, ny, nz, nt, nm = 288, 145, 37, 12, 6
    t = np.memmap("anl_p125_tmp.bin",      # ファイル名
                  dtype = "f4",            # 4バイト浮動小数点型
                  mode = "r",              # "r"は読み込み専用
                  shape = (nz,ny,nx),      # アレイを指定
                  offset = nm*nz*ny*nx*4)  # 読み込み開始位置のバイト数
    return t
t = read_binary_data()

これに続けて、下記を実行すると、バイナリデータとして書き込まれます。

def write_binary_data(t):
    import numpy as np
    nx, ny, nz, nt, nm = 288, 145, 37, 12, 6
    np.memmap("anl_p125_tmp_Jun.bin",  # ファイル名
              dtype = "f4",            # 4バイト浮動小数点型
              mode = "w+",             # "w+"は書き込みのみ
              shape = (nz,ny,nx))[:,:,:] = t[:,:,:] # 左辺の[:,:,:]は必須
    return 
write_binary_data(t)

(2) Python「専用」形式への読み書きです。npyは他言語からでも読めますが、pythonでしか見たことがないです。アレイの形が保存されるなど便利なので、 データの仮置きにお薦めです。下記、例では(1)で定義したtを使います。

import numpy as np
np.save("anl_p125_tmp.npy",t)         # ファイルにデータを書く

これに続けて、下記を実行すると、アレイの形が保存されていることがわかります。

import numpy as np
t = np.load("anl_p125_tmp.npy")       # ファイルからデータを読む
print(t.shape)

(3) Excelデータの読み書きです。pandasを使うのが正解です。AMeDASデータをエクセルにしたものCSVにしたものを利用します。 pandasは縦がindex横がcolumnsと呼びます。下図のようなイメージです。

import pandas as pd   # エクセル操作関係はpandasを使う
df = pd.read_excel("Amedas.47412.202001.xlsx", index_col=0) # index_col = 0を書かないとUnnamed: 0というコラムが出てしまう
# df = pd.read_csv("Amedas.47412.202001.csv", index_col=0) # csvも同様に読むことができる
display(df) # displayがエラーになる場合があるらしい

上記に続けて、同じ内容を別のエクセルファイルとして書き出します。

import pandas as pd   # エクセル操作関係はpandasを使う
df.to_excel("test.xlsx") # excelはto_excel
df.to_csv("test.csv")    # csvはto_csvで書き込む

気温の高い順のような整列も簡単にできます。

df1 = df.sort_values("日平均気温") # 昇順(気温が低い方から高い方へ並べる)
df2 = df.sort_values("日平均気温", ascending=False) # 昇順(気温が高い方から低い方へ並べる)

行(index)方向の切り出しと列(column)方向の切り出し。

temp = df["日平均気温"].values # .valuesと書くとDataFrameからアレイに変換できる
print(temp)
day1 = df.iloc[0].values # .iloc[行数]が便利
day1_3 = df.iloc[0:3,0:10].values # .ilocなら[行範囲,列範囲]で範囲指定して切り出すこともできる。
print(df.index)  # indexの値を確認
day1 = df.loc[1] # .ilocに代えてloc[index名]でも可

降水データのように観測ナシに--記号が入っている場合に、欠損値に変換する必要があります。

import numpy as np
undef = -999.0  # 欠損値を設定
crain = df["日降水量"].values # 文字型で読まれている
nrain = len(crain)           # アレイの大きさをnrainとする
rain = np.ones((nrain),dtype="f4") * undef # nrainの大きさの浮動小数点数のアレイを用意し、
                                           # np.onesで1をすべていれておいて、欠損値を乗じて、
                                           # すべて欠損値にしておく
for irain in range(nrain):
    try:
        rain[irain] = float(crain[irain])  # floatに変換できない場合は欠損値のままとなる例外処理
    except:
        pass

風向データを角度に変換し、さらに東西風・南北風に変換するプログラムです。

def wind_to_vector(wdir,wspeed):
    import numpy as np
    undef = -999.9
    cdir = ["北","北北東","北東","東北東","東","東南東","南東","南南東",
            "南","南南西","南西","西南西","西","西北西","北西","北北西"]
    if wdir in cdir:
        angle = cdir.index(wdir) * 22.5        # wdirと完全一致するものを見つけます
        u = -np.sin(angle*np.pi/180) * wspeed  # 北から時計まわりの角度であり
        v = -np.cos(angle*np.pi/180) * wspeed  # 風向が風が「吹き込む」方角であることに注意
    else:                           # wdirが風向を表す文字列ではないときはundefを返す
        u, v = undef, undef
    return u, v
wind_to_vector("西北西",3.0)
wind_to_vector("--",3.0)

(4) 気象データによくあるgrib(グリブ)形式の読み書きです。気象庁再解析データJRA55もgribで提供されています。 この例はJRA55に特化したものです。

def inquire_grib_data(path):
    import pygrib                  # gribファイルの中身を見たい場合はinquire_grib_data(path)を実行                
    grbs = pygrib.open(path)
    for grb in grbs:
        print(grb)
    return
def read_grib_data(path,name=None,level=None):
    import numpy as np
    import pygrib                   # pygribは!pip3 install pygrib --userでインストール
    grbs = pygrib.open(path)

    if name != None:                # anl_surf125に対しては変数名を与える
        alines = grbs.select(name=name)
    elif level != None:             # anl_p125に対しては気圧面を与えるとその水平面データ
        alines = grbs.select(level=level)
    else:                           #                  気圧面を与えないと全3次元データ
        alines = grbs.select()

    lat, lon = alines[0].latlons()  # lonは経度、latは緯度データ: (ny,nx)の2次元格子です
    ny, nx = lat.shape
    nline = len(alines)
    gdata = np.empty( (nline,ny,nx), dtype = "f4" )
    levels = np.empty( (nline), dtype = "f4" )
    for iline, aline in enumerate(alines):
        gdata[iline,:,:] = aline.values[::-1,:]
        levels[iline] = aline["level"]

    return lon, lat[::-1], level, gdata 
inquire_grib_data("anl_surf125.2020010100") # これを実行すると、この下のnameに何を与えるかがわかります
lon, lat, _, SLP = read_grib_data("anl_surf125.2020010100",name="Mean sea level pressure")
T850 = read_grib_data("anl_p125_tmp.2020010100",level=850)[3]

(5) 気象データによくあるnetcdf(ネットシーディーエフ)形式の読み書きです。 この例で利用するデータは降水の月平均データGPCPです。

def read_netcdf_data(path):
    import netCDF4                  # netCDF4は!pip3 install netCDF --userでインストール
    import numpy as np
    import datetime as dt           # 日付関係の制御ライブラリ
    nc = netCDF4.Dataset(path, "r")
    var = nc.variables["precip"][:,:,:]  # 格納されているデータ
    lon = nc.variables["lon"][:]         # 経度データ(1次元)
    lat = nc.variables["lat"][:]         # 緯度データ(1次元)
    time = nc.variables["time"][:]       # 時間データ(ただし1800年元日からの日数)
    dtime = [dt.datetime(1800,1,1)+dt.timedelta(days=t) for t in time]
    print(nc.variables.keys())
    print(nc.dimensions.keys())
    nc.close()
    return dtime, lon, lat, var
dtime, lon, lat, prcp = read_netcdf_data("precip.mon.mean.nc")

(6) HTML形式のAMeDASデータの読み書きです。BeautifulSoupを使います。ここでは 2020年1月の札幌の観測値を例に挙げて説明します。まず、HTMLファイルとして、ダウンロードします。過度な自動ダウンロードは先方サーバに負荷を 与えるので、控えてください(つまり下記のプログラムは一回実行すればよく、何度も実行してはいけません)。

def download_amedas_daily(pid,gid,year,month):
    from urllib.request import urlopen
    from bs4 import BeautifulSoup            # !pip3 install bs4 --userでインストール
    path = "https://www.data.jma.go.jp/obd/stats/etrn/view/daily_s1.php?"
    req = "prec_no={0:d}&block_no={1:05d}&year={2:04d}&month={3:02d}&day=&view=p1".format(pid,gid,year,month)
                                             # 文字列.format(変数)で変数に書式を定めて文字列内に反映させる
                                             # {0:d}はpidを整数型で文字にするの意味
                                             # {1:05d}はgidを整数型(zero padding 5桁)で文字にするの意味
    html = urlopen(path+req).read()          
    soup = BeautifulSoup(html, 'html.parser')
    with open("Amedas.{0:d}.{1:04d}{2:02d}.html".format(gid,year,month),"w") as fw:
        fw.write(soup.prettify())
    return
download_amedas_daily(14,47412,2020,1) # AMeDASの都道府県名や地点名は実際にウェブページに行くと確認できます。

次に表形式としてpandasデータフレームを通じて、エクセルに書き出します。AMeDASデータは欠損などあるので、とりあえず文字列として保存するのが安全です。

def read_amedas_daily_html(html_file,excel_file):
    import numpy as np
    import pandas as pd
    from urllib.request import urlopen
    from bs4 import BeautifulSoup            
    html = open( html_file, mode = 'r', encoding = 'utf8').read()
    soup = BeautifulSoup(html, 'html.parser')
    table = soup.find_all("table")[5]
    values = []               # リストだと逐次追加が便利で、サイズ未定のときに使いやすい
    for itd in range(31):
        try:   # 小の月はitdが30のところで下記の試みが失敗することを想定して、ここで回避している。
            values.append([td.get_text().replace("\n","").replace(" ","") \
                           for td in table.find_all("td")[itd*21:(itd+1)*21]]) # \は継続行の意味
        except:
            pass
    values = np.array(values)  # リストのままだと成分表示が面倒なので、アレイにする
    columns=["現地気圧","海面気圧","日降水量","時間最大降水量","10分最大降水量", # AMeDASのウェブページが込み入った表なので
            "日平均気温","日最高気温","日最低気温","日平均湿度","日最低湿度",    # このように作り直した方がきれいな表になる
            "平均風速","最大風速","最大風向","瞬間最大風速","瞬間最大風向",
            "日照時間","降雪量","最大積雪深","天気概況(昼)","天気概況(夜)"]
    df = pd.DataFrame(values[:,1:],index=values[:,0],columns=columns)  # indexは縦、columnsは横
    df.to_excel(excel_file)
    return
read_amedas_daily_html("Amedas.47412.202001.html","Amedas.47412.202001.xlsx")
この結果、下記のようなエクセルファイルが出力されます。

(7) ワイオミング州立大学のサウンディングデータ(アスキーデータ)を読んでみます。これもpandasに格納するのが美しいですかね。

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from urllib.request import urlopen
from bs4 import BeautifulSoup

### この部分は何回も繰り返さないように!
url = "http://weather.uwyo.edu/cgi-bin/sounding?"
req = "region=pac&TYPE=TEXT%3ALIST&YEAR=2021&MONTH=01&FROM=0900&TO=0900&STNM=47412"
soup = BeautifulSoup(urlopen(url+req).read(), 'html.parser')

html = soup.find_all("pre")[0].text.strip().splitlines()  # HTMLを行ごとに文字列リストに格納
print(html)

dat = [[html[1][7*i:7*(i+1)] for i in range(11)]]         # 1行目の文字列をcolumns用に確保
for s in html[4:]:
    dat.append([float(s[7*i:7*(i+1)]) if s[7*i:7*(i+1)].strip()!="" else np.nan for i in range(11)])
    ### 内包表現と呼ばれるもの
    ### X if A else Y はpython的な変な順番だが、もしAが真ならXを実行し、偽ならYを実行するの意味
    ### P for i in range(11)でiを0から10(11じゃない!)までをループする    
df = pd.DataFrame(dat[1:],columns=dat[0])
display(df)

こんな感じの出力となります。2021年1月1日9時、札幌における高層気象観測データです。

	PRES 	HGHT 	TEMP 	DWPT 	RELH 	MIXR 	DRCT 	SKNT 	THTA 	THTE 	THTV
0 	1004.0 	26.0 	-6.5 	-13.5 	58.0 	1.35 	250.0 	7.0 	266.4 	270.2 	266.6
1 	1003.0 	31.0 	-6.5 	-13.7 	56.0 	1.33 	260.0 	4.0 	266.5 	270.2 	266.7
2 	1000.0 	46.0 	-6.3 	-14.3 	53.0 	1.27 	265.0 	4.0 	266.9 	270.5 	267.1
3 	926.0 	642.0 	-11.0 	-17.1 	61.0 	1.09 	290.0 	16.0 	267.9 	271.1 	268.1
4 	925.0 	650.0 	-11.1 	-17.1 	61.0 	1.09 	290.0 	16.0 	267.9 	271.1 	268.1
... 	... 	... 	... 	... 	... 	... 	... 	... 	... 	... 	...
77 	20.0 	25970.0 	-41.1 	NaN 	NaN 	NaN 	205.0 	43.0 	709.6 	NaN 	709.6
78 	17.6 	26834.0 	-43.3 	NaN 	NaN 	NaN 	225.0 	41.0 	729.0 	NaN 	729.0
79 	17.0 	27069.0 	-42.8 	NaN 	NaN 	NaN 	230.0 	41.0 	738.0 	NaN 	738.0
80 	16.0 	27478.0 	-41.8 	NaN 	NaN 	NaN 	250.0 	23.0 	754.0 	NaN 	754.0
81 	15.5 	27693.0 	-41.3 	NaN 	NaN 	NaN 	NaN 	NaN 	762.5 	NaN 	762.5

(8) 気象庁で配信する独自形式nusdasの読み書きです。データ書式の本質は、 このマニュアルの付録Aだけで理解できます。




Ⅱ.地図上への描画関係

(1) 天気図っぽいランベルト正角円錐図法。規準となる緯度経度の情報が不明なので、気象庁の天気図っぽくなるように適当に与えています。 axes.set_extentで指定するようですが、この範囲が何を表しているのかわからず、場当たり的に指定するしかないかもしれません。 figure座標から緯度経度座標への変換とその逆もコードに書きました。図法によって縦横比が決まっているので、後にfigsizeをぴったりはまるように指定し直すと、 無駄な余白がなくなります。

# figure座標: axesで指定した範囲を(0,1)x(0,1)とする座標
# pixel座標: plt.figureで指定した大きさxDPIに合わせ、左下を原点とするpixelで測った座標
# 距離座標: 地図投影における中心を原点とし、距離で測った座標
# 緯度経度座標
### ここではfigure座標から緯度経度座標へと変換しています ###
def transform_figure_to_lonlat(p_on_fig,ax,proj,proj_cart,flush):
    transform = proj_cart._as_mpl_transform(ax) 
    p_on_pix = ax.transAxes.transform(p_on_fig) # figure座標-->pixel座標
    p_on_dis = ax.transData.inverted().transform(p_on_pix)# pixel座標-->距離座標
    p_on_lonlat = proj_cart.transform_point(*p_on_dis, src_crs=proj)# 距離座標-->緯度経度座標 
    if flush:
        print("\nFigure coordinate:(A,B)=({0:6.3f},{1:6.3f})".format(*p_on_fig))
        print("Pixel  coordinate:(W,H)=({0:6.1f}px,{1:6.1f}px)".format(*p_on_pix))
        print("Distance coordinate:(X,Y)=({0:8.1f}km,{1:8.1f}km)".format(*p_on_dis/1000))
        print("Lon/Lat coordinate:(λ,φ)=({0:6.1f}deg,{1:6.1f}deg)".format(*p_on_lonlat))
    return p_on_lonlat, p_on_dis, p_on_pix

### ここでは緯度経度座標からfigure座標へと変換しています ###
def transform_lonlat_to_figure(p_on_lonlat,ax,proj,proj_cart,flush):
    transform = proj_cart._as_mpl_transform(ax) 
    p_on_dis = proj.transform_point(*p_on_lonlat,proj_cart)# 距離座標<--緯度経度座標 
    p_on_pix = ax.transData.transform(p_on_dis)# pixel座標<--距離座標
    p_on_fig = ax.transAxes.inverted().transform(p_on_pix)# figure座標<--pixel座標
    if flush:
        print("\nLon/Lat coordinate:(λ,φ)=({0:6.1f}deg,{1:6.1f}deg)".format(*p_on_lonlat))
        print("Distance coordinate:(X,Y)=({0:8.1f}km,{1:8.1f}km)".format(p_on_dis[0]/1000,p_on_dis[1]/1000))
        print("Pixel  coordinate:(W,H)=({0:6.1f}px,{1:6.1f}px)".format(*p_on_pix))
        print("Figure coordinate:(A,B)=({0:6.3f},{1:6.3f})".format(*p_on_fig))
    return p_on_fig, p_on_pix, p_on_dis

def draw_synoptic_chart(ax_range):   # このルーチンは天気図っぽい地図を書くためのものです。
    import numpy as np
    import matplotlib.pyplot as plt
    import matplotlib.ticker as mticker    
    import cartopy.crs as ccrs             # !pip3 install cartopy --userでインストール
    import cartopy.feature as cfeature     # cartopyはprojが最新でないと入らないなど入れるのは難しい

    proj = ccrs.LambertConformal(central_longitude=140,     # ランベルト正角図法の中心経度
                                 central_latitude=35,       # 同中心緯度
                                 standard_parallels=(25,45))# 同標準緯度2本
    proj_cart = ccrs.PlateCarree()     

    ax = plt.axes(ax_range,projection=proj)               # projection=projで図法指定

# 赤枠の範囲が収まるように指定している。
    ax.set_extent([105,180,0,65],crs=proj_cart)    

# もしも東西5000km, 南北4000kmの範囲を指定したい場合は、上に代えて下記と指定する。
#     ax.set_xlim([-5e6,5e6])
#     ax.set_ylim([-4e6,4e6])    
    
    ax.coastlines(resolution='50m',linewidth=0.2)    # 海岸線を書く
    ax.add_feature(cfeature.LAND.with_scale('50m'), color='#bbbbbb') # 陸を塗潰す
    g = ax.gridlines(crs=proj_cart, linewidth=0.2,color="k")  # 緯経線を書く
    g.xlocator = mticker.FixedLocator(np.arange(-180,190,10)) # 経線は10度ごと(書かない範囲をむやみに指定しない)
    g.ylocator = mticker.FixedLocator(np.arange(-10,80,10))    # 緯線は10度ごと
    
    return ax, proj, proj_cart

def draw_synoptic_chart_example():
    import matplotlib.pyplot as plt   
    fig = plt.figure(figsize=(105/25.4,74.2/25.4),dpi=300)  # figsizeはinch指定なのでmm/25.4する
    ax, proj, proj_cart = draw_synoptic_chart((0,0,1,1))                     
    for p_on_fig in [(0,0),(1,1)]:
        _, _, _ = transform_figure_to_lonlat(p_on_fig,ax,proj,proj_cart,True) # 戻り値不要なら_にするお約束
    for p_on_lonlat in [(120,20),(160,60)]:
        _, _, _ = transform_lonlat_to_figure(p_on_lonlat,ax,proj,proj_cart,True)       
        
    import numpy as np
    ax.plot(np.arange(105,181,1),0*np.ones((76)),color="r",lw=1,transform=proj_cart)
    ax.plot(np.arange(105,181,1),65*np.ones((76)),color="r",lw=1,transform=proj_cart)
    ax.plot(105*np.ones((66)),np.arange(0,66,1),color="r",lw=1,transform=proj_cart)
    ax.plot(180*np.ones((66)),np.arange(0,66,1),color="r",lw=1,transform=proj_cart)
        
    fig.savefig("synoptic_chart.png",bbox_inches="tight",pad_inches=0.1)
    return
draw_synoptic_chart_example()

出力を抜粋します。この出力から図の縦横比は1240.2:876.4であることがわかります。

Figure coordinate:(A,B)=( 0.000, 0.000)
Pixel  coordinate:(W,H)=(   0.0px,   0.0px)
Distance coordinate:(X,Y)=( -4480.8km, -3984.5km)
Lon/Lat coordinate:(λ,φ)=( 106.8deg,  -5.1deg)

Figure coordinate:(A,B)=( 1.000, 1.000)
Pixel  coordinate:(W,H)=(1240.2px, 876.4px)
Distance coordinate:(X,Y)=(  5088.3km,  4069.0km)
Lon/Lat coordinate:(λ,φ)=(-139.7deg,  52.0deg)

Lon/Lat coordinate:(λ,φ)=( 120.0deg,  20.0deg)
Distance coordinate:(X,Y)=( -2117.0km, -1442.2km)
Pixel  coordinate:(W,H)=( 306.3px, 276.7px)
Figure coordinate:(A,B)=( 0.247, 0.316)

Lon/Lat coordinate:(λ,φ)=( 160.0deg,  60.0deg)
Distance coordinate:(X,Y)=(  1219.2km,  2958.4km)
Pixel  coordinate:(W,H)=( 738.7px, 755.5px)
Figure coordinate:(A,B)=( 0.596, 0.862)

完成した図はこのような感じです。

(2) ベクトル図をなるべく等間隔にはっきり描く。axes.quiverは見にくいので、仕方がないからaxes.annotateを「濫用」します。 東西風南北風のデータ(気候値です)をダウンロードしてください。

###  矢印を書くのは千鳥格子(正方形の頂点および対角線の交点からなる模様)のような方が矢印が干渉しずらくてよいので
### この関数ではそれを定義します。
def stagger(m,n):
    import numpy as np
    px, py = [], []
    for kax in [0,0.05]:
        for jax in np.arange(1/float(2*n)+kax,1,1/float(n)):
            for iax in np.arange(1/float(2*m)+kax,1,1/float(m)):
                px.append(iax)
                py.append(jax)
    return px, py

### 指定された格子上の点にベクトルをプロットします。
def plot_vectors(ax,proj,proj_cart,    # cartopyを使って定義したもの(synoptic_chart参照)
                px,py,                # figure座標におけるベクトルを書く場所
                x,y,u,v,              # 緯度経度座標で定義された(x,y)と(u,v)
                factor,               # figure座標でよい塩梅になるように長さを調整する
                ulimit):              # ベクトルを書く長さの下限(これ以下の風速はmask out)
    import numpy as np
    import math
    for iax, jax in zip(px,py):
        p_a = (iax,jax)               # 与えられたfigure座標におけるplot位置
        lon, lat = transform_figure_to_lonlat(p_a,ax,proj,proj_cart,False)[0]
                                      # figure座標の位置を緯度経度に変換
        ix = np.argmin(np.abs(x-lon%360))  # その緯度経度に近いデータ点を探索
        iy = np.argmin(np.abs(y-lat))      # データ点座標indexを(ix,iy)とする
        q_a = transform_lonlat_to_figure((x[ix],y[iy]),ax,proj,proj_cart,False)[0]
                                      # 緯度経度(x[ix],y[iy])点をfigure座標に変換(ほぼp_aと同じはず)
        R1 = transform_figure_to_lonlat((q_a[0]-0.01,q_a[1]),ax,proj,proj_cart,False)[0]
        R2 = transform_figure_to_lonlat((q_a[0]+0.01,q_a[1]),ax,proj,proj_cart,False)[0]
                                      # figure座標点は東西南北の方向が左右上下とは異なるので
                                      # figure座標位置の東西方向が緯度経度座標の東西方向から
                                      # どのくらいずれているかを確認
        angle = math.atan2(R2[1]-R1[1],R2[0]%360-R1[0]%360)  # 上記を角度に変換

        U = u[iy,ix] / 2 * factor
        V = v[iy,ix] / 2 * factor
        uabs = np.sqrt(u[iy,ix]*u[iy,ix]+v[iy,ix]*v[iy,ix])  # 風速の計算

        if uabs >= ulimit:
            U, V = U * np.cos(angle) + V * np.sin(angle), -U * np.sin(angle) + V *np.cos(angle) 
                                         # U, Vを上記で計算したangleで回転変換(U,V=をばらして2行で書いてはいけない)
            x1, y1 = q_a[0] - U, q_a[1] - V 
            x2, y2 = q_a[0] + U, q_a[1] + V                    
            ax.annotate("",xy=(x2,y2),   # ベクトルの終点
                        xytext=(x1,y1),  # ベクトルの始点(annotateは矢印を積極的に書く目的ではないので、こんな仕様に)
                        xycoords="axes fraction", # figure座標で書くの意味
                        arrowprops=dict(arrowstyle="->", lw=0.5, color="k",
						shrinkA=0,shrinkB=0))  # これがないと始点・終点が2ptだけ短くなる       
    return

def plot_unit_vector(ax,proj,proj_cart,
                     px,py,              #参照ベクトルの位置(figure座標)
                     u,v,factor,         #参照ベクトルの大きさ
                     px_text,py_text,    #参照ベクトルに添える注意書きの場所(figure座標, py_textは負値許容)
                     title):             #参照ベクトルに添える注意書き
    u = u/2*factor
    v = v/2*factor
    x1, y1 = px - u, py - v 
    x2, y2 = px + u, py + v    
    ax.annotate("",xy=(x2,y2),xytext=(x1,y1),xycoords="axes fraction",
                 arrowprops=dict(arrowstyle="->", lw=0.5, color="k",shrinkA=0,shrinkB=0))    
    ax.text( px_text, py_text, title, transform=ax.transAxes, ha = "center", va = "center" )    

    return

def plot_vector_example():
    import numpy as np
    import matplotlib.pyplot as plt   
    
    # 東西風と南北風を読み込む
    nx, ny, nz = 288, 145, 37
    x = np.linspace(0,360,nx,endpoint=False)
    y = np.linspace(-90,90,ny,endpoint=True)
    u = np.memmap("anl_p125_ugrd.bin", dtype = "f4", mode = "r", shape = (nz,ny,nx) )[6] 
    v = np.memmap("anl_p125_vgrd.bin", dtype = "f4", mode = "r", shape = (nz,ny,nx) )[6]
    
    fig = plt.figure(figsize=(110/25.4,81.8/25.4),dpi=300)  # figsizeはinch指定なのでmm/25.4する
    ax, proj, proj_cart = draw_synoptic_chart((0,0.05,1,1))
    factor, ulimit = 0.01, 3
    plot_vectors(ax,proj,proj_cart,*stagger(12,8),  # 縦横比1240:876にあわせた千鳥格子でベクトルプロット
                 x,y,u,v,factor,ulimit)
    plot_unit_vector(ax,proj,proj_cart,0.9,-0.05,10,0, #図下に西風10 m/s参照ベクトルをプロット
                     factor,0.8,-0.05,"10 m/s")

    fig.savefig("vector_plot.png",bbox_inches="tight",pad_inches=0.1)

    return
plot_vector_example()

完成した図はこのような感じです。

(3) シェイプファイルを使った都道府県・市町村区境界の描画。 このページからjapan_ver81.shpをダウンロードして解凍したものを同ディレクトリに格納してください。

def draw_ward_border(ax):
    import shapefile     # !pip3 install shapefile --userでインストールしてください
    import cartopy.crs as ccrs

    # shape fileは https://www.esrij.com/products/japan-shp/からダウンロードしましょう
    src = shapefile.Reader("./japan_ver81/japan_ver81.shp",encoding="SHIFT-JIS")

    # 場所にもよりますが、いま書こうとしている石狩振興局はファイルの上100にあります。
    pos_ward = []
    for n in range(100):
        rec = src.record(n)
        if rec[2] == "石狩振興局":
            shp = src.shape(n)   
            pos_ward.append(shp.points)

    # shape file から緯度経度情報を抜き出します 
    for n in range(len(pos_ward)):
        lon = [p[0] for p in pos_ward[n]]
        lat = [p[1] for p in pos_ward[n]]
        
        # 市町村区境界を書きます
        ax.plot( lon,lat, color = "k", linewidth = 0.3, transform=ccrs.PlateCarree() )
        
    return

def draw_ward_border_example():
    import numpy as np
    import matplotlib.pyplot as plt
    import matplotlib.ticker as mticker
    from matplotlib.colors import LogNorm
    import cartopy.crs as ccrs
    import cartopy.feature as cfeature
    
    fig = plt.figure(figsize=(4,4.5),dpi=300)
    lonll, lonur, latll, latur = 141.222, 141.495, 42.955, 43.175
    lonc, latc, lat1, lat2 = 141.3566, 43.0627, 30.0, 60.0
    proj = ccrs.LambertConformal( central_longitude = lonc, 
                                 central_latitude = latc, 
                                 standard_parallels = (lat1, lat2) )

    ax = plt.axes((0,0,1,1),projection=proj)
    ax.set_extent([lonll,lonur,latll,latur])    
    ax.coastlines(resolution='50m',linewidth=0.2)    
    g1 = ax.gridlines(crs=ccrs.PlateCarree(), linewidth=0.2,color="k") 
    g1.xlocator = mticker.FixedLocator(np.arange(140,143,0.1))
    g1.ylocator = mticker.FixedLocator(np.arange(42,44,0.1))
    draw_ward_border(ax)
    fig.savefig("ward_border_example.png",bbox_inches="tight",pad_inches=0.1)
    return
draw_ward_border_example()

完成した図はこのような感じです。これは札幌市周辺です(念のため)。

(4) 緯度経度座標が1次元で与えられたときの描画。 ここからtopo.example.binをダウンロードしてください。

def draw_japan_map(ax_range,map_range):   # このルーチンは日本地図を書くためのものです。
    import numpy as np
    import matplotlib.pyplot as plt
    import matplotlib.ticker as mticker    
    import cartopy.crs as ccrs             # !pip3 install cartopy --userでインストール
    import cartopy.feature as cfeature     # cartopyはprojが最新でないと入らないなど入れるのは難しい

    proj = ccrs.LambertConformal(central_longitude=140,     # ランベルト正角図法の中心経度
                                 central_latitude=35,       # 同中心緯度
                                 standard_parallels=(25,45))# 同標準緯度2本
    proj_cart = ccrs.PlateCarree()     

    ax = plt.axes(ax_range,projection=proj)               # projection=projで図法指定
    ax.set_extent(map_range)    # これが意味不明、左下隅と右上隅の緯度経度ってわけでもない
    
    ax.coastlines(resolution='10m',linewidth=0.2)    # 海岸線を書く
    g = ax.gridlines(crs=proj_cart, linewidth=0.2,color="k")  # 緯経線を書く
    g.xlocator = mticker.FixedLocator(np.arange(-180,190,10)) # 経線は10度ごと(書かない範囲をむやみに指定しない)
    g.ylocator = mticker.FixedLocator(np.arange(-10,80,10))    # 緯線は10度ごと
    
    return ax, proj, proj_cart

def draw_on_synoptic_chart():
    import matplotlib.pyplot as plt   
    fig = plt.figure(figsize=(105/25.4,105.2/25.4),dpi=300)  # figsizeはinch指定なのでmm/25.4する
    ax, proj, proj_cart = draw_japan_map((0.05,0.05,0.90,0.90),[122,148,23,46])              
    nland = 1311
    lon, lat, zs = np.memmap("topo.example.bin",dtype="f4",mode="r",shape=(3,nland))
    C=ax.tricontourf(lon,lat,zs,cmap="terrain",transform=proj_cart)
    plt.colorbar(C,cax=plt.axes((0.05,0.0,0.90,0.02)),orientation="horizontal")
    fig.savefig("synoptic_chart_on.png",bbox_inches="tight",pad_inches=0.1)
    return

draw_on_synoptic_chart()

完成した図はこのような感じです。

(5) 図の「座標変換」。

def convert_fig_xy():  # 日付軸とY軸で構成される図のfigure座標、pixel座標、日付・Y座標の変換
    import numpy as np
    import datetime as dt
    import matplotlib.pyplot as plt
    x = np.array([dt.datetime(1970,1,1,0)+dt.timedelta(days=i) for i in range(20)])
    y = np.linspace(20,60,20)
    fig = plt.figure(figsize=(5.5,3),dpi=300) # 300px/inchなので、X=5.5x300=1650px, Y=3x300=900pxのpixel座標となる。
    ax = plt.axes((0.1,0.1,0.8,0.8))          # figure座標は常に[0,1]x[0,1]
    ax.set_xlim([x[0],x[-1]])                 # 距離座標のX軸は日付軸にする
    ax.set_ylim([y[0],y[-1]])                 # 距離座標のY軸は緯度軸にする
    p_on_fig = (0.5,0.5)
    p_on_pix = ax.transAxes.transform(p_on_fig)            # figure座標-->pixel座標
    p_on_dis = ax.transData.inverted().transform(p_on_pix) # pixel座標-->距離座標

    print(p_on_fig)
    print(p_on_pix)
    print(p_on_dis)
    print(dt.datetime.utcfromtimestamp(p_on_dis[0]*86400))  # utcfromtimestampで1970.1.1 00Zからの秒数から日付を返す

    dt_aware = dt.datetime(1970, 1, 11, 0, tzinfo=dt.timezone.utc)  # 1970.1.1 00Zからの秒数を返す
    p_on_dis = (dt_aware.timestamp()/86400,35.0)
    p_on_pix = ax.transData.transform(p_on_dis)            # pixel座標 <-- 距離座標
    p_on_fig = ax.transAxes.inverted().transform(p_on_pix) # figure座標 <-- pixel座標
    
    print(p_on_dis)
    print(p_on_pix)
    print(p_on_fig)
    
    return
convert_fig_xy()
以下は出力。
(0.5, 0.5)
[825. 450.]
[ 9.5 40. ]
1970-01-10 12:00:00
[0.52631579 0.375     ]
[859.73684211 360.        ]
(10.0, 35.0)

(6) 棒グラフ

自然災害が多いこの世、風水害と雪害とで死者・行方不明者が同程度であることがわかります、という図の書き方です。内閣府から令和2年度防災白書に付随したエクセルファイルをダウンロードします。
def draw_bar_graph():
    import numpy as np
    import matplotlib.pyplot as plt 
    import pandas as pd # excelファイルはpandasで処理するのがやりやすい
    bousai_csv = "http://www.bousai.go.jp/kaigirep/hakusho/r02/csv/siryo08.csv"
    dead_miss_number = pd.read_csv(bousai_csv,header=1,encoding="shift-jis")
    
    dead_wind = dead_miss_number.iloc[:26,1].values  # .valuesで値をたぶんarrayで返す
    dead_snow = dead_miss_number.iloc[:26,4].values  # .iloc(line,column)で指定できる(便利!)
    time = np.arange(1993,2019) # 1993,2019-1までの1刻みの数列

    fig = plt.figure(figsize=(11/2.54,5/2.54),dpi=300)  # figure領域を 11cm x 5cm に指定
    ax = plt.axes((0.1,0.1,0.8,0.8))                    # axes領域を0.1-0.9, 0.1-0.9に指定
    
    p1 = ax.bar(time, dead_wind, color="k")
    p2 = ax.bar(time, dead_snow, bottom=dead_wind, color="gray") # 棒グラフに乗せる
    ax.legend((p1[0], p2[0]), ("Wind/Water", "Snow"))   # 凡例をつける
    
    ax.set_xlim([time[0]-1,time[-1]+1])  # X軸の範囲指定
    ax.set_ylim([0,400])                 # Y軸の範囲指定
    ax.set_title("Dead and missing")     # タイトル
    
    plt.savefig("Bar_graph.jpg",bbox_inches="tight",pad_inches=0)  # jpegで保存
    fig.show()
    
    return
draw_bar_graph()

完成した図はこのような感じです。