RuntaScience diary

気象系データを扱う学生　旅が好きです

【Numpy&Matplotlib】相関プロットでエラーバーと回帰直線を表示しよう

Python Python-NumPy Python-Matplotlib

こんにちは！ エラーバーはデータの不確実性を示すため、データを示し際には重要です。
今回は、エラーバーを表示する方法をpythonを用いて実装したいと思います。

相関プロット

相関プロットは、高校数学のデータの分析においてよく見たグラフだと思います。
相関係数によって、比較対象の2つの要素の関係がわかります。

データの準備

まずは、以下の2つをインポートします。 ''' import numpy as np import matplotlib.pyplot as plt ''' numpy: 計算ツール
matplotlib: 描写ツール

xとyのデータと、それぞれの標準偏差の
データを用意します。
''' x = np.arange(1,12,1) #1 y = np.array([2,3,6,5,10,10,12,14,18,20,20]) #2 x_err = np.array([0.6, 0.6, 0.6, 0.6, 0.7, 0.9, 0.7, 0.9, 0.8, 0.2, 0.5]) y_err = np.array([0.6, 0.7, 0.5, 0.8, 0.5, 0.3, 0.2, 0.3, 0.6, 0.9, 0.7]) '''

1) np.arange(開始の数, 終りの数, 間隔)
2) np.array([list])　リストをArrayにする。

エラーバー

今回、エラーバーはデータの標準偏差(Standard Deviation)を使います。

#相関係数計算
corr = np.corrcoef(x,y)

相関係数の計算は、np.corrcoef(xのデータ,yのデータ)
でできます。xとyの要素数は揃えないと計算ができません。

corrは配列で出力されます。
⇒array([[1. , 0.9852236],
　　[0.9852236, 1. ]])

したがって、相関係数は、

corr[0, 1]
corr[1, 0]

のどちらかとします。

グラフ作成

グラフを書いてみます。

fig = plt.figure(figsize=(15,5))
plt.rcParams["font.size"] = 18
plt.suptitle("Errobar")

ax1 = plt.subplot(131)
ax2 = plt.subplot(132)
ax3 = plt.subplot(133)

ax1.errorbar(x, y, xerr=x_err, yerr=y_err, fmt = "o"
             ,markersize = 10,color="k", markerfacecolor="w",capsize=8)
ax2.errorbar(x, y, xerr=x_err, fmt = "o"
             ,markersize = 10,color="k", markerfacecolor="w",capsize=8)
ax3.errorbar(x, y, yerr=y_err, fmt = "o"
             ,markersize = 10,color="k", markerfacecolor="w",capsize=8)

axes = [ax1, ax2, ax3]
for ax in axes:
    #範囲の設定
    ax.set_xlim(0, 12)
    ax.set_ylim(0, 22)

    #メモリの設定
    ax.minorticks_on() #補助メモリの描写
    ax.tick_params(axis="both", which="major",direction="in",length=5,width=2,top="on",right="on")
    ax.tick_params(axis="both", which="minor",direction="in",length=2,width=1,top="on",right="on")

    #ラベルの設定
    ax.set_xlabel("X-axis")
    ax.set_ylabel("Y-axis")

    #テキストの貼り付け
    ax.text(0.5, 18, "R={:.3f}".format(corr[0,1]))
    
    ax.label_outer()

plt.subplots_adjust(wspace=0.1)
plt.show()

#保存
fig.savefig("XXX.png",format="png", dpi=330)

f:id:RuntaScience:20200519141253p:plain

回帰直線

1次の回帰なので、y=ax+bが求める回帰直線です。 以下の操作で計算できます。

#回帰直線
p = np.polyfit(x, y, 1)
y_reg = x*p[0]+p[1]

np.polyfit(xのデータ, yのデータ, 次元)
p[0]=傾き
p[1]=切片

グラフ作成

グラフを書いてみます。

fig = plt.figure(figsize=(15,5))
plt.rcParams["font.size"] = 18
plt.suptitle("FittingLine")

ax1 = plt.subplot(131)
ax2 = plt.subplot(132)
ax3 = plt.subplot(133)

ax1.plot(x, y_reg, "--" ,color="r")
ax2.plot(x, y_reg, "-" ,color="r")
ax3.plot(x, y_reg, ".-" ,color="r")

axes = [ax1, ax2, ax3]
for ax in axes:
    #範囲の設定
    ax.set_xlim(0, 12)
    ax.set_ylim(0, 22)

    #メモリの設定
    ax.minorticks_on() #補助メモリの描写
    ax.tick_params(axis="both", which="major",direction="in",length=5,width=2,top="on",right="on")
    ax.tick_params(axis="both", which="minor",direction="in",length=2,width=1,top="on",right="on")

    #ラベルの設定
    ax.set_xlabel("X-axis")
    ax.set_ylabel("Y-axis")

    #テキストの貼り付け
    ax.text(0.5, 20, "y={:.2f}".format(p[0])+"x{0:+.2f}".format(p[1]))
    
    ax.label_outer()

plt.subplots_adjust(wspace=0.1)
plt.show()

#保存
fig.savefig("XXX.png",format="png", dpi=330)

f:id:RuntaScience:20200519141147p:plain

相関グラフまとめ

def main():
    fig = plt.figure(figsize=(10,10))
    plt.rcParams["font.size"] = 18

    ax = plt.subplot(111)

    ax.errorbar(x, y, xerr=x_err, yerr=y_err, fmt = "o"
                 ,markersize = 10,color="k", markerfacecolor="w",capsize=8)
    ax.plot(x, y_reg, color="r")

    #範囲の設定
    ax.set_xlim(0, 12)
    ax.set_ylim(0, 22)

    #メモリの設定
    ax.minorticks_on() #補助メモリの描写
    ax.tick_params(axis="both", which="major",direction="in",length=5,width=2,top="on",right="on")
    ax.tick_params(axis="both", which="minor",direction="in",length=2,width=1,top="on",right="on")
    
    #ラベルの設定
    ax.set_title("Correlation")
    ax.set_xlabel("X-axis")
    ax.set_ylabel("Y-axis")

    #テキストの貼り付け
    ax.text(0.5, 20, "y={:.2f}".format(p[0])+"x{0:+.2f}".format(p[1]))
    ax.text(0.5, 18, "R={:.3f}".format(corr[0,1]))

    plt.show()
    
    #保存
    fig.savefig("XXX.png",format="png", dpi=330)

if __name__ == "__main__":
    main()

f:id:RuntaScience:20200519135551p:plain

参考

matplotlib.org

それでは 🌏

【Matplotlib】Python Matplotlib時系列プロット(気温ー日平均)

Python 気象 Python-Matplotlib Python-NumPy Python-Pandas

こんにちは！ 今日は時系列プロットを紹介します。

時系列プロット

データの取得

縦軸は気温データで、横軸は時間です。
データは気象庁の気温データ(時間平均値)を使用しました。
気象庁|過去の気象データ・ダウンロード

データフレーム(df)の扱い

私のanacondaは日本語を表示する設定をしていないので、
まずはダウンロードしたデータの日本語を英語に変換しました。
(時間→Time, 気温→Temperature ,.... ) それでは必要なモジュールをいれて、データをdfに入れていきます。

import numpy as np
df = pd.read_csv("data/tokyo_temp_2019.csv", header=4)
df.head()
Time Temperature Info num
0 2019/1/1 1:00 1.4 8 1
1 2019/1/1 2:00 2.1 8 1
2 2019/1/1 3:00 1.5 8 1
3 2019/1/1 4:00 1.4 8 1
4 2019/1/1 5:00 1.1 8 1

ファイルタイプがcsvだったのでread_csvを使います。
(csv=Comma Separated Value カンマ区切りの値)という意味
headerは読み飛ばす行を設定します。

それぞれのデータを見てみます。

print(len(df))
print(df.dtypes)

⇒
8760
Time object
Temperature float64
Info int64
num int64
dtype: object

気温がobjectだったので、タイプを日付に変換し、dfの右に貼り付けます。
名前はTime(dt)とします

df["Time(dt)"] = pd.to_datetime(df["Time"])
df
Time Temperature Info num Time(dt)
0 2019/1/1 1:00 1.4 8 1 2019-01-01 01:00:00
1 2019/1/1 2:00 2.1 8 1 2019-01-01 02:00:00
2 2019/1/1 3:00 1.5 8 1 2019-01-01 03:00:00
3 2019/1/1 4:00 1.4 8 1 2019-01-01 04:00:00
4 2019/1/1 5:00 1.1 8 1 2019-01-01 05:00:00

これで準備は完成です！

グラフ描画

まずは、描写に必要なモジュールをインポートします。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import dates as mdates *1
from matplotlib.dates import DateFormatter *1
import datetime *1

1.日付の設定に必要なモジュール

描写していきます！

#データの選択
x = df["Time(dt)"]
y = df["Temperature"]

#グラフ作成
fig = plt.figure(figsize=(20,5))
#fontsizeを一括管理
plt.rcParams["font.size"] = 16

ax = plt.subplot(111)

#データのプロット
ax.plot(x, y, linewidth=1)

#ラベル設定
ax.set_title("Temperature 2019 Tokyo")
ax.set_ylabel("Temperature ($\mathrm{^\circ C}$)")
ax.set_xlabel("Time")

#軸の設定
ax.minorticks_on()  #補助目盛も設定に必要*2
#*3
ax.tick_params(axis="both", which="major",direction="in",length=7,width=2,top="on",right="on")
ax.tick_params(axis="both", which="minor",direction="in",length=4,width=1,top="on",right="on")
#*4-1,2,3
ax.xaxis.set_major_formatter(DateFormatter("%b-%d"))
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=2)) 
ax.xaxis.set_minor_locator(mdates.MonthLocator(interval=1)) 

#範囲設定
ax.set_ylim(y.min()*1.2, y.max()*1.2) #気温の最小値、最大値の1.2倍を範囲とした
ax.set_xlim(datetime.datetime(2019, 1, 1), datetime.datetime(2019, 12, 31)) #*5

plt.show()

fig.savefig("XXX.png",format="png", dpi=330)

f:id:RuntaScience:20200518114013p:plain

2) x軸の後に入れると、時間軸に変な補助線が入るので、必ず前に！
3) axis: "both", "y", "x"を入れる
　which: majorは主目盛、minorは補助目盛
　direction: メモリの方向"in", "out"
　length: メモリの長さ
　width: メモリの太さ
　top, right, left, bottom: どの方向を利用するか
4-1) 主目盛のフォーマット
　%b-%dで、月(英語表示)-日
4-2) 主目盛の間隔
　MonthLocaterで月間隔
　intervalで間隔を設定
4-3) 補助目盛の間隔
　MonthLocaterで月間隔
　intervalで間隔を設定
5) datetime.datetimeで範囲設定

参考

時間軸の設定
Date tick labels — Matplotlib 3.3.3 documentation

matplotlib.dates — Matplotlib 3.3.3 documentation

それでは 🌏

【Jupyter Notebook】Python Jupyter Notebook-単位表示のお話

Python Python-Basic Python-Matplotlib

こんにちは！

Jupyternotebookを用いた解析でグラフを作成する際、ラベルの表示に気を使います。
具体的には、単位と物理量です。
それぞれよく使うものを載せます。

単位の話

単位表示は
単位は立体、ローマン体で
物理量は斜体で

立体

立体は傾いていない表示です。 JupyterNoteでそのままラベルを表示させると、斜体になってしまう場合があります。
それに注意して、例えばμgは
⇒○μg
⇒×μg

斜体

斜体は、イタリック体のことで傾いている表示です
たとえば、F=ma は
⇒ F=ma
と表します。

具体例

気象で使う単位を具体例として載せます (2020/5現在)

斜体は$\it{文字}$
立体は$\mathrm{文字}$
です！

種類 表示 書き方
Temperature a $\mathrm{^\circ C}$
WindSeed b $\mathrm{(m \cdot s^{-1})}$
Concentration c $\mathrm{(\mu g \, m^{-3})}$
Pressure d $\mathrm{(hPa)}$
Flux e $\mathrm{(W \, m^{-2})}$
F=ma 1 $\it{F} = \it{ma}$

f:id:RuntaScience:20200518121811p:plain

参考

LaTeXと、matplotlibのレファレンスを参考にして
数式を書くといいでしょう

LaTeXコマンド集
https://matplotlib.org/gallery/text_labels_and_annotations/mathtext_examples.html?highlight=mathrm%20italic

それでは 🌏

【Cartopy】Python Cartopyを使ったMapping

Python Python-Cartopy Python-Matplotlib

こんにちは！ 私は気象のデータを用いて研究を行っています。

Pythonを用いて、Jupyternotebookでコーディングを行っています 今回はCartopyの使い方の基本をメモします。

Cartopyについて

今まではBasemapを使って解析をしていたのですが、今回Cartopyに切り替えることにしました。

Cartopyの基本的な使い方は、ホームページに掲載されていますので参照ください！

Using cartopy with matplotlib — cartopy 0.19.0rc2.dev8+gd251b2f documentation

グラフの作成

アジア～日本をマッピング

地図は正距円筒図法で書きます。 まず、これらをインポートします

#*1
import matplotlib.pyplot as plt       
#*2
import cartopy.crs as ccrs               
import cartopy.feature as cfeature

Keys
1. グラフを描くツール
2. Cartopyでの描写で用いるツール

def main():
    fig = plt.figure(figsize=(15,5))
    plt.rcParams["font.size"] = 18 #*3
    plt.suptitle("Map") #全体のタイトル設定
    
    #図の用意
    ax1 = fig.add_subplot(1,3,1, projection=ccrs.PlateCarree()) 
    #*4
    ax1.set_extent([90, 150, 0, 60], crs=ccrs.PlateCarree()) 
    #タイトルの設定
    ax1.set_title("Asia")
    
    ax2 = fig.add_subplot(1,3,2, projection=ccrs.PlateCarree())
    ax2.set_extent([128, 148, 30, 50], crs=ccrs.PlateCarree())
    ax2.set_title("Japan")
    
    ax3 = fig.add_subplot(1,3,3, projection=ccrs.PlateCarree())
    ax3.set_extent([139, 141, 34.5, 36.5], crs=ccrs.PlateCarree())
    ax3.set_title("Kanto")
    
    #まとめて設定するものはforループで!
    axes = [ax1, ax2, ax3]
    for ax in axes:
        #海岸線の解像度を10 mにする*5
        ax.coastlines(resolution='10m') 
        #国境線を入れる*6
        ax.add_feature(cfeature.BORDERS, linestyle=':')

    plt.show()
    
    #図の保存
    fig.savefig('XXX.png', format='png', dpi=360)

if __name__ == '__main__':
    main()

Keys
3. 一括でフォントサイズ設定
4. 図の範囲[経度(min), 経度(max), 緯度(min), 緯度(max)]
5. resolutionは、"110", "50","10"のどれか
6. linestyleは、":","--","-"など

f:id:RuntaScience:20200516152747p:plain 以下のコードを参考にしました。 https://scitools.org.uk/cartopy/docs/latest/gallery/features.html#sphx-glr-gallery-features-py

県庁所在地のプロット

試しに県庁所在地をマップにプロットします。 県庁所在地のデータは以下のＨＰを利用いたしました。

【みんなの知識 ちょっと便利帳】都道府県庁所在地 緯度経度データ - 各都市からの方位地図 - 10進数/60進数での座標・世界測地系(WGS84)

まずは使うライブラリをインポートします。

import pandas as pd

次にエクセルファイルからデータをdf(データフレーム)に入れます

df = pd.read_excel("latlng_data.xls", skiprows=4, skipfooter=6, index_col=0)
df.head() #dfのはじめのみの表示

Keys
1. skiprows: はじめの行数をスキップ
2. skipfooter: 終わりの行数をスキップ
3. index_col: インデックスにする行を指定する

データを用意したら、緯度(latitude)と経度(longitude)を切りとります

lat = df["緯度"]
lon = df["経度"]

グラフを描きます 今回はScatter(散布)で点を打っていきます

def main():
    fig = plt.figure(figsize=(5,5))
    plt.rcParams["font.size"] = 18
    
    ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree())
    ax.set_extent([128, 148, 30, 50], crs=ccrs.PlateCarree())
    ax.set_title("Japan")

    ax.coastlines(resolution='10m') #海岸線の解像度を10 mにする
    ax.add_feature(cfeature.BORDERS, linestyle=':')
    ax.stock_img() #地図の色を塗る
    
    #データのプロット
    ax.scatter(lon, lat, color="r", marker="o", s = 5)

    plt.show()
    
    #保存
    fig.savefig('XXX.png', format='png', dpi=360)

if __name__ == '__main__':
    main()

f:id:RuntaScience:20200516154956p:plain

今日はここまで

それでは 🌏

プライバシーポリシー
お問い合わせ