【データ解析】Pythonでデータ解析[基礎]ー数と文字・計算・配列・行列
こんにちは
今日はpythonでデータ解析の基礎です
数と文字
主に用いられる型(type)
整数(int)
a = 4
type(a)
>>>int浮動小数点(float)
b = 4.4
type(b)
>>>float
b2 = 1/3
type(b2)
>>>float真偽(True or False)
c = True
type(c)
>>>True
文字列(str)
d = "NAME"
type(d)
>>>str
日付(datetime)
import pandas as pd
e = pd.to_datetime("2019-01-01 12:30")
print(e)
type(e)
>>>2019-01-01 12:30:00
pandas._libs.tslibs.timestamps.Timestamp
リスト(list)
f = [1, 2, 3, 4]
type(f)
>>>list
辞書
g = {"A1":1, "A2":2, "A3":3, "A5":4}
type(g)
>>>dict
タプル
h = (1, 2, 3, 4)
type(h)
>>>tuple
セット
i = {1, 2, 3, 4}
type(i)
>>>set
配列(array)
import numpy as np
j = np.arange(1,5,1)
type(j)
>>>numpy.ndarray
関数
def my_name(name):
print(name)
my_name("Alex")
type(my_name)
>>>Alex
function
計算
基本的な計算から、行列計算までです
加法・減法・乗法・除法
1 + 4
>>>5
4 - 1
>>>3
11 * 12
>>>132
33 / 3
>>>11.0
11 % 2
>>>1
配列
np.arange(x, y, n)
xからyまでのn間隔の配列の作成
注意するのは、y-1しか数が並ばないこと。
a = np.arange(1, 10, 1)
a
>>>array([1, 2, 3, 4, 5, 6, 7, 8, 9])
a2 = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9])
a2
>>>array([1, 2, 3, 4, 5, 6, 7, 8, 9])
.reshape(行, 列)で配列の次元を変更できます。
#2次元配列作成
b = np.arange(1, 51, 1).reshape(5, 10)
b
>>>array([[ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10],
[11, 12, 13, 14, 15, 16, 17, 18, 19, 20],
[21, 22, 23, 24, 25, 26, 27, 28, 29, 30],
[31, 32, 33, 34, 35, 36, 37, 38, 39, 40],
[41, 42, 43, 44, 45, 46, 47, 48, 49, 50]])
.reshape(-1,)で1次元配列に戻せます
b2 = b.reshape(-1,)
b2
>>>array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34,
35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50])
最後に配列の情報を確認しておきましょう。
print(b.shape) #配列の形
print(b.ndim) #配列の次元
print(b.size) #配列のサイズ
>>>(5, 10)
>>>2
>>>50
行列計算
基礎
A = np.array([[2, 2], [2, 4]])
A
>>>array([[2, 2],
[2, 4]])
A + 10
>>>array([[12, 12],
[12, 14]])
A + A
>>>array([[4, 4],
[4, 8]])
A * 3
>>>array([[ 6, 6],
[ 6, 12]])
※行列の掛け算ではなく、要素の掛け算である
A * A
>>>array([[ 4, 4],
[ 4, 16]])
行列の積(3種類の計算方法)
A @ A
np.dot(A, A)
A.dot(A)
>>>array([[ 8, 12],
[12, 20]])
転置行列・逆行列・単位行列
np.transpose(A)
>>>array([[2, 2],
[2, 4]])
inv_A = np.linalg.inv(A)
inv_A
>>>array([[ 1. , -0.5],
[-0.5, 0.5]]) A.dot(inv_A)
# np.dot(A, inv_A)
>>>array([[1., 0.],
[0., 1.]])
B = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
B
>>>array([[1, 2, 3],
[4, 5, 6],
[7, 8, 9]])
#逆行列
inv_B = np.linalg.inv(B)
inv_B
>>>array([[ 3.15251974e+15, -6.30503948e+15, 3.15251974e+15],
[-6.30503948e+15, 1.26100790e+16, -6.30503948e+15],
[ 3.15251974e+15, -6.30503948e+15, 3.15251974e+15]])#元の行列×逆行列
B.dot(inv_B)
>>>array([[ 0. , 1. , -0.5],
[ 0. , 2. , -1. ],
[ 0. , 3. , 2.5]])⇑単位行列が求まっていない
それは、元の行列の行列式≠0を満たしていないからである
int(np.linalg.det(B))
>>>0行列式
np.linalg.det(A)
>>>4.0
固有値と固有ベクトル
A = np.array([[2, 2], [2, 4]])
A
>>>array([[2, 2],
[2, 4]])
w, v = np.linalg.eig(A)
print(w) #固有値
print(v) #固有ベクトル
>>>[0.76393202 5.23606798]
[[-0.85065081 -0.52573111]
[ 0.52573111 -0.85065081]]
それでは🌏
【Pandas】 データをファイルに格納(保存; WRITE)
こんにちは。今日は作成したデータを保存する方法です。
データフレーム
まずはモジュールをインポートします。
import pandas as pd
使用するデータは、以下のデータフレームです。
print(df)>>>
| Year | Month | Day | Data | |
|---|---|---|---|---|
| 0 | 2019 | 1 | 1 | 12 |
| 1 | 2019 | 2 | 1 | 14 |
| 2 | 2019 | 3 | 1 | 13 |
| 3 | 2019 | 4 | 1 | 12 |
| 4 | 2019 | 5 | 1 | 13 |
| 5 | 2019 | 6 | 1 | 11 |
| 6 | 2019 | 7 | 1 | 14 |
| 7 | 2019 | 8 | 1 | 12 |
| 8 | 2019 | 9 | 1 | 14 |
| 9 | 2019 | 10 | 1 | 14 |
| 10 | 2019 | 11 | 1 | 11 |
| 11 | 2019 | 12 | 1 | 15 |
テキストファイル
indexまで保存されるので、読みこみの時は、indexを指定しましょう
保存時にindexを除きたい場合は、
df.to_csv("XXX.txt", index=False)
としましょう。
カンマ区切り
#データの保存
df.to_csv("XXX.txt")
#確認
pd.read_csv("XXX.txt", index_col=0)
タブ区切り
#データの保存
df.to_csv("XXX.txt",sep="\t")
#確認
pd.read_table("XXX.txt",index_col=0)
エクセルファイル
#データの保存
df.to_excel("XXX.xlsx",sheet_name="X")
#確認
pd.read_excel("XXX.xlsx",index_col=0,sheet_name="X")
複数データを同時に
複数のデータフレームを1つのファイルにシートを分けて保存するには、ExcelWriterを使います。
例として、データフレームをdf1, df2, df3とします
path = "XXX.xlsx"
with pd.ExcelWriter(path) as writer:
df1.to_excel(writer, sheet_name = "sheet1")
df2.to_excel(writer, sheet_name = "sheet2")
df3.to_excel(writer, sheet_name = "sheet3")
バイナリデータ(配列;Array)
まずはモジュールをインポートします。
import numpy as np今回は1次元の配列データと2次元の配列データを保存します。
lon = np.arange(-180,180,0.1)
lat = np.arange(-90, 90, 0.1)
two_dim = np.arange(0,150,1).reshape(10,15)
print(lon, lat, two_dim)
>>>[-180. -179.9 -179.8 ... 179.7 179.8 179.9]
>>>[-90. -89.9 -89.8 ... 89.7 89.8 89.9]
>>>[[ 0 1 ... 12 13 14][ 15 16 ... 27 28 29]...
...[120 121 ... 132 133 134] [135 136 ... 147 148 149]]
保存します。
np.save("XXX1.npy", lon)
np.save("XXX2.npy", lat)
np.save("XXX3.npy", two_dim)読み込んでデータがしっかり入っていることを確認しましょう
np.load("XXX1.npy")
>>>[-180. -179.9 -179.8 ... 179.7 179.8 179.9]
読み込みに関しては前の記事で紹介しています。
それでは🌏
【Pandas】ファイルをデータフレームに格納(読み込み; READ)
こんにちは。今日はファイルの読み込みについてです。
テキストファイル
まずはモジュールをインポートします。
import numpy as np
カンマ区切り
pd.read_csv("XXX.txt")
カンマで区切ってあるテキストファイルをインポートします。
データ⇓
Year, Month, Day, Data 2019, 1, 1, 12 2019, 2, 1, 14 2019, 3, 1, 13 2019, 4, 1, 12 2019, 5, 1, 10 2019, 6, 1, 10 2019, 7, 1, 15 2019, 8, 1, 12 2019, 9, 1, 10 2019, 10, 1, 12 2019, 11, 1, 10 2019, 12, 1, 12
df = pd.read_csv("XXX.txt")
df
| Year | Month | Day | Data | |
|---|---|---|---|---|
| 0 | 2019 | 1 | 1 | 12 |
| 1 | 2019 | 2 | 1 | 14 |
| 2 | 2019 | 3 | 1 | 13 |
| 3 | 2019 | 4 | 1 | 12 |
| 4 | 2019 | 5 | 1 | 10 |
| 5 | 2019 | 6 | 1 | 10 |
| 6 | 2019 | 7 | 1 | 15 |
| 7 | 2019 | 8 | 1 | 12 |
| 8 | 2019 | 9 | 1 | 10 |
| 9 | 2019 | 10 | 1 | 12 |
| 10 | 2019 | 11 | 1 | 10 |
| 11 | 2019 | 12 | 1 | 12 |
タブ区切り
pd.read_table("XXX.txt")
タブ区切りのテキストファイルも同様にデータフレームに格納できます
Year Month Day Data 2019 1 1 12 2019 2 1 14 2019 3 1 13 2019 4 1 12 2019 5 1 13 2019 6 1 11 2019 7 1 14 2019 8 1 12 2019 9 1 14 2019 10 1 14 2019 11 1 11 2019 12 1 15
1)
df = pd.read_table("XXX.txt")
df
2)
df = pd.read_csv("XXX.txt",sep="\t")
df
| Year | Month | Day | Data | |
|---|---|---|---|---|
| 0 | 2019 | 1 | 1 | 12 |
| 1 | 2019 | 2 | 1 | 14 |
| 2 | 2019 | 3 | 1 | 13 |
| 3 | 2019 | 4 | 1 | 12 |
| 4 | 2019 | 5 | 1 | 10 |
| 5 | 2019 | 6 | 1 | 10 |
| 6 | 2019 | 7 | 1 | 15 |
| 7 | 2019 | 8 | 1 | 12 |
| 8 | 2019 | 9 | 1 | 10 |
| 9 | 2019 | 10 | 1 | 12 |
| 10 | 2019 | 11 | 1 | 10 |
| 11 | 2019 | 12 | 1 | 12 |
エクセルファイル
エクセルファイルも同様のデータを使います。
基本
| Year | Month | Day | Data |
| 2019 | 1 | 1 | 12 |
| 2019 | 2 | 1 | 14 |
| 2019 | 3 | 1 | 13 |
| 2019 | 4 | 1 | 12 |
| 2019 | 5 | 1 | 13 |
| 2019 | 6 | 1 | 11 |
| 2019 | 7 | 1 | 14 |
| 2019 | 8 | 1 | 12 |
| 2019 | 9 | 1 | 14 |
| 2019 | 10 | 1 | 14 |
| 2019 | 11 | 1 | 11 |
| 2019 | 12 | 1 | 15 |
df = pd.read_excel("XXX.xlsx")
dfシート名・インデックスを指定したい場合
sheet ="b"でインデックスをidxで指定する場合は、
df = pd.read_excel("XXX.xlsx", sheet_name="b",index_col=0)
df
| Year | Month | Day | Data | |
|---|---|---|---|---|
| idx | ||||
| 1 | 2019 | 1 | 1 | 12 |
| 2 | 2019 | 2 | 1 | 14 |
| 3 | 2019 | 3 | 1 | 13 |
| 4 | 2019 | 4 | 1 | 12 |
| 5 | 2019 | 5 | 1 | 13 |
| 6 | 2019 | 6 | 1 | 11 |
| 7 | 2019 | 7 | 1 | 14 |
| 8 | 2019 | 8 | 1 | 12 |
| 9 | 2019 | 9 | 1 | 14 |
| 10 | 2019 | 10 | 1 | 14 |
| 11 | 2019 | 11 | 1 | 11 |
| 12 | 2019 | 12 | 1 | 15 |
データの一部のみ使いたい
| Year | Month | Day | Data_x1 | Data_x2 | Data_y1 | Data_y2 |
| 2019 | 1 | 1 | 12 | 12 | 12 | 12 |
| 2019 | 2 | 1 | 14 | 14 | 14 | 14 |
| 2019 | 3 | 1 | 13 | 13 | 13 | 13 |
| 2019 | 4 | 1 | 12 | 12 | 12 | 12 |
| 2019 | 5 | 1 | 13 | 13 | 13 | 13 |
| 2019 | 6 | 1 | 11 | 11 | 11 | 11 |
| 2019 | 7 | 1 | 14 | 14 | 14 | 14 |
| 2019 | 8 | 1 | 12 | 12 | 12 | 12 |
| 2019 | 9 | 1 | 14 | 14 | 14 | 14 |
| 2019 | 10 | 1 | 14 | 14 | 14 | 14 |
| 2019 | 11 | 1 | 11 | 11 | 11 | 11 |
| 2019 | 12 | 1 | 15 | 15 | 15 | 15 |
1)データを含む
include=["Data_x1"]
df = pd.read_excel("XXX.xlsx", sheet_name="c",usecols=include)
df
| Data_x1 | |
|---|---|
| 0 | 12 |
| 1 | 14 |
| 2 | 13 |
| 3 | 12 |
| 4 | 13 |
| 5 | 11 |
| 6 | 14 |
| 7 | 12 |
| 8 | 14 |
| 9 | 14 |
| 10 | 11 |
| 11 | 15 |
2)データを抜く
exclude=["Data_x2", "Data_y2"] #抜く列
df = pd.read_excel("XXX.xlsx", sheet_name = "c",usecols = lambda x: x not in exclude)
df| Year | Month | Day | Data_x1 | Data_y1 | |
|---|---|---|---|---|---|
| 0 | 2019 | 1 | 1 | 12 | 12 |
| 1 | 2019 | 2 | 1 | 14 | 14 |
| 2 | 2019 | 3 | 1 | 13 | 13 |
| 3 | 2019 | 4 | 1 | 12 | 12 |
| 4 | 2019 | 5 | 1 | 13 | 13 |
| 5 | 2019 | 6 | 1 | 11 | 11 |
| 6 | 2019 | 7 | 1 | 14 | 14 |
| 7 | 2019 | 8 | 1 | 12 | 12 |
| 8 | 2019 | 9 | 1 | 14 | 14 |
| 9 | 2019 | 10 | 1 | 14 | 14 |
| 10 | 2019 | 11 | 1 | 11 | 11 |
| 11 | 2019 | 12 | 1 | 15 | 15 |
3)header飛ばし
| data:2019/1/1 | |||
| Site:Japan | |||
| Year | Month | Day | Data |
| 2019 | 1 | 1 | 12 |
| 2019 | 2 | 1 | 14 |
| 2019 | 3 | 1 | 13 |
| 2019 | 4 | 1 | 12 |
| 2019 | 5 | 1 | 13 |
| 2019 | 6 | 1 | 11 |
| 2019 | 7 | 1 | 14 |
| 2019 | 8 | 1 | 12 |
| 2019 | 9 | 1 | 14 |
| 2019 | 10 | 1 | 14 |
| 2019 | 11 | 1 | 11 |
| 2019 | 12 | 1 | 15 |
df = pd.read_excel("XXX.xlsx", sheet_name = "a",header = 3)
df
| Year | Month | Day | Data | |
|---|---|---|---|---|
| 0 | 2019 | 1 | 1 | 12 |
| 1 | 2019 | 2 | 1 | 14 |
| 2 | 2019 | 3 | 1 | 13 |
| 3 | 2019 | 4 | 1 | 12 |
| 4 | 2019 | 5 | 1 | 13 |
| 5 | 2019 | 6 | 1 | 11 |
| 6 | 2019 | 7 | 1 | 14 |
| 7 | 2019 | 8 | 1 | 12 |
| 8 | 2019 | 9 | 1 | 14 |
| 9 | 2019 | 10 | 1 | 14 |
| 10 | 2019 | 11 | 1 | 11 |
| 11 | 2019 | 12 | 1 | 15 |
netCDFファイル
まずはモジュールをインポート
from netCDF4 import Dataset
nc = Dataset("XXX.nc4", mode="r")
nc
<データの情報が出てくる>・・・・データの情報を基に、データをarrayにする。
data = nc.variables["dataname"][:].data #1次元データの場合
data
>>>array([100,120,122,111,92,222,…, 128, 100])
それでは🌏
Lightroomで動物の背景を黒くする方法
こんにちは！
今回は、写真の編集方法に関してです。
撮影
まずは編集前です。
シャッタースピード：1/800、f：5.6、ISO：1000
背景を黒くすることを考えて、主役の動物が映えるように撮ります。
動物がぶれないようにシャッタースピードを上げます。
編集
ライト
コントラストを高めて、ハイライト、シャドウ、黒レベルをマイナスに落とすと
主役が映えるようになります。
ハイライトを下げると、羽などの陰がシャープになります。
シャドウを下げると主に背景が黒みがかります。
黒レベルを下げると、全体的に黒くなり、白が映えます。
カラー
カラーは生き物によって変えましょう。
過度な彩度変化は施さなくても大丈夫だと思います。
より美しく見えるように調節します。
作品例
対象が白で無い場合も編集可能です。
また、Rawで撮らなくても上記のように編集ができます。
大切なこと
一番大切なことは、編集後をイメージすることだと思います。
背景を黒くするために撮った写真は、編集すると美しいですが、編集前は色が薄い写真となります。したがって、撮る段階が一番大切だと感じます。
それでは🌏
【Matplotlib】Python Format文～必要なものだけ厳選(理系)
こんにちは！
今回は、pythonでよく用いるformat文についてです
使い方
基本的な使い方は簡単です。
"{}"の中身を.format()で書き入れましょう
"{}".format(2019)
>>>'2019'
1つだけではなく、複数の文字や数を入れられます。
0番から始まる連番で指定すると、その数が入ります。
"I'm {0} {1}.".format("Noah", "Saito")
>>>"I'm Noah Saito."
具体例
次によく使うものを具体例で見ていきましょう。
小数点
小数点は「f」で表現します。fの前に何桁か指定できます。デフォルトは6桁です。
a = 0.2183
"{0:f} || {0:.1f}".format(a)
>>>'0.218300 || 0.2'
同様に指数表記のＥもできます。こちらもデフォルトは6桁です。eが大文字か小文字かも指定できます。
a = 0.2183
"{0:e} || {0:.1e} || {0:.1E}".format(a)
>>>'2.183000e-01 || 2.2e-01 || 2.2E-01'
日にち
それぞれの日付のデータを使って以下のように表すことができます。
year, month, day = 2019, 1, 1
"{0}-{1}-{2}".format(year, month, day)
>>>'2019-1-1'
日付を'2019-01-01'のような桁数指定をしたい場合は、以下のように
「位置:0+桁数d」としましょう。(dは10整数である)
year, month, day = 2019, 1, 1
"{0:04d}-{1:02d}-{2:02d}".format(year, month, day)
>>>'2019-01-01'
さらに、モジュールのdatetimeを用いて表現できます。
import datetime
time = datetime.datetime(2019, 1, 1, 1, 1, 30)
"{:%Y-%m-%d %H:%M:%S}".format(time)
'2019-01-01 01:01:30'
ちなみに、月をJan, Feb...のようにしたい場合は、「%b」と指定します。
import datetime
for i in range(12):
time = datetime.datetime(2019, i+1, 1)
print("{:%b}".format(time))
>>>Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
n進数
2進数(binary number)、8進数(octal number)、10進数(decimal number)、16進数(hexadecimal number)があります。
binary, octal, decimal, hexadecimal = 30, 30, 30, 30
"{0:b}-{1:o}-{2:d}-{3:x}(or{3:X})".format(binary, octal, decimal, hexadecimal)
>>>'11110-36-30-1e(or1E)'
％
最後にパーセント表示です。100倍されることに注意してください。
そして、桁数も変更できます。デフォルトは6桁です。
a = 0.2183
"{0:%} || {0:.0%} || {0:.2%}".format(a)
>>>'21.830000% || 22% || 21.83%'
さらに詳しく
今回の記事はこちらを参考にして、理系が用いるものを厳選して紹介しました。
それでは
【Numpy&Matplotlib】 Python 軸を日付フォーマットに変更
こんにちは今日はx軸を日にちのフォーマットに変更する方法です。
DataFlame
まずは使用するデータフレームをファイルから読み込んでください。
今回私は適当に作成したデータを使います。
print(df)| Year | Month | Day | Hour | Minute | Data | |
|---|---|---|---|---|---|---|
| 0 | 2019 | 1 | 1 | 0 | 1 | 1.683375 |
| 1 | 2019 | 1 | 1 | 0 | 6 | 1.732797 |
| 2 | 2019 | 1 | 1 | 0 | 11 | 0.985744 |
| 3 | 2019 | 1 | 1 | 0 | 16 | 1.522657 |
| 4 | 2019 | 1 | 1 | 0 | 21 | 0.726131 |
| ... | ... | ... | ... | ... | ... | ... |
| 105115 | 2019 | 12 | 31 | 23 | 36 | 1.442926 |
| 105116 | 2019 | 12 | 31 | 23 | 41 | 1.532434 |
| 105117 | 2019 | 12 | 31 | 23 | 46 | 1.513058 |
| 105118 | 2019 | 12 | 31 | 23 | 51 | 0.594683 |
| 105119 | 2019 | 12 | 31 | 23 | 56 | 0.767968 |
105120 rows × 6 columns
1日当たりのDataの平均値を時系列プロットしてみます
#日平均
import calendar
xdata = []
year=2019
for j in range(12):
month=j+1
x=calendar.monthrange(year, month)[1] #*1
df2 = df.loc[df["Month"]==month]
for i in range(x):
day = i + 1
xdata0 = df2.loc[df2["Day"]==day]["Data"].mean()
xdata.append(xdata0)Keys
1) year年month月の日数を取ることができる。
例えばyear=2019 month=1 の時、x=31となる。
それではグラフを書いてみます。横軸はDN(Day Number)です。
import matplotlib.pyplot as plt
#グラフ作成
fig = plt.figure(figsize=(20,5))
plt.rcParams["font.size"] = 16
ax = plt.subplot(111)
ax.plot(np.arange(1,366), xdata)
ax.set_xlabel("DN")
ax.set_ylabel("Data")
ax.set_xlim(0,365)
plt.show()
fig.savefig("XXX.png",format="png", dpi=330)
DN⇒Datetime
次に横軸を日付にします。
time = []
year = 2019
for i in range(12):
month = i + 1
x=calendar.monthrange(year, month)[1]
for j in range(x):
day = j + 1
t = "{0:04d}-{1:02d}-{2:02d}".format(year, month, day) #1*
time.append(pd.to_datetime(t)) #*2
Keys
1) format文⇓
2)pd.to_datetime("2019-01-01 01:30:00:00")⇒Timestamp("2019-01-01 01:30:00:00")
オブジェクト⇒datetimeに変換されます
グラフ
#データの選択
x = time
y = xdata
#グラフ作成
fig = plt.figure(figsize=(20,5))
plt.rcParams["font.size"] = 16
ax = plt.subplot(111)
ax.plot(x, y)
#ラベル設定
ax.set_title(" ")
ax.set_ylabel("data")
ax.set_xlabel("Time")
#軸の設定 *1
ax.minorticks_on()
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")
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_xlim(datetime.datetime(2019, 1, 1), datetime.datetime(2019, 12, 31))
plt.show()
fig.savefig("XXX.png",format="png", dpi=330)
Keys
1)軸の設定
それでは 🌏
【Science】 Matplotlibで振り子ー二重振り子 カオスを知る
こんにちは！
今日は二重振り子でカオスを見てみたいと思います！
単振り子に基礎的なコードを載せています。
二重振り子
ドキュメント
モジュール
まずは必要なモジュールです。
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate
import matplotlib.animation as animation
軌跡
まずは二重振り子の軌跡を描いてみます。
%matplotlib nbagg
#物理量設定
G = 9.8 #重力加速度(m/s^2)
L1 = 1.0 #振り子(上)の長さ(m)
L2 = 1.0 #振り子(下)の長さ(m)
M1 = 1.0 #振り子(上)の重さ(kg)
M2 = 1.0 #振り子(下)の重さ(kg)
def derivs(state, t):
dydx = np.zeros_like(state)
dydx[0] = state[1]
delta = state[2] - state[0]
den1 = (M1+M2) * L1 - M2 * L1 * np.cos(delta) * np.cos(delta)
dydx[1] = ((M2 * L1 * state[1] * state[1] * np.sin(delta) * np.cos(delta)
+ M2 * G * np.sin(state[2]) * np.cos(delta)
+ M2 * L2 * state[3] * state[3] * np.sin(delta)
- (M1+M2) * G * np.sin(state[0]))
/ den1)
dydx[2] = state[3]
den2 = (L2/L1) * den1
dydx[3] = ((- M2 * L2 * state[3] * state[3] * np.sin(delta) * np.cos(delta)
+ (M1+M2) * G * np.sin(state[0]) * np.cos(delta)
- (M1+M2) * L1 * state[1] * state[1] * np.sin(delta)
- (M1+M2) * G * np.sin(state[2]))
/ den2)
return dydx
#初期値
dt = 0.05 #時間間隔
t = np.arange(0, 20, dt) #時間間隔分
th1 = 120.0 #振り子(上)初期角度
w1 = 0.0 #振り子(上)角速度
th2 = -10.0 #振り子(下)初期角度
w2 = 0.0 #振り子(下)角速度
state1 = np.radians([th1, w1, th2, w2]) # 初期状態
#微分方程式を解く
y01 = integrate.odeint(derivs, state1, t)
x1 = L1 * np.sin(y01[:, 0])
y1 = -L1 * np.cos(y01[:, 0])
x2 = L2 * np.sin(y01[:, 2]) + x1
y2 = -L2 * np.cos(y01[:, 2]) + y1
#--------------------グラフ-----------------------
fig = plt.figure()
ax = fig.add_subplot(111, autoscale_on=False, xlim=(-2, 2), ylim=(-2, 2))
ax.set_aspect("equal")
ax.grid()
line, = ax.plot([], [], "o-", lw=2)
locus, = ax.plot([], [], "r-", linewidth=2) #軌跡用
time_template = "time = %.1fs"
time_text = ax.text(0.05, 0.9, "", transform=ax.transAxes)
xlocus, ylocus = [], []
def init():
line.set_data([], [])
locus.set_data([], [])
time_text.set_text("")
return line,locus, time_text
def animate(i):
thisx1 = [0, x1[i], x2[i]]
thisy1 = [0, y1[i], y2[i]]
xlocus.append(x2[i])
ylocus.append(y2[i])
line.set_data(thisx1, thisy1)
locus.set_data(xlocus, ylocus)
time_text.set_text(time_template % (i*dt))
return line,locus, time_text
ani = animation.FuncAnimation(fig, animate, range(1, len(y01)),
interval=dt*1000, blit=True, init_func=init)
plt.show()
# ani.save("XXX.gif", writer="pillow", fps=20)
2個
次にいきなり複数の振り子にいくのはなんなので、2つの振り子を描写してみます。
%matplotlib nbagg
#物理量設定
G = 9.8 #重力加速度(m/s^2)
L1 = 1.0 #振り子(上)の長さ(m)
L2 = 1.0 #振り子(下)の長さ(m)
M1 = 1.0 #振り子(上)の重さ(kg)
M2 = 1.0 #振り子(下)の重さ(kg)
def derivs(state, t):
dydx = np.zeros_like(state)
dydx[0] = state[1]
delta = state[2] - state[0]
den1 = (M1+M2) * L1 - M2 * L1 * np.cos(delta) * np.cos(delta)
dydx[1] = ((M2 * L1 * state[1] * state[1] * np.sin(delta) * np.cos(delta)
+ M2 * G * np.sin(state[2]) * np.cos(delta)
+ M2 * L2 * state[3] * state[3] * np.sin(delta)
- (M1+M2) * G * np.sin(state[0]))
/ den1)
dydx[2] = state[3]
den2 = (L2/L1) * den1
dydx[3] = ((- M2 * L2 * state[3] * state[3] * np.sin(delta) * np.cos(delta)
+ (M1+M2) * G * np.sin(state[0]) * np.cos(delta)
- (M1+M2) * L1 * state[1] * state[1] * np.sin(delta)
- (M1+M2) * G * np.sin(state[2]))
/ den2)
return dydx
dt = 0.05
t = np.arange(0, 20, dt)
th1 = 120.0
th12 = 120.1
w1 = 0.0
th2 = -10.0
w2 = 0.0
#初期状態
state1 = np.radians([th1, w1, th2, w2])
state2 = np.radians([th12, w1, th2, w2])
#微分方程式を解く
y01 = integrate.odeint(derivs, state1, t)
y02 = integrate.odeint(derivs, state2, t)
x1 = L1 * np.sin(y01[:, 0])
y1 = -L1 * np.cos(y01[:, 0])
x2 = L2 * np.sin(y01[:, 2]) + x1
y2 = -L2 * np.cos(y01[:, 2]) + y1
x3 = L1 * np.sin(y02[:, 0])
y3 = -L1 * np.cos(y02[:, 0])
x4 = L2 * np.sin(y02[:, 2]) + x3
y4 = -L2 * np.cos(y02[:, 2]) + y3
fig = plt.figure()
ax = fig.add_subplot(111, autoscale_on=False, xlim=(-2, 2), ylim=(-2, 2))
ax.set_aspect("equal")
ax.grid()
line, = ax.plot([], [], "o-", lw=2)
line2, = ax.plot([], [], "o-", lw=2)
time_template = "time = %.1fs"
time_text = ax.text(0.05, 0.9, "", transform=ax.transAxes)
def init():
line.set_data([], [])
line2.set_data([], [])
time_text.set_text("")
return line,line2, time_text
def animate(i):
thisx1 = [0, x1[i], x2[i]]
thisy1 = [0, y1[i], y2[i]]
thisx2 = [0, x3[i], x4[i]]
thisy2 = [0, y3[i], y4[i]]
line.set_data(thisx1, thisy1)
line2.set_data(thisx2, thisy2)
time_text.set_text(time_template % (i*dt))
return line,line2, time_text
ani = animation.FuncAnimation(fig, animate, range(1, len(y01)),
interval=dt*1000, blit=True, init_func=init)
plt.show()
#ani.save("XXX.gif", writer="pillow", fps=20)
カオス
最後にカオスを描写してみます。
先ほどよりも振り子数が一気に増えます。
%matplotlib nbagg
#物理量設定
G = 9.8 #重力加速度(m/s^2)
L1 = 1.0 #振り子(上)の長さ(m)
L2 = 1.0 #振り子(下)の長さ(m)
M1 = 1.0 #振り子(上)の重さ(kg)
M2 = 1.0 #振り子(下)の重さ(kg)
def derivs(state, t):
dydx = np.zeros_like(state)
dydx[0] = state[1]
delta = state[2] - state[0]
den1 = (M1+M2) * L1 - M2 * L1 * np.cos(delta) * np.cos(delta)
dydx[1] = ((M2 * L1 * state[1] * state[1] * np.sin(delta) * np.cos(delta)
+ M2 * G * np.sin(state[2]) * np.cos(delta)
+ M2 * L2 * state[3] * state[3] * np.sin(delta)
- (M1+M2) * G * np.sin(state[0]))
/ den1)
dydx[2] = state[3]
den2 = (L2/L1) * den1
dydx[3] = ((- M2 * L2 * state[3] * state[3] * np.sin(delta) * np.cos(delta)
+ (M1+M2) * G * np.sin(state[0]) * np.cos(delta)
- (M1+M2) * L1 * state[1] * state[1] * np.sin(delta)
- (M1+M2) * G * np.sin(state[2]))
/ den2)
return dydx
dt = 0.05
t = np.arange(0, 20, dt)
th1s = np.arange(120.0,120.01,0.0001)
w1 = 0.0
th2 = -10.0
w2 = 0.0
#初期状態
states=[]
for th1 in th1s:
state = np.radians([th1, w1, th2, w2])
states.append(state)
#微分方程式を解く
ys=[]
for state in states:
y = integrate.odeint(derivs, state, t)
ys.append(y)
x1s,x2s,y1s,y2s=[],[],[],[]
for y in ys:
x1 = L1 * np.sin(y[:, 0])
y1 = -L1 * np.cos(y[:, 0])
x2 = L2 * np.sin(y[:, 2]) + x1
y2 = -L2 * np.cos(y[:, 2]) + y1
x1s.append(x1)
y1s.append(y1)
x2s.append(x2)
y2s.append(y2)
fig = plt.figure()
ax = fig.add_subplot(111, autoscale_on=False, xlim=(-2, 2), ylim=(-2, 2))
ax.set_aspect("equal")
ax.grid()
lines=[]
for i in range(len(th1s)):
line, = ax.plot([], [], "o-", lw=2)
lines.append(line,)
time_template = "time = %.1fs"
time_text = ax.text(0.05, 0.9, "", transform=ax.transAxes)
def init():
for line in lines:
line.set_data([], [])
time_text.set_text('')
return lines, time_text
def animate(i):
for j in range(len(lines)):
line = lines[j]
thisx1 = [0, x1s[j][i], x2s[j][i]]
thisy1 = [0, y1s[j][i], y2s[j][i]]
line.set_data(thisx1, thisy1)
time_text.set_text(time_template % (i*dt))
return lines, time_text
ani = animation.FuncAnimation(fig, animate, range(1, len(y)),
interval=dt*1000, blit=True, init_func=init)
plt.show()
# ani.save("XXX.gif", writer="pillow", fps=20)
それでは🌏
【Science】 Matplotlibで振り子ー単振り子
こんにちは
今日は関数を使って、単振り子をアニメーションを作成したいと思います！
アニメーションの基礎は前の記事を参照してください！
微分方程式
単振り子を描写するのに関数が必要なので書き方を確認しましょう
モジュール
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import numpy as np
解く
それでは、solve_ivpで指数関数的減衰を解いてみます。
指数関数的減衰は半減期の計算で使われますね
https://www.nli-research.co.jp/report/detail/id=58700?site=nli
# 指数関数的減衰
def exponential_decay(t, y):
return -0.5 * y
t_span = [0, 10] #tの範囲
y0 = [2, 4, 8] #y初期値
t_eval = np.arange(t_span[0], t_span[-1], 0.1) #tの間隔
sol = solve_ivp(exponential_decay, t_span, y0, t_eval=t_eval) #方程式を解く
解(sol)の中身を確認します。
print(sol.t)
⇒[0. 0.1 0.2 0.3…9.7 9.8 9.9]
print(sol.y)
⇒[[2. 1.90245885 1.809668…0.01492886 0.01420166]
[4. 3.8049177 3.619336…0.02985772 0.02840332]
[8. 7.6098354 7.238672…0.05971543 0.05680663]]tに対するyが作成されていることがわかります。
描写
それでは描写してみます
fig, ax = plt.subplots()
ax.plot(sol.t, sol.y[0])
ax.plot(sol.t, sol.y[1])
ax.plot(sol.t, sol.y[2])
ax.set_xlabel("T")
ax.set_ylabel("Y-axis")
ax.set_xlim(t_span[0],t_span[-1])
ax.set_ylim(0, 9)
plt.show()
#保存
#fig.savefig("XXX.png", format="png", dpi=330)
単振り子
それでは単振り子のアニメーション作成してみます。
以下を参考に作成しました。
モジュール
まずはモジュールをインポートします
from scipy.integrate import solve_ivp
from matplotlib.animation import FuncAnimation
import matplotlib.pyplot as plt
import numpy as np
単振り子
%matplotlib nbagg
G = 9.8 # 重力加速度(m/s^2)
L = 1.0 # 振り子の長さ(m)
# 運動方程式を解く
def eom(t, state):
dydt = np.zeros_like(state)
dydt[0] = state[1]
dydt[1] = -(G/L)*np.sin(state[0])
return dydt
#設定
t_span = [0,20] #時間 (s)
dt = 0.05# 時間の間隔 (s)
t = np.arange(t_span[0], t_span[-1], dt)# 時間の間隔のarray作成
#初期状態設定
th1 = 45.0# 角度 [deg]
w1 = 0.0 # 初期角速度 [deg/s]
state = np.radians([th1, w1])# 初期状態
#微分方程式を解く
sol = solve_ivp(eom, t_span, state, t_eval=t)
theta = sol.y[0] #y[0]＝θ、y[1]=ω
x = L * np.sin(theta)
y = -L * np.cos(theta)
#ーーーーーーグラフーーーーーーー
fig, ax = plt.subplots()
line, = ax.plot([], [], "o-", linewidth=2)
time_template = "time = %.1fs"
time_text = ax.text(0.05, 0.9, " ", transform=ax.transAxes)
def init():
line.set_data([], [])
time_text.set_text("")
return line, time_text
def animate(i):
thisx = [0, x[i]]
thisy = [0, y[i]]
line.set_data(thisx, thisy)
time_text.set_text(time_template % (i*dt))
return line, time_text
ani = FuncAnimation(fig, animate, frames=np.arange(0, len(t))
, interval=25, blit=True,init_func=init)
ax.set_xlim(-L-0.5,L+0.5)
ax.set_ylim(-L-0.5,L+0.5)
ax.set_aspect("equal")
ax.grid()
plt.show()
#ani.save("XXX.gif", writer="pillow", fps=20)
時間変化
次に単振り子の角度、x, yの時間変化をプロットしてみます
fig = plt.figure(figsize=(10,5))
ax1 = plt.subplot(111)
ax2 = ax1.twinx()
ax1.plot(sol.t, sol.y[0], label="$\\theta$")
ax2.plot(sol.t, x , color="r", linestyle="--",label="x")
ax2.plot(sol.t, y , color="g", linestyle="--",label="y")
ax1.set_yticks(np.array([-np.pi, -(1/3)*np.pi, -(1/4)*np.pi, -(1/6)*np.pi, 0,(1/6)*np.pi,(1/4)*np.pi,(1/3)*np.pi,np.pi]))
ax1.set_yticklabels(["-$\pi$","-1/3$\pi$","-1/4$\pi$","-1/6$\pi$", "0","1/6$\pi$","1/4$\pi$","1/3$\pi$","$\pi$"])
ax1.set_ylabel("$\\theta$")
ax2.set_yticks(np.array([-1,-(1/2)*np.sqrt(3),-(1/2)*np.sqrt(2),-(1/2)*np.sqrt(1),0,(1/2)*np.sqrt(1),(1/2)*np.sqrt(2),(1/2)*np.sqrt(3),1]))
ax2.set_yticklabels(["$-$1","$-\sqrt{\\frac{3}{2}}$","$-\sqrt{\\frac{2}{2}}$","$-{\\frac{1}{2}}$", "0","${\\frac{1}{2}}$","$\sqrt{\\frac{2}{2}}$","$\sqrt{\\frac{3}{2}}$","1"])
ax2.set_ylabel("X, Y")
ax1.set_xlabel("T")
ax1.set_xlim(t_span[0],t_span[-1])
ax1.set_ylim(-(1/3)*np.pi,(1/3)*np.pi)
#凡例
h1, l1 = ax1.get_legend_handles_labels()
h2, l2 = ax2.get_legend_handles_labels()
ax1.legend(h1+h2, l1+l2, ncol=3, fontsize=12)
plt.show()
#保存
#fig.savefig("XXX.png", format="png", dpi=330)
単振り子の軌跡
%matplotlib nbagg
G = 9.8 # 重力加速度(m/s^2)
L = 1.0 # 振り子の長さ(m)
# 運動方程式を解く
def eom(t, state):
dydt = np.zeros_like(state)
dydt[0] = state[1]
dydt[1] = -(G/L)*np.sin(state[0])
return dydt
#設定
t_span = [0,20] #時間 (s)
dt = 0.05# 時間の間隔 (s)
t = np.arange(t_span[0], t_span[-1], dt)# 時間の間隔のarray作成
#初期状態設定
th1 = 45.0# 角度 [deg]
w1 = 0.0 # 初期角速度 [deg/s]
state = np.radians([th1, w1])# 初期状態
#微分方程式を解く
sol = solve_ivp(eom, t_span, state, t_eval=t)
theta = sol.y[0] #y[0]＝θ、y[1]=ω
x = L * np.sin(theta)
y = -L * np.cos(theta)
#ーーーーーーグラフーーーーーーー
fig, ax = plt.subplots()
line, = ax.plot([], [], "o-", linewidth=2)
locus, = ax.plot([], [], 'r-', linewidth=2)
xlocus, ylocus = [], []
time_template = "time = %.1fs"
time_text = ax.text(0.05, 0.9, " ", transform=ax.transAxes)
def init():
line.set_data([], [])
time_text.set_text("")
return line, time_text
def animate(i):
thisx = [0, x[i]]
thisy = [0, y[i]]
line.set_data(thisx, thisy)
xlocus.append(x[i])
ylocus.append(y[i])
locus.set_data(xlocus, ylocus)
time_text.set_text(time_template % (i*dt))
return line,xlocus, ylocus, time_text
ani = FuncAnimation(fig, animate, frames=np.arange(0, len(t))
, interval=25, blit=True,init_func=init)
ax.set_xlim(-L-0.5,L+0.5)
ax.set_ylim(-L-0.5,L+0.5)
ax.set_aspect("equal")
ax.grid()
plt.show()
#ani.save("XXX.gif", writer="pillow", fps=20)
複数の単振り子
次に複数の単振り子です。
%matplotlib nbagg
G = 9.8 # 重力加速度(m/s^2)
L = 1.0 # 振り子の長さ(m)
# 運動方程式を解く
def eom(t, state):
dydt = np.zeros_like(state)
dydt[0] = state[1]
dydt[1] = -(G/L)*np.sin(state[0])
return dydt
#設定
t_span = [0,20] #時間 (s)
dt = 0.05# 時間の間隔 (s)
t = np.arange(t_span[0], t_span[-1], dt)# 時間の間隔のarray作成
#初期状態設定
th1s = np.arange(40.0,90.0,10.0)
w1 = 0.0 # 初期角速度 [deg/s]
states=[]
for th1 in th1s:
state = np.radians([th1, w1])# 初期状態
states.append(state)
thetas=[]
for state in states:
sol = solve_ivp(eom, t_span, state, t_eval=t)
theta = sol.y[0] #y[0]＝θ、y[1]=ω
thetas.append(theta)
xs,ys=[],[]
for theta in thetas:
x = L * np.sin(theta)
y = -L * np.cos(theta)
xs.append(x)
ys.append(y)
#ーーーーーーグラフーーーーーーー
fig, ax = plt.subplots()
lines=[]
for i in range(len(th1s)):
line, = ax.plot([], [], 'o-', lw=2, label=th1s[i])
lines.append(line,)
time_template = "time = %.1fs"
time_text = ax.text(0.05, 0.9, " ", transform=ax.transAxes)
def init():
for line in lines:
line.set_data([], [])
time_text.set_text("")
return line, time_text
def animate(i):
for j in range(len(lines)):
line = lines[j]
thisx = [0, xs[j][i]]
thisy = [0, ys[j][i]]
line.set_data(thisx, thisy)
time_text.set_text(time_template % (i*dt))
return line, time_text
ani = FuncAnimation(fig, animate, frames=np.arange(0, len(t))
, interval=25, blit=True,init_func=init)
ax.set_xlim(-L-0.5,L+0.5)
ax.set_ylim(-L-0.5,L+0.5)
ax.legend(title="init $\\theta$")
ax.set_aspect("equal")
ax.grid()
plt.show()
3ani.save("XXX.gif", writer="pillow", fps=20)
次は二重振り子によるカオスをアニメーションで描写したいと思います！
それでは🌏
【Matplotlib】アニメーション(FuncAnimation)
こんにちは
今日はMatplotlibのFuncAnimationを使ってアニメーションを作成したいと思います。
グラフ作成
Step1
まずはドキュメントの3個目のsin関数を描写してみましょう。
Step2
ドキュメントを参考にして
「y=ax²」のグラフを描いてみます。
モジュール
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
グラフ
%matplotlib nbagg
fig, ax = plt.subplots()
#y = ax^2
xdata, ydata = [], []
line, = plt.plot([], [], "b-") #*1
a = 1/2
def init():
ax.set_title("Basic")
ax.set_xlim(-10,10)
ax.set_ylim(-100,100)
#*2
ax.axhline(y=0, xmin=-10, xmax=10, color="k", linestyle="--")
ax.axvline(x=0, ymin=-100, ymax=100, color="k", linestyle="--")
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.grid() #*3
return line,
def update(frame):
#y=ax^2
xdata.append(frame)
ydata.append(a*(frame)*(frame))
line.set_data(xdata, ydata)
return line,
ani = FuncAnimation(fig, update, frames=np.linspace(-10, 10, 40),
init_func=init, blit=True) #*4
plt.show()
#ani.save("XXX.gif", writer="pillow", fps=15)
Keys
1) "b-"⇒colorが「b(blue)」でlineが「-」です
2) axhline(y=Y, xmin=最小値, xmax=最大値)⇒y=Yで直線を引く
axvline(x=X, ymin=最小値, ymax=最大値)⇒x=Xで直線を引く
3) 図にグリッド欄を拭く
4) np.linspace(a, b, n)⇒初項a、末項bで、項数nの等差数列を作成
(e.g.)
np.linspace(0,20,20)
⇒array([ 0. , 1.05263158, 2.10526316, 3.15789474, 4.21052632,
5.26315789, 6.31578947, 7.36842105, 8.42105263, 9.47368421,
10.52631579, 11.57894737, 12.63157895, 13.68421053, 14.73684211,
15.78947368, 16.84210526, 17.89473684, 18.94736842, 20. ])
Step3
最後に複数の線をプロットしてみます
モジュール
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
グラフ
%matplotlib nbagg
fig, ax = plt.subplots()
#sin
xdata, ydata = [], []
ln, = plt.plot([], [], "b-")
#cos
xdata2, ydata2 = [], []
ln2, = plt.plot([], [], "r-")
#tan
xdata3, ydata3 = [], []
ln3, = plt.plot([], [], "go",alpha=0.5)
def init():
ax.set_title("SIN & COS & TAN")
ax.set_xlim(-2*np.pi, 2*np.pi)
ax.set_ylim(-1, 1)
ax.axhline(y=0, xmin=-2*np.pi, xmax=2*np.pi, color="k", linestyle="--")
ax.axvline(x=0, ymin=-1, ymax=1, color="k", linestyle="--")
ax.set_xticks(np.arange(-2*np.pi,2*np.pi+0.01,(1/2)*np.pi)) #*2
ax.set_xticklabels(["-2$\pi$","-3/2$\pi$","-$\pi$","-1/2$\pi$", "0","1/2$\pi$","$\pi$","3/2$\pi$", "2$\pi$"]) #*1, 2
ax.set_xlabel("$\\theta$") #*2
ax.grid()
return ln,
def update(frame):
#sin
xdata.append(frame)
ydata.append(np.sin(frame))
ln.set_data(xdata, ydata)
#cos
xdata2.append(frame)
ydata2.append(np.cos(frame))
ln2.set_data(xdata2, ydata2)
#tan
xdata3.append(frame)
ydata3.append(np.tan(frame))
ln3.set_data(xdata3, ydata3)
return ln,ln2,ln3
ani = FuncAnimation(fig, update, frames=np.linspace(-2*np.pi, 2*np.pi, 300),
init_func=init, blit=True)
plt.show()
# ani.save("XXX.gif", writer="pillow", fps=15)
Keys
1)ラベルをラジアンに変換
2)ギリシャ文字
いかがだったでしょうか。
それでは🌏
【Matplotlib】2軸のグラフを作成しよう
こんにちは！2軸のグラフ作成をします
データの用意ができていて、図の書き方だけ知りたいという方は、描写 2軸へ！！
モジュール
必要なモジュールをインポートします
import numpy as np
import matplotlib.pyplot as plt
データ
それではまず初めにデータを用意します。
今回はランダムなデータを使用します。
x = np.arange(0,100,1)
y1 = np.random.rand(100)
y2 = np.random.rand(100) * 100
描写
1軸
2軸の描写の前に、それぞれのデータをプロットして確認しておきましょう。
fig = plt.figure(figsize=(10, 5))
plt.rcParams["font.size"] = 18 #fontsize一括設定
ax1 = plt.subplot(211)
ax2 = plt.subplot(212)
ax1.plot(x, y1)
ax2.plot(x, y2)
#軸設定
ax1.set_ylabel("Y1-axis")
ax2.set_ylabel("Y2-axis")
ax1.label_outer()
plt.show()
#保存
fig.savefig("XXX.png" ,format="png", dpi=330) 
2軸
それでは本題の2軸グラフです。
今回の目的は2軸グラフの描写なので、細かい設定はしません。
fig = plt.figure(figsize=(10, 3))
plt.rcParams["font.size"] = 18
ax1 = plt.subplot(111)
ax2 = ax1.twinx() #*1
ax1.plot(x, y1, color="b", label = "Y1")
ax2.plot(x, y2, color="r", label = "Y2")
#軸設定
ax1.set_ylabel("Y1-axis")
ax2.set_ylabel("Y2-axis")
#範囲設定
ax1.set_xlim(0,100)
ax1.set_ylim(0,1.2)
ax2.set_ylim(0,120)
#*2
h1, l1 = ax1.get_legend_handles_labels()
h2, l2 = ax2.get_legend_handles_labels()
ax1.legend(h1+h2, l1+l2, ncol=2, fontsize=12)
plt.show()
fig.savefig("XXX.png" ,format="png", dpi=330)
Keys
1) ax2は、ax1に.twinx()で結合します
2)凡例(legend)の表示です。
通常に
ax1.legend()
ax2.legend()としても良いのですが、それぞれが独立してしまうため表示が不自然になります。
したがって、.get_legend_handles_labels()から得られたhとlを足し算することによって統一の凡例が表示できます。
参考
ちなみにbarで違うデータを表示させると以下のようになります。
このようにすると、2つのデータの対比がしやすくなりますね！
それでは 🌏
【Science】水平線までどのくらいの距離？Part2ーマップで表示
こんにちは！
前回に引き続き、Cartopyでマップに表示したいと思います
[前回の記事]
水平線までどのくらいの距離なのか？
今回はCartopyを使ってどこまで見られるのか表示してみました！
都道府県などの境を表示するためにGADMのシェープファイルを利用させていただきました。
[ シェープファイル利用(Japan⇒Shapefile)]
地図
大気スケール
まずは大気スケールです。
日本列島の端から端まではだいたい2000 kmです
したがって、半径が1000 kmだとだいたい日本列島全体が見れますね
飛行機の高度の10 kmくらいの時はいちばん小さい円です。
それでも関東地方全体は軽く眺められますね。
| スケール | H | D | |
| 名前 | m | km | |
| 大気 | 熱圏 | 10000 | 1133.23 |
| 中間層 | 8000 | 1012.8 | |
| 成層圏 | 5000 | 799.75 | |
| 対流圏 | 1000 | 357.099 | |
山スケール
次に山スケールです。
今回は富士山とエベレストのみにします。(他はそこまで高さが変わらないので)
赤の点は富士山頂です。
富士山とエベレストで見られる範囲を比べると100 km程しか変わらないのでそこまで違いが見られませんね。それにしても、山に登るだけでここまで見れるのはすごいですよね。
建物スケール
次に建物スケールです。
赤い点は東京スカイツリーです。(634 m)
世界1位のブルジュハリファとスカイツリーを比較してもあまり見える距離は変わらないんですね。意外です。
しかし、東京タワーと比べると差が大きくなりますね。
都内の大きな高層マンションでは水平距離で24区見れそうですね、障害物がなければですが(笑
| スケール | H | D | |
| 名前 | m | km | |
| 建物 | ブルジュ ハリファ | 828 | 102.7 |
| 東京スカイツリー | 634 | 89.9 | |
| 東京タワー | 333 | 65.1 | |
| 福岡タワー | 234 | 54.6 | |
| 40階 | 150 | 43.7 | |
| 20階 | 60 | 27.7 | |
| 10階 | 30 | 19.6 | |
| 5階 | 12 | 12.4 | |
| 3階 | 9 | 10.7 | |
| 2階 | 6 | 8.7 | |
人スケール
最後に人スケールです。
平均身長あたりだと150ー170 cmなので見られる距離は4キロほどです。
東京ディズニーリゾートのある舞浜からです。
実際東京湾に立つと、羽田空港から千葉県が見られますし、千葉県から東京や川崎が見られます。
それは、水平線ではなく建物の高さがあるためや地上の高低差ならびに大気の屈折があるためと考えられます。
| スケール | H | D | |
| 名前 | m | km | |
| 身長 | 10cm | 0.1 | 1.1 |
| 30cm | 0.3 | 2 | |
| 50cm | 0.5 | 2.5 | |
| 70cm | 0.7 | 3 | |
| 90cm | 0.9 | 3.4 | |
| 110cm | 1.1 | 0.7 | |
| 130cm | 1.3 | 4.1 | |
| 150cm | 1.5 | 4.4 | |
| 170cm | 1.7 | 4.7 | |
| 190cm | 1.9 | 4.9 | |
| 210cm | 2.1 | 5.2 | |
| 230cm | 2.3 | 5.4 | |
| 250cm | 2.5 | 5.6 | |
| 270cm | 2.7 | 5.9 |
コード
まずは描写に必要なモジュールをインポートします
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
#shape file
import cartopy.io.shapereader as shpreader
次に、地球状のどこかかの点から円を描くので
半径(r km)と中心の緯度経度から、円の緯度経度座標に変換します
# 円
def circle(r,S,T):
#r:radius(km), (S,T):centerDegree
s,t = S*111, T*111
points = 1000
theta = np.linspace(0, 2 * np.pi, points)
x = r * np.sin(theta) + s
y = r * np.cos(theta) + t
#x,y=>Degree
X, Y = x/111, y/111
return (X,Y)Keys
「km⇄緯度経度」は赤道半径を利用した値で変換しました
km = Degree ✕ 111
それでは描写してみます
def main():
fig = plt.figure(figsize=(5,5))
plt.rcParams["font.size"] = 18
# TokyoStation
ts_lon = 139.767125
ts_lat = 35.681236
ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree())
ax.set_extent([ts_lon-1, ts_lon+1, ts_lat-1, ts_lat+1], crs=ccrs.PlateCarree())
ax.set_title("Title")
ax.add_feature(cfeature.BORDERS, linestyle=":")
#都道府県境と市町村境
fname1 = "data/gadm36_JPN_shp/gadm36_JPN_1.shp"
fname2 = "data/gadm36_JPN_shp/gadm36_JPN_2.shp"
adm1_shapes = list(shpreader.Reader(fname1).geometries())
adm2_shapes = list(shpreader.Reader(fname2).geometries())
#データのプロット
ax.add_geometries(adm1_shapes, ccrs.PlateCarree(),
edgecolor="k",facecolor="w",alpha=0.5)
ax.add_geometries(adm2_shapes, ccrs.PlateCarree(),
edgecolor="gray", facecolor="w", alpha=0.5)
#TokyoStation
ax.scatter(ts_lon,ts_lat,s=20,c="r",alpha=1, zorder=10)
c = circle(100, ts_lon,ts_lat)
ax.plot(c[0],c[1],c="b",alpha=1, zorder=10)
plt.show()
if __name__ == "__main__":
main()
参考コード
それでは🌏
関連記事
【Science】水平線までどのくらいの距離？Part1
こんにちは！
皆さん水平線までの距離ってどのくらいなんだろうって考えたことありませんか？
今回はそれを計算してみました。
計算
半径をrとし、地表面からの高さをhとします。
そして地球(あるいは天体)を完全な球体と仮定し、hと円との接線を引き、そこまでの距離をdとします。三平方の定理より、
⇔
∵d>0
※簡易的な計算であるため、dまでには障害物がないと仮定し、大気による光の屈折はしないと仮定します。地域で気温や大気質に違いがあるため、考慮しないです。
※実際見える対象物は平面ではなく、高さがあるのでdまでが見える範囲とは限りません。
高さと距離の関係
試しにエベレストで距離dを調べてみました。エベレストの標高は8848 mですので、
地球の赤道半径6378 kmを利用すると、約340 kmでした。
意外と見れますね。
これを踏まえていろいろな高さから見られる距離を求めてみました。
グラフ
大気スケール
上から熱圏、中間層、成層圏、対流圏です。
高度が高くなればなるほど見える距離は大きくなりますが、その変化の幅は小さくなります。
それはdの式がh²のルートに比例するからでしょう。
山スケール
世界の山ランキングと日本の山ランキングからのデータをプロットしました。
世界トップ10の山はどれも8000 m級で、それほど差が無いので見える距離もほどんど変わらず、300ー350 km
日本のトップテンも同様に3000 m急なので、200 km
という結果になりました。
Highest Mountain Peaks of the World
建物スケール
次はさらに低い建物で見ていきます。
世界一のブルジュハリファだと102 m
日本一の東京スカイツリーだと90 m
でした。
他にも1階あたりを3 mとして、40 m, 20 m, 10m, 5m, 3m, 2mと計算してみました
人スケール
そして、身長によってどのくらい見える距離が変わるのかを調べてみました。
身長170 cmのとき4.7 km
身長150 cmのとき4.4 km
身長が20 cm違っても300 m位しか差がありませんね(笑
太陽系
先ほどまでは地球の半径で計算してきましたが、次に太陽系の惑星の半径で計算してみました。
高さは人スケールにしました。
太陽、水星、金星、地球、火星、木星、土星、天王星、海王星、冥王星、エリス、月
太陽は圧倒的に半径が大きいので差が明瞭に見られますね。
まとめ
人はだいたい4キロくらいですね
| スケール | H | D | スケール | H | D | ||
| 名前 | m | km | 名前 | m | km | ||
| 大気 | 熱圏 | 10000 | 1133.23 | 建物 | ブルジュ ハリファ | 828 | 102.7 |
| 中間層 | 8000 | 1012.8 | 東京スカイツリー | 634 | 89.9 | ||
| 成層圏 | 5000 | 799.75 | 東京タワー | 333 | 65.1 | ||
| 対流圏 | 1000 | 357.099 | 福岡タワー | 234 | 54.6 | ||
| 山(世界) | エベレスト | 8850 | 335.9 | 40階 | 150 | 43.7 | |
| K2 | 8611 | 331.4 | 20階 | 60 | 27.7 | ||
| カンチェンジュンガ | 8586 | 330.9 | 10階 | 30 | 19.6 | ||
| ローツェ | 8516 | 329.5 | 5階 | 12 | 12.4 | ||
| マカルー | 8463 | 328.5 | 3階 | 9 | 10.7 | ||
| チョ・オユ－ | 8201 | 323.4 | 2階 | 6 | 8.7 | ||
| ダウラギリ | 8167 | 322.7 | 身長 | 10cm | 0.1 | 1.1 | |
| マナスル | 8163 | 322.6 | 30cm | 0.3 | 2 | ||
| ナンガ・パルバット | 8125 | 321.9 | 50cm | 0.5 | 2.5 | ||
| アンナプルナ | 8091 | 321.2 | 70cm | 0.7 | 3 | ||
| 山(日本) | 富士山 | 3776 | 219.4 | 90cm | 0.9 | 3.4 | |
| 北岳 | 3193 | 201.7 | 110cm | 1.1 | 0.7 | ||
| 奥穂高岳 | 3190 | 201.6 | 130cm | 1.3 | 4.1 | ||
| 間ノ岳 | 3190 | 201.6 | 150cm | 1.5 | 4.4 | ||
| 槍ヶ岳 | 3189 | 201.3 | 170cm | 1.7 | 4.7 | ||
| 東岳 | 3141 | 200.1 | 190cm | 1.9 | 4.9 | ||
| 赤石岳 | 3121 | 199.4 | 210cm | 2.1 | 5.2 | ||
| 涸沢岳 | 3110 | 199.1 | 230cm | 2.3 | 5.4 | ||
| 北穂高岳 | 3106 | 199 | 250cm | 2.5 | 5.6 | ||
| 大喰岳 | 3101 | 198.8 | 270cm | 2.7 | 5.9 | ||
是非皆さんも半径や高度を変えていろいろと調べてみてください！
それでは🌏
【Jupyter Notebook】記号
こんにちは！
科学(気象系)で主に使う記号を紹介します。
LaTeXと少し違って、「\」を2つ付けるものもあるので注意しましょう。
| Symbol_1 | Symbol_2 | ||||
| Symbol Name | Letter | Jupyter | Symbol Name | Letter | Jupyter |
| イコール | = | = | 大なり | > | > |
| ノットイコール | ≠ | $\\neq$ | 小なり | < | < |
| ニアリーイコール | ≃ | $\simeq$ | 十分大 | ≫ | $\gg$ |
| プラス | ＋ | + | 十分小 | ≪ | $\ll$ |
| マイナス | － | $-$ | 以上 | ≥ | $\geq$ |
| 掛ける | × | $\\times$ | 以上 | ≧ | $\geqq$ |
| 割る | ÷ | $\div$ | 以下 | ≤ | $\leq$ |
| プラスマイナス | ± | $\pm$ | 以下 | ≦ | $\leqq$ |
| マイナスプラス | ∓ | $\mp$ | 無限 | ∞ | $\infty$ |
Keys
neq = not equal
simeq = similar equal
div = division
pm = plus minus
mp = minus plus
geq = greater or equal
leq = less or equal
infty = infinity
参考(LaTeX)
それでは🌏
類似記事
【Jupyter Notebook】ギリシャ文字
こんにちは
今回は、Jupyternotebookでグラフ作成時によく使う
ギリシャ文字についてです。
小文字のαやβなどの文字は、「\」を二つ(「\\」)付けないといけないので
注意しましょう。
| Capital | Small | |||
| Greek Alphabet | Letter | Jupyter | Letter | Jupyter |
| アルファ | Α | $A$ | α | $\\alpha$ |
| ベータ | B | $B$ | β | $\\beta$ |
| ガンマ | Γ | $\Gamma$ | γ | $\gamma$ |
| デルタ | Δ | $\Delta$ | δ | $\delta$ |
| イプシロン | Ε | $E$ | ε | $\epsilon$ |
| ゼータ | Z | $Z$ | ζ | $\zeta$ |
| エータ | Η | $H$ | η | $\eta$ |
| シータ | Θ | $\\Theta$ | θ | $\\theta$ |
| イオタ | Ι | $I$ | ι | $\iota$ |
| カッパ | Κ | $K$ | κ | $\kappa$ |
| ラムダ | Λ | $\Lambda$ | λ | $\lambda$ |
| ミュー | M | $M$ | μ | $\mu$ |
| ニュー | N | $N$ | ν | $\\nu$ |
| クサイ | Ξ | $\\Xi$ | ξ | $\\xi$ |
| オミクロン | Ο | $O$ | ο | $\o$ |
| パイ | Π | $\Pi$ | π | $\pi$ |
| ロー | Ρ | $P$ | ρ | $\\rho$ |
| シグマ | Σ | $\Sigma$ | σ | $\sigma$ |
| タウ | Τ | $T$ | τ | $\\tau$ |
| ウプシロン | Υ | $\\Upsilon$ | υ | $\\upsilon$ |
| ファイ | Φ | $\Phi$ | φ | $\phi$ |
| キー | Χ | $X$ | χ | $\chi$ |
| プサイ | Ψ | $\Psi$ | ψ | $\psi$ |
| オメガ | Ω | $\Omega$ | ω | $\omega$ |
参考
それでは🌏
類似記事
【Matplotlib】Gridspecーグラフの柔軟な分割(3:1でもで4:1でも)
こんにちは。
グラフの分割をすると思うのですが、普通に分割すると1:1とか一つ一つが同じ大きさのグラフになりますよね。
そこで今回はMatplotlibのGridspecを用いて、3:1やさらに複雑なグラフの分割をしていきたいと思います！
★事前に必要な知識
グラフの基本的な書き方(matplotlib)
リストのスライス
基礎
まずは今回使うGridspecをインポートしましょう・
import matplotlib.pyplot as plt
from matplotlib import gridspec
例1
縦方向に3:1で分割してみます
fig = plt.figure(figsize=(4,3))
#分割 *1
ncols=2
nrows=1
gs = gridspec.GridSpec(ncols, nrows, height_ratios=(3, 1))
ax1 = fig.add_subplot(gs[:1, 0])
ax2 = fig.add_subplot(gs[-1, 0])
ax1.text(0.5,0.5, "ax1", fontsize=30,
horizontalalignment='center',
verticalalignment='center',)
ax2.text(0.5,0.5, "ax2", fontsize=30,
horizontalalignment='center',
verticalalignment='center',)
#グラフ間の隙間
plt.subplots_adjust(hspace=0.25)
fig.savefig("XXX.png", format='png', dpi=330)
Key
1) gridspec.GridSpec(縦方向分割数, 横方向分割数, height_ratios=(3, 1))
height_ratioはその比率です。
例2
横方向に1:3で分割してみます
fig = plt.figure(figsize=(4,3))
ncols=1
nrows=2
gs = gridspec.GridSpec(ncols, nrows, width_ratios=(1,3)) #*1
ax1 = fig.add_subplot(gs[0, 0])
ax2 = fig.add_subplot(gs[0, 1:])
ax1.text(0.5,0.5, "ax1", fontsize=30,
horizontalalignment='center',
verticalalignment='center',)
ax2.text(0.5,0.5, "ax2", fontsize=30,
horizontalalignment='center',
verticalalignment='center',)
#グラフ間の隙間
plt.subplots_adjust(wspace=0.25)
fig.savefig("XXX.png", format='png', dpi=330)
Key
1) gridspec.GridSpec(縦方向分割数, 横方向分割数, width_ratios=(1, 3))
width_ratioはその比率です。
応用
さらに柔軟にグラフを描いてみましょう
例3
縦方向に1:3で、横方向に3:1に分割してみましょう
fig = plt.figure(figsize=(6, 6))
ncols=2
nrows=2
gs = gridspec.GridSpec(ncols, nrows, height_ratios=(1, 3), width_ratios=(3, 1))
ax1 = fig.add_subplot(gs[:1, :1])
ax2 = fig.add_subplot(gs[-1, :1])
ax3 = fig.add_subplot(gs[:1, -1])
ax4 = fig.add_subplot(gs[-1, -1])
ax1.text(0.5,0.5, "ax1", fontsize=30,
horizontalalignment='center',
verticalalignment='center',)
ax2.text(0.5,0.5, "ax2", fontsize=30,
horizontalalignment='center',
verticalalignment='center',)
ax3.text(0.5,0.5, "ax3", fontsize=30,
horizontalalignment='center',
verticalalignment='center',)
ax4.text(0.5,0.5, "ax4", fontsize=30,
horizontalalignment='center',
verticalalignment='center',)
#グラフ間の隙間
plt.subplots_adjust(hspace = 0.20, wspace=0.25)
fig.savefig("XXX.png", format='png', dpi=330)
例4
次は比をばらばらにして分割してみましょう
fig = plt.figure(figsize=(6, 6))
ncols=3
nrows=2
gs = gridspec.GridSpec(ncols, nrows, height_ratios=(1, 1, 1), width_ratios=(3, 1))
ax1 = fig.add_subplot(gs[:2, :1])
ax2 = fig.add_subplot(gs[-1, :1])
ax3 = fig.add_subplot(gs[:1, -1])
ax4 = fig.add_subplot(gs[1:, -1])
ax1.text(0.5,0.5, "ax1", fontsize=30,
horizontalalignment='center',
verticalalignment='center',)
ax2.text(0.5,0.5, "ax2", fontsize=30,
horizontalalignment='center',
verticalalignment='center',)
ax3.text(0.5,0.5, "ax3", fontsize=30,
horizontalalignment='center',
verticalalignment='center',)
ax4.text(0.5,0.5, "ax4", fontsize=30,
horizontalalignment='center',
verticalalignment='center',)
#グラフ間の隙間
plt.subplots_adjust(hspace = 0.3,wspace=0.25)
fig.savefig("XXX.png", format='png', dpi=330)
例5
さらに細かく変形してみましょう
fig = plt.figure(figsize=(6, 6))
ncols=3
nrows=3
gs = gridspec.GridSpec(ncols, nrows, height_ratios=(1, 1, 1), width_ratios=(2, 2, 1))
ax1 = fig.add_subplot(gs[:2, :2])
ax2 = fig.add_subplot(gs[-1, :1])
ax3 = fig.add_subplot(gs[-1, 1:2])
ax4 = fig.add_subplot(gs[:1, -1])
ax5 = fig.add_subplot(gs[1:, -1])
ax1.text(0.5,0.5, "ax1", fontsize=30,
horizontalalignment='center',
verticalalignment='center',)
ax2.text(0.5,0.5, "ax2", fontsize=30,
horizontalalignment='center',
verticalalignment='center',)
ax3.text(0.5,0.5, "ax3", fontsize=30,
horizontalalignment='center',
verticalalignment='center',)
ax4.text(0.5,0.5, "ax4", fontsize=30,
horizontalalignment='center',
verticalalignment='center',)
ax5.text(0.5,0.5, "ax5", fontsize=30,
horizontalalignment='center',
verticalalignment='center',)
#グラフ間の隙間
plt.subplots_adjust(hspace = 0.3, wspace=0.4)
fig.savefig("XXX.png", format='png', dpi=330)
いかがだったでしょうか。
Grodspecを用いると、柔軟にグラフを分割することが可能です。
是非試してみてください！
それでは🌏
