【Matplotlib】円と楕円の描画

円（Circle）の描画

原点 $\mathrm{O}$ に中心をもつ円周上の点 $\mathrm{P}$ の座標は三角関数を使って、$(r\cos\theta,\ r\sin\theta)$ で表されることが知られています。$r$ は円の半径、$\theta$ は 線分 $OP$ が $x$ 軸となす角度です。中心が $(x_c,\ y_c)$ にあるときは、円周上の点は次のように表せます。
 \[(x_c+r\cos\theta,\ y_c+r\sin\theta)\]
NumPy で角度の連続データを作成して数学的に円をプロットすることができます。半径と角度、中心座標を指定して円を描く関数を定義してみます。

# PYTHON_DRAW_CIRCLE

# In[1]

import numpy as np
import matplotlib.pyplot as plt

# 円をプロットする関数
# radiusで半径、nで角度360°の分割数を指定
# center_markerにTrueを指定したときは円の中心点を表示する
def draw_circle(ax, radius, center=(0,0), n=129,
                color="blue", center_marker=True):

    # 角度θの刻みデータ
    theta = np.linspace(0, 2*np.pi, n)

    # 円周上のx座標とy座標
    x = center[0] + radius * np.cos(theta)
    y = center[1] + radius * np.sin(theta)

    # x軸とy軸のラベルを表示
    ax.set_xlabel("x", fontsize = 16)
    ax.set_ylabel("y", fontsize = 16)

    # データをプロット
    ax.plot(x, y, color=color)
    
    # 原点Oをプロット
    ax.scatter(0, 0, color="black")

    # 円の中心座標をプロット
    if center_marker == True:
        ax.scatter(center[0], center[1], color=color)

draw_circle() 関数を使って、中心座標 $(1,\ 1)$、半径 $3$ の円を描いてみます。

# In[2]

# 描画領域の設定
fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111)
ax.grid()
ax.set_title("Circle", fontsize = 16)
ax.set_xlim(-5, 5)
ax.set_ylim(-5, 5)

# 中心(1,1),半径3の円を描く
draw_circle(ax, radius=3, center=(1, 1))

draw_circle() を使って単位円 (半径 $1$ の円) を描き、いくつかラベルを添えてみます。コードライブラリ の pointer() 関数の実装が前提となります。

# in[3]

# ここにpointer()関数を実装する

# 描画領域の設定
fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111)
ax.grid()
ax.set_title("Unit Circle", fontsize = 16)
ax.set_xlabel("x", fontsize = 16)
ax.set_ylabel("y", fontsize = 16)
ax.set_xlim(-1.5, 1.5)
ax.set_ylim(-1.5, 1.5)

# 角度t=50°におけるx座標とy座標
px = np.cos(np.deg2rad(50))
py = np.sin(np.deg2rad(50))

# 単位円を描画
draw_circle(ax, radius=1, color="red", center_marker=False)

# 円周上の点P
pointer(ax, px, py, "P", angle=45,
        pcolor="red", textcolor = "red",
        textsize=16, pad=0.05)

# cosθ
pointer(ax, px, 0, "cosθ", angle=250,
        pcolor="black", textcolor = "black",
        textsize=16, pad=0.18)

# sinθ
pointer(ax, 0, py, "sinθ", angle=190,
        pcolor="black", textcolor = "black",
        textsize=16, pad=0.34)

# x軸とy軸に沿った直線
ax.plot([0,1],[0,0], color="gray", linewidth=1.5)
ax.plot([0,0],[0,1], color="gray", linewidth=1.5)

# 原点と単位円を結ぶ直線をプロット
ax.plot([0,px],[0,py], color="black")

# Pからx軸とy軸に下ろした直線を点線で表示
ax.plot([px,px],[0,py], color="black", linestyle="--")
ax.plot([0,px],[py,py], color="black", linestyle="--")

# 角度θを表記
ax.text(0.12, 0.03, "θ", size=15)

matplotlib.patches.Circle クラスは Axes に円形として描画されるインスタンス（オブジェクト）を作成します。引数には中心座標、半径、塗り潰しの色などを指定します。

# PYTHON_MATPLOTLIB_PATCHES_CIRCLE

# In[1]

# matplotlib.pyplotをインポート
import matplotlib.pyplot as plt

# matplotlib.patchesをインポート
import matplotlib.patches as pat

# Figureを作成
fig = plt.figure(figsize=(5, 5))

# FigureにAxes(サブプロット)を追加
ax = fig.add_subplot(111)

# matplotlib.patches.Circleのインスタンスを作成
# 中心位置の座標は左下を基準に全体の幅を1として設定
# radiusは半径、colorは塗り潰しの色
c = pat.Circle(xy = (0.5, 0.5), radius = 0.25, color = "blue")

# Axesに円cを追加
ax.add_patch(c)

中心位置の座標は Axes (サブプロット) の左下を基準に全体の幅を1として設定します。サンプルコードでは (0.5, 0.5) で指定しているので、Axes の中央に表示されることになります。同様に半径も Axes の大きさを 1 とした長さを指定します。サンプルコードでは radius = 0.25 なので、Axes の幅の 1/4 の長さの半径を指定していることになります。
 
次のサンプルコードは、黒い線で描かれるシンプルな円を表示させます。

# In[2]

# FigureとAxesを準備する
fig = plt.figure(figsize=(5, 5))
ax = fig.add_subplot(111)

# matplotlib.patches.Circleのインスタンスを作成
# xy座標は左下を基準に全体の幅を1として設定
# fcはface color(塗り潰しの色), ecはedge color（縁の色）
c = pat.Circle(xy = (0.5, 0.5), radius = 0.25,
               fc = "white", ec = "black")

# Axesに円cを追加
ax.add_patch(c)

fc は face color (塗り潰しの色), ec は edge color（縁の色）の略記です。サンプルコードでは fc に white (白)、ec に black (黒) を指定しています。fc = “white” の代わりに color = “white” と記述しても同じです。

楕円（Ellipse）の描画

楕円 は平面上の 2 定点からの距離の和が一定である点の集合として定義されています。2 定点のことを焦点 (focus) といいます。両端が楕円に接するように、楕円内部に 2 焦点を結ぶ直線を引いたとき、この線分を長軸とよび、その長さを長径、長径の半分の長さを長軸半径といいます。また、この長軸を二等分するように垂線を楕円内部に引くとき、この線分を短軸とよび、その長さを短径、短径の半分の長さを短軸半径とよびます。
 

上図において、長径は $2a$、長軸半径は $a$ です。また、短径は $2b$、短軸半径は $b$ です。楕円の重心 $O$ から 2 つの焦点 $F,\ F’$ の距離 $c$ は $\sqrt{a^2-b^2}$ で与えられます。直交座標系における楕円の代数方程式 (標準形) は
 \[\frac{x^2}{a^2}+\frac{y^2}{b^2}=1\tag{1}\]
で表されます。$a=b$ のときは、半径 $a$ の円の方程式 $x^2+y^2=a^2$ となるので、楕円は円を特別な形として含みます。直交座標系の楕円方程式を $F$ を原点とする極座標 $(r,\ \varphi)$ で書き直すと、
 \[r=\frac{l}{1+e\cos\varphi}\tag{2}\]
となります。$e$ は離心率、$l$ は半直弦とよばれる定数で、それぞれ以下のように定義されています。
 \[e=\frac{\sqrt{a^2-b^2}}{a},\quad \ l=\frac{b^2}{a}\tag{3}\]
離心率 $e$ は楕円の形を決定するパラメータです。後でコードを実行して確認したいと思います。プログラミングで実装しやすいのは、離心近点角 $\theta=\angle QOF$ を媒介変数とする
 \[x=a\cos\theta, \quad y=b\sin\theta\tag{4}\]
の形で表される方程式です。本格的なコード (楕円クラス) は後で実装するとして、最初は簡単なコードで楕円を描いてみます。

# PYTHON_DRAW_ELLIPSE

# In[1]

import numpy as np
import matplotlib.pyplot as plt

# 長軸半径と短軸半径
a = 2
b = 1

# 離心率
e = np.sqrt(a**2 - b**2) / a

# 離心近点角の連続データ
theta = np.linspace(0, 2*np.pi, 65)

# 楕円のデータ
x = a * np.cos(theta)
y = b * np.sin(theta)

# 描画領域の設定
fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111)
ax.grid()
ax.set_aspect('equal')
ax.set_title("Ellipse (e={:.3f})".format(e), fontsize=16)
ax.set_xlabel("x", fontsize=16)
ax.set_ylabel("y", fontsize=16)
ax.set_xlim(-2.5, 2.5)
ax.set_ylim(-1.5, 1.5)

# 楕円をプロット
ax.plot(x, y, color="blue")

plt.show()

第二種完全楕円積分と楕円の弧長

　楕円の弧長 (周の長さ) の計算式はかなり複雑です：
 \[L=4a\int_{0}^{\pi/2}\sqrt{1-e^2\sin^2\theta}d\theta\]
右辺の定積分
 \[E(e)=\int_{0}^{\pi/2}\sqrt{1-e^2\sin^2\theta}d\theta\]
は第二種完全楕円積分とよばれる難解な積分です。この積分は普通の方法で求めることはできないので、入り組んだ無限級数を収束するまで計算する必要がありますが、幸いにも SciPy の特殊関数パッケージには第二種完全楕円積分を求める ellipe() 関数が用意されています。引数には離心率を渡します。In[1] で描いた楕円 ($a=2,\ b=1$) の弧長を計算してみます。

# In[2]

from scipy.special import ellipe

# 楕円(a=2,b=1)の弧長
L = 4 * a * ellipe(e)

print("{:.3f}".format(L))
# 9.052

楕円の面積

楕円の面積は $\pi ab$ で与えられます。$a=b$ のときは $\pi a^2$ となって、半径 $a$ の円の面積を表します。In[1] で描いた楕円の面積を計算してみましょう。

# In[3]

#楕円(a=2,b=1)の面積
S = a * b * np.pi

print("{:.3f}".format(S))
# 6.283

楕円クラスの実装

楕円に関する各種データを保持するオブジェクトを生成するために、楕円クラス Ellipse を実装しておきましょう。

# In[4]

# 楕円クラスを定義
class Ellipse:
    # 円周率の近似値をクラス変数として定義
    pi = np.pi

    # インスタンス変数の定義
    def __init__(self, a, b, center=(0,0), n=129):
        self.major_radius = a  # 長軸半径
        self.minor_radius = b  # 短軸半径
        self.center = center  # 中心点の座標
        self.f_distance = np.sqrt(a**2 - b**2)  # 原点から焦点までの距離
        
        # 右側焦点の座標と左側焦点の座標
        self.focus_1 = (self.center[0] + self.f_distance, self.center[1])
        self.focus_2 = (self.focus_1[0] - 2 * self.f_distance, self.center[1])
        
        self.eccentricity = self.f_distance / a  # 離心率
        self.area =self.pi * self.major_radius * self.minor_radius  # 面積
        self.perimeter = 4 * ellipe(self.eccentricity)  # 周の長さ
        
        t = np.linspace(0, 2*np.pi, n)
        self.x = self.center[0] + self.major_radius * np.cos(t)
        self.y = self.center[1] + self.minor_radius * np.sin(t)

    # 回転メソッド
    # angleで回転角度(ラジアン)を指定
    # r_axisで回転軸の座標を指定
    def rotate(self, angle, r_axis):

        # 回転行列
        R = np.array([[np.cos(angle), -np.sin(angle)],
                      [np.sin(angle),  np.cos(angle)]])

        # 焦点の回転
        p = np.array(r_axis)
        self.focus_1 = np.dot(R, self.focus_1 - p) + p
        self.focus_2 = np.dot(R, self.focus_2 - p) + p

        # 楕円の回転
        q = np.array([[r_axis[0]],[r_axis[1]]])
        arr = np.vstack((self.x, self.y))
        arr_2 = np.dot(R, arr - q) + q
        self.x = arr_2[0]
        self.y = arr_2[1]

    # データ一覧表示メソッド
    def show_data(self):
        print("長軸半径　{}\n".format(self.major_radius),
              "短軸半径　{}\n".format(self.minor_radius),
              "中心点　{}\n".format(self.center),
              "第1焦点　{}\n".format(self.focus_1),
              "第2焦点　{}\n".format(self.focus_2),
              "離心率　{}\n".format(self.eccentricity),
              "面積　{}\n".format(self.area),
              "弧長　{}\n".format(self.perimeter))

    # 楕円の描画メソッド
    # oval_colorで楕円の線の色を指定
    # display_focusをTrueに設定すると焦点を表示
    # display_centerをTrueに設定すると中心点を表示
    def draw(self, ax, n=129, linestyle="-", oval_color="black",
             display_focus=False, focus_color="black",
             display_center=False, center_color="black"):

        ax.plot(self.x, self.y, color=oval_color, linestyle=linestyle)

        if display_focus == True:
            ax.scatter(self.focus_1[0], self.focus_1[1], color=focus_color)
            ax.scatter(self.focus_2[0], self.focus_2[1], color=focus_color)
        
        if display_center == True:
            ax.scatter(self.center[0], self.center[1], color=center_color)

座標 $(1,\ 1)$ を中心とする、長軸半径 $a=2$、短軸半径 $b=1$ の楕円オブジェクトを生成してみます。

# In[5]

# 楕円クラスのインスタンスを生成
oval = Ellipse(2, 1, center=(1, 1))

　楕円オブジェクトは自身のデータをインスタンス変数という形で常に保持しています (このへんがオブジェクト指向設計の便利なところです)。たとえば .eccentricity で楕円の離心率を取得できます。

# In[6]

# 楕円の離心率を取得
e = oval.eccentricity

print("{:.3f}".format(e))
# 0.866

　楕円のデータ一覧を表示するメソッドも備えておきました。

# In[7]

# 楕円のデータ一覧
oval.show_data()

'''
長軸半径　2
短軸半径　1
中心点　(1, 1)
第1焦点　(2.732050807568877, 1)
第2焦点　(-0.7320508075688772, 1)
離心率　0.8660254037844386
面積　6.283185307179586
弧長　4.52587616787211
'''

　draw() メソッドで平面座標上に楕円を描画できます。

# In[8]

# 描画領域の設定
fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111)
ax.grid()
ax.set_aspect('equal')
ax.set_title("Ellipse", fontsize=16)
ax.set_xlabel("x", fontsize=16)
ax.set_ylabel("y", fontsize=16)
ax.set_xlim(-4, 4)
ax.set_ylim(-3, 3)

# 楕円をプロット
oval.draw(ax, display_focus=True, display_center=True)

plt.show()

デフォルトでは長軸を $x$ 軸、短軸を $y$ 軸に合わせています。rotate() メソッドで楕円を任意の中心点の軸周りに回転させることができます。

# In[9]

# 楕円を(1,1)を軸として反時計回りにpi/4だけ回転させる
oval.rotate(np.pi/4, (1, 1))

# 回転した楕円をプロット
oval.draw(ax, display_focus=True, display_center=True,
          oval_color="red", center_color="red", focus_color="red")

# Figureの再表示
display(fig)

離心率 $e=\sqrt{a^2-b^2}/a$ は楕円の形を決定するパラメータです。$a=b$ のときは $e=0$ なので、円は離心率 $0$ の楕円といえます。$a$ と $b$ の差が大きくなるほど、楕円の形は細長くなっていきます。短軸半径 $b$ を $1$ に固定して、長軸半径 $a$ を $1$ から $4$ まで変化させて、楕円の形の変化を眺めてみましょう。

# In[10]

# 描画領域の設定
fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111)
ax.grid()
ax.set_aspect('equal')
ax.set_title("Ellipse", fontsize = 16)
ax.set_xlabel("x", fontsize = 16)
ax.set_ylabel("y", fontsize = 16)
ax.set_xlim(-6, 6)
ax.set_ylim(-6, 6)

# 長軸半径aを変化させながらプロット
for k in range(1, 5):
    oval = Ellipse(k, 1, center=(0, 6.5-2.5*k))
    oval.draw(ax)
    ax.text(-0.75, 6.4-2.5*k,
            "e={:.2f}".format(oval.eccentricity), size=12, color="red")

matplotlib.patches.Ellipse

Matplotlib には楕円オブジェクトを生成する matplotlib.patches.Ellipse クラスが用意されています。引数には中心座標、楕円の横幅と縦幅、水平軸からの回転角、塗り潰しの色などを指定します。

# PYTHON_MATPLOTLIB_PATCHES_ELLIPSE

import matplotlib.pyplot as plt
import matplotlib.patches as pat

# FigureとAxes
fig = plt.figure(figsize=(5, 5))
ax = fig.add_subplot(111)

# matplotlib.patches.Ellipsクラスのインスタンスを作成
# 水平軸からの回転角(angle)を 0, 45, 90 に設定
e1 = pat.Ellipse(xy=(0.5, 0.5), width=0.8, height=0.4, alpha=0.6,
                 angle=0, color="red", label="angle 00")

e2 = pat.Ellipse(xy=(0.5, 0.5), width=0.8, height=0.4, alpha=0.6,
                 angle=45, color="blue", label="angle 45")

e3 = pat.Ellipse(xy=(0.5, 0.5), width=0.8, height=0.4, alpha=0.6,
                 angle=90, color="green", label="angle 90")

# Axesに楕円を追加
ax.add_patch(e1)
ax.add_patch(e2)
ax.add_patch(e3)

# 凡例の表示位置を左上に設定
ax.legend(loc = "upper left")

コード PYTHON_MATPLOTLIB_PATCHES_ELLIPSE では、angle (回転角) の異なる 3 つの楕円を作成して重ねています。左上の凡例に、それぞれの図形に指定した angle が表示されています。実行結果をみると、angle は水平軸から反時計回りの角度で指定することがわかります。また、alpha引数を指定すると図形を半透明で表示させることができます。

扇形（くさび型）の描画

matplotlib.patches.Wedge クラスから扇形を作ることができます。
英語で wedge は「くさび型」を意味しますが、「円をくさび型に切る」ことで扇形を作ることから、このようなクラス名がついています。切り取り位置（角度）は引数 theta1 と theta2 で指定します。

# PYTHON_MATPLOTLIB_PATCHES_WEDGE

import matplotlib.pyplot as plt
import matplotlib.patches as pat

# FigureとAxes
fig = plt.figure(figsize=(5, 5))
ax = fig.add_subplot(111)

# patches.Wedgeクラスのインスタンスを作成
# 中心座標(0.5, 0.5), 半径0.4, 切込み位置0°, 60°
w = pat.Wedge(center=(0.5, 0.5), r=0.4,
             theta1=0, theta2=60, color="green")

# Axesに扇形を追加
ax.add_patch(w)

theta1 と theta2 は水平軸から反時計回りの角度で指定します。サンプルコードでは、緑色の円を 0°、60°の位置で切り取って扇形を作っています。