3D curved surfaces

2021-09-03

2025-12-30

3D, matplotlib, python

자주 있는 일은 아니지만 3차원 곡면을 그릴 때가 있습니다.
어떤 분은 원자를 표현하느라, 또는 쇠구슬을 표현하느라 구가 필요할지도 모릅니다.
저는 업무상 태양이 하늘에 떠 있는 지점을 고민할 때가 많아서 반구가 필요합니다.
과거에는 원자의 3차원 에너지를 표현하느라 이런 그림이 필요했습니다.

1. 데이터 준비

wikipedia: Spherical coordinate system

wikipedia: 구면 좌표계

구나 반구를 그리려면 구면 좌표계로 정의된 데이터가 필요합니다.
간단한 삼각함수를 사용해서 구면 데이터를 만듭니다.
방위각(azimuthal angle, $\theta$)과 극고도각(polar angle, $\varphi$)을 나열한 뒤 여기에 아래 수식을 적용합니다.

$$
\begin{aligned}
x &= r \cos{\varphi} \sin{\theta}\\
y &= r \sin{\varphi} \sin{\theta}\\
z &= r \cos{\theta}
\end{aligned}$$

코드로는 이렇게 정리됩니다.

# angles
polars = np.linspace(0, 180, 19)
azimuths = np.linspace(0, 360, 37)

# points
from itertools import product

df = pd.DataFrame(product(polars, azimuths), columns=["azi", "polar"])
df["x"] = df.apply(lambda x: np.cos(np.deg2rad(x[1]))*np.sin(np.deg2rad(x[0])), axis=1)
df["y"] = df.apply(lambda x: np.sin(np.deg2rad(x[1]))*np.sin(np.deg2rad(x[0])), axis=1)
df["z"] = df.apply(lambda x: np.cos(np.deg2rad(x[0])), axis=1)

df

실행 결과

2. 3D 곡면 도형 그리기

matplotlib: 3D surface (colormap)
matplotlib: mpl_toolkits.mplot3d.axes3d.Axes3D

주어진 데이터로 만들어지는 3D 표면 형상을 만드는데는 .plot_surface()가 적격입니다.
앞서 만든 데이터로 구(sphere)와 반구(hemisphere)를 만들어봅시다.

2.1. sphere

우리가 만든 좌표가 이미 구의 좌표입니다.
방위각 $\varphi$를 0도에서 360도, 극고도각 $\theta$를 0도에서 180도까지 변화시켰고, 이는 3차원 공간의 모든 방향에 해당되기 때문입니다.
그렇다면 남은 일은, DataFrame에 1차원으로 들어있는 각각의 좌표를 차원에 맞게 2차원으로 바꿔주는 것 뿐입니다.
극고도각과 방위각의 경우의 수가 각각 19가지, 37가지이므로 .reshape((19, 37)을 적용합니다.
matplotlib은 구면좌표계가 아니라 직교좌표계를 사용합니다. x, y, z로 변환합니다.

fig, ax = plt.subplots(figsize=(5, 5), 
                       subplot_kw={"projection":"3d"})
ax.plot_surface(df["x"].values.reshape((19, 37)), 
                df["y"].values.reshape((19, 37)), 
                df["z"].values.reshape((19, 37)), 
                alpha=0.3)

구가 그려졌습니다.
그런데 z축 방향으로 조금 찌그러져서 납작한 느낌이 나네요.
aspect ratio를 1:1:1로 맞춰줍시다.
matplotlib 3차원 Axes(Axes3D)에서 aspect ratio를 지정하는 명령은 .set_box_aspect()입니다.

fig, ax = plt.subplots(figsize=(5, 5), 
                       subplot_kw={"projection":"3d"})
ax.plot_surface(df["x"].values.reshape((19, 37)), 
                df["y"].values.reshape((19, 37)), 
                df["z"].values.reshape((19, 37)), 
                alpha=0.3)
ax.set_box_aspect((1, 1, 1))

2.2. partial sphere

극좌표계를 사용해서 구의 일부를 쪼개볼 수도 있습니다.
plot_surface()안에 넣는 x, y, z에다 slicing을 추가하면 됩니다.
극고도각은 전체 범위를 사용하고 방위각은 마지막 9개를 사용하지 않겠다고 하면 이런 그림이 그려집니다.

fig, ax = plt.subplots(figsize=(5, 5), 
                       subplot_kw={"projection":"3d"})
ax.plot_surface(df["x"].values.reshape((19, 37))[:,:-9], 
                df["y"].values.reshape((19, 37))[:,:-9], 
                df["z"].values.reshape((19, 37))[:,:-9], 
                alpha=0.3)
ax.set_box_aspect((1, 1, 1))

2.3. 선 색깔 조정

우리는 지금 면을 그리고 있습니다.
따라서 edgecolor (ec), linestyle (ls), linewidth (lw)가 적용됩니다.
불필요한 axis도 지워버립시다. 명령은 똑같이 ax.axis(False)입니다.

fig, ax = plt.subplots(figsize=(5, 5), 
                       constrained_layout=True,
                       subplot_kw={"projection":"3d"})
ax.plot_surface(df["x"].values.reshape((19, 37))[:,:-9], 
                df["y"].values.reshape((19, 37))[:,:-9], 
                df["z"].values.reshape((19, 37))[:,:-9], 
                ec="k", lw=0.2, ls=":",
                alpha=0.3)
ax.set_box_aspect((1, 1, 1))
ax.axis(False)

아까보다 한결 깔끔해 보입니다.

2.4. 면 색깔 조정

면 색깔을 지정하는 매개변수는 color입니다.
colors="g"처럼 색상을 의미하는 약어나 숫자를 넣으면 적용됩니다.

fig, ax = plt.subplots(figsize=(5, 5), 
                       constrained_layout=True,
                       subplot_kw={"projection":"3d"})
ax.plot_surface(df["x"].values.reshape((19, 37))[:,:-9], 
                df["y"].values.reshape((19, 37))[:,:-9], 
                df["z"].values.reshape((19, 37))[:,:-9], 
                ec="k", lw=0.2, ls=":",
                color="g",
                alpha=0.3)
ax.set_box_aspect((1, 1, 1))
ax.axis(False)

조명 방향도 바꿀 수 있습니다.
matplotlib.colors의 LightSource 클래스를 사용해 조명 방향을 지정합니다.
두 개의 인자로 방위각과 극고도각을 도(degree) 단위로 입력합니다.

from matplotlib.colors import LightSource

fig, axs = plt.subplots(ncols=3,
                        figsize=(15, 5), 
                        constrained_layout=True,
                        subplot_kw={"projection":"3d"})

lightsources = [LightSource(30, 30), LightSource(-30, 60), LightSource(0, 0)]
titles = ["azi: 30, polar: 30", "azi: -30, polar: 60", "azi: 0, polar: 0"]
for ax, ls, title in zip(axs, lightsources, titles):
    ax.plot_surface(df["x"].values.reshape((19, 37))[:,:-9], 
                    df["y"].values.reshape((19, 37))[:,:-9], 
                    df["z"].values.reshape((19, 37))[:,:-9], 
                    ec="none", lw=0, ls="-",
                    color="w", lightsource=ls,
                    alpha=1)
    ax.set_title(title, fontsize="large")
    ax.set_box_aspect((1, 1, 1))
    ax.axis(False)

멋진 그라데이션을 입히고 싶다면 cmap 매개변수를 사용할 수 있습니다.
z축 방향 값에 따른 그라데이션이 매겨집니다.
Matplotlib이 제공하는 gradation 이름을 집어넣습니다.
다만 cmap과 lightsource는 함께 쓰일 수 없으니 유의해야 합니다.

fig, ax = plt.subplots(figsize=(5, 5), 
                       constrained_layout=True,
                       subplot_kw={"projection":"3d"})
ax.plot_surface(df["x"].values.reshape((19, 37))[:,:-9], 
                df["y"].values.reshape((19, 37))[:,:-9], 
                df["z"].values.reshape((19, 37))[:,:-9], 
                ec="k", lw=0.2, ls=":",
                cmap="inferno",
                alpha=0.3)
ax.set_box_aspect((1, 1, 1))
ax.axis(False)

3. 응용

3.1. 구와 반구

stackoverflow: Axes3D 여백 없애기

구와 반구를 동시에 그려봅니다.
구는 위에서 만든 코드를 그대로 사용하면 되고,
반구는 극고도각을 절반만 사용하면 됩니다.

fig, axs = plt.subplots(ncols=2, 
                        figsize=(10, 5), 
                        constrained_layout=True, 
                        subplot_kw={"projection":"3d"})

# sphere
axs[0].plot_surface(df["x"].values.reshape((19, 37)), 
                    df["y"].values.reshape((19, 37)), 
                    df["z"].values.reshape((19, 37)), 
                    cmap="viridis", 
                    alpha=0.3)
axs[0].set_box_aspect((1, 1, 1))
axs[0].set_title("sphere", pad=12)

# hemisphere
axs[1].plot_surface(df["x"].iloc[:370].values.reshape((10, 37)), 
                    df["y"].iloc[:370].values.reshape((10, 37)), 
                    df["z"].iloc[:370].values.reshape((10, 37)), 
                    cmap="viridis", 
                    alpha=0.3)

axs[1].set_box_aspect((1, 1, 0.5))
axs[1].set_title("hemisphere", pad=12)

for ax in axs:
    ax.set_xlabel("X")
    ax.set_ylabel("Y")
    ax.set_zlabel("Z")

X, Y, Z축 범위가 -1~1을 약간 벗어나 있습니다.
matplotlib 기본 설정 문제입니다.
이를 해결하려면 위 코드에서 fig, axs를 생성한 뒤에 짧은 코드를 추가합니다.

from mpl_toolkits.mplot3d.axis3d import Axis
if not hasattr(Axis, "_get_coord_info_old"):
    def _get_coord_info_new(self, renderer):
        mins, maxs, centers, deltas, tc, highs = self._get_coord_info_old(renderer)
        mins += deltas / 4
        maxs -= deltas / 4
        return mins, maxs, centers, deltas, tc, highs
    Axis._get_coord_info_old = Axis._get_coord_info  
    Axis._get_coord_info = _get_coord_info_new

동서남북, 그리고 상하를 잇는 수직선을 그립니다.
좌표 원점이라고 봐도 됩니다.
여기서 한 가지 주의할 점이 있습니다.
matplotlib이 3D 도형과 2D 도형의 위치관계를 깔끔하기 처리하지 못합니다.
3D 도형에 가려질 부분과 가려지지 않을 부분을 따로 그리고 표현해줘야 합니다.
동서남북을 나타내는 N, E, W, S도 적당한 위치에 넣어줍니다.

fig, axs = plt.subplots(ncols=2, 
                        figsize=(10, 5), 
                        subplot_kw={"projection":"3d"})

# Axes3D 딱 붙이기
# https://stackoverflow.com/questions/16488182/removing-axes-margins-in-3d-plot
from mpl_toolkits.mplot3d.axis3d import Axis
if not hasattr(Axis, "_get_coord_info_old"):
    def _get_coord_info_new(self, renderer):
        mins, maxs, centers, deltas, tc, highs = self._get_coord_info_old(renderer)
        mins += deltas / 4
        maxs -= deltas / 4
        return mins, maxs, centers, deltas, tc, highs
    Axis._get_coord_info_old = Axis._get_coord_info  
    Axis._get_coord_info = _get_coord_info_new

# sphere
axs[0].plot_surface(df["x"].values.reshape((19, 37)), 
                    df["y"].values.reshape((19, 37)), 
                    df["z"].values.reshape((19, 37)), 
                    cmap="viridis", ec="k", lw=0.1, 
                    alpha=0.3)
axs[0].set_box_aspect((1, 1, 1))
axs[0].set_title("sphere", pad=12)

# hemisphere
axs[1].plot_surface(df["x"].iloc[:370].values.reshape((10, 37)), 
                    df["y"].iloc[:370].values.reshape((10, 37)), 
                    df["z"].iloc[:370].values.reshape((10, 37)), 
                    cmap="viridis", ec="k", lw=0.1, 
                    alpha=0.3)

axs[1].set_box_aspect((1, 1, 0.5))
axs[1].set_title("hemisphere", pad=12)

for ax in axs:
    ax.set_xlabel("X")
    ax.set_ylabel("Y")
    ax.set_zlabel("Z")
    ax.view_init(azim=235)
    # make margins
    ax.margins(0)
    # center axis
    zmin = -1 if ax == axs[0] else 0
    ax.plot([-1.1, 1.1], [0, 0], [0, 0], c="k", zorder=-10)
    ax.plot([0, 0, 0], [-1.1, 1.1, 0], c="k", zorder=-10)
    ax.plot([0, 0], [0, 0], [zmin, 1], c="k", zorder=-10)
    ax.plot([0, 0], [0, 0], [1, 1.1], c="k", zorder=3)
    
    if ax == axs[0]:
        ax.plot([-1, -1.1], [0, 0], [0, 0], c="k", zorder=3)
        ax.plot([0, 0], [-1, -1.1], [0, 0], c="k", zorder=3)
     
    # 동서남북 표시
    font_text = {"fontweight":"bold", "fontsize":15, "ha":"center", "va":"center"}
    ax.text(-1.4, 0, 0, "W", transform=ax.transData, fontdict=font_text)
    ax.text(1.4, 0, 0, "E", transform=ax.transData, fontdict=font_text)
    ax.text(0, 1.4, 0, "N", transform=ax.transData, fontdict=font_text, zorder=-10)
    ax.text(0, -1.4, 0, "S", transform=ax.transData, fontdict=font_text)
    ax.axis(False)

3.2. magnetocrystalline anisotropy

wikipedia: magnetocrystalline anisotropy

결정자기이방성(magnetocrystalline anisotropy)이라는 말을 들어보신 분은 거의 없으실 겁니다.
자성(magnetism)을 공부할 때 나오는 용어인데, 방위에 따라 다른 에너지라고 대충 넘어가셔도 좋습니다.
중요한 것은 이번엔 구나 반구가 아니라 울퉁불퉁한 모양을 만들 것이라는 점입니다.
원리는 간단합니다.
이 글의 처음에서 구면좌표계는 다음과 같은 식으로 표현된다고 했습니다.

$$
\begin{aligned}
x &= r \cos{\varphi} \sin{\theta}\\
y &= r \sin{\varphi} \sin{\theta}\\
z &= r \cos{\theta}
\end{aligned}$$

앞에서는 $r = 1$로 두었지만, 여기서는 $r = f(\varphi, \theta)$입니다.
방위각과 극고도각을 이용해서 물질마다 다르게 r을 정의합니다.
모양이 울퉁불퉁한만큼 아까 10도 단위로 쪼갠 공간을 이번엔 5도 단위로 쪼갭니다.
자세한 설명은 생략하고 결과 그림과 코드만 간략히 보여드리겠습니다.
5도 단위로 공간 분할

# angles
polars = np.linspace(0, 179, 37)
azimuths = np.linspace(0, 360, 73)

# points
from itertools import product

df = pd.DataFrame(product(polars, azimuths), columns=["azi", "polar"])
df["x"] = df.apply(lambda x: np.cos(np.deg2rad(x[1]))*np.sin(np.deg2rad(x[0])), axis=1)
df["y"] = df.apply(lambda x: np.sin(np.deg2rad(x[1]))*np.sin(np.deg2rad(x[0])), axis=1)
df["z"] = df.apply(lambda x: np.cos(np.deg2rad(x[0])), axis=1)

방향에 따른 $r$ 계산: 자기이방성 에너지

# 자기이방성 계수 K1, K2
K1_Co = 45
K2_Co = 15

K1_Fe = 4.8
K2_Fe = 0.5

K1_Ni = -0.5 
K2_Ni = -0.2

# 방향에 따른 자기이방성 에너지 계산.
# 그림을 그리는 것이 목적이므로 보기 좋도록 적절히 스케일링
def calc_uni(K1, K2, df):
    return K1*(df["x"]**2 + df["y"]**2)

def calc_cubic(K1, K2, df):
    return K1*(df["x"]**2 * df["y"]**2 + \
               df["y"]**2 * df["z"]**2 + \
               df["z"]**2 * df["x"]**2) + \
           K2*(df["x"]**2 * df["y"]**2 * df["z"]**2)

df["E_Co"] = df.apply(lambda x: calc_uni(K1_Co, K2_Co, x), axis=1) + 1
df["E_Fe"] = df.apply(lambda x: calc_cubic(K1_Fe, K2_Fe, x), axis=1) + 1
df["E_Ni"] = df.apply(lambda x: calc_cubic(K1_Ni, K2_Ni, x), axis=1)*3 +1

# 극좌표계를 직교좌표계로 변환
df["x_Co"] = df["E_Co"] * df["x"]
df["y_Co"] = df["E_Co"] * df["y"]
df["z_Co"] = df["E_Co"] * df["z"]

df["x_Fe"] = df["E_Fe"] * df["x"]
df["y_Fe"] = df["E_Fe"] * df["y"]
df["z_Fe"] = df["E_Fe"] * df["z"]

df["x_Ni"] = df["E_Ni"] * df["x"]
df["y_Ni"] = df["E_Ni"] * df["y"]
df["z_Ni"] = df["E_Ni"] * df["z"]

세 가지 물질(Co, Fe, Ni)의 자기이방성 에너지 시각화

fig, axs = plt.subplots(ncols=3, 
                        figsize=(15, 6), 
                       constrained_layout=True,
                        subplot_kw={"projection":"3d"})
for ax, mat in zip(axs, ["Co", "Fe", "Ni"]):
    ax.plot_surface(df[f"x_{mat}"].values.reshape((37, 73))[:,:-19], 
                    df[f"y_{mat}"].values.reshape((37, 73))[:,:-19], 
                    df[f"z_{mat}"].values.reshape((37, 73))[:,:-19], 
                    ec="k", lw=0.2,
                    color="w", lightsource=LightSource(0, 10),
                    alpha=1)
    ax.plot(df[f"x_{mat}"].values.reshape((37, 73))[:,-20], 
            df[f"y_{mat}"].values.reshape((37, 73))[:,-20], 
            df[f"z_{mat}"].values.reshape((37, 73))[:,-20], 
            c="k", alpha=1)
    ax.plot(df[f"x_{mat}"].values.reshape((37, 73))[:,0], 
            df[f"y_{mat}"].values.reshape((37, 73))[:,0], 
            df[f"z_{mat}"].values.reshape((37, 73))[:,0], 
            c="k", alpha=1)
    ax.set_box_aspect((1, 1, 1))
    ax.axis(False)
    ax.set_title(mat, fontsize="xx-large")
    
fig.suptitle("crystalline anisotropy energy surface\n", 
             color="blue", fontsize="xx-large")

제 박사학위 주제와 밀접한 그림입니다.
교과서 그림을 모방해 학위논문(2011)에도 같은 그림을 그려서 넣었습니다.
mathematica를 사용해서 그린 그림입니다.

도움이 되셨나요? 카페인을 투입하시면 다음 포스팅으로 변환됩니다

PythonData Science

1. 데이터 준비

2. 3D 곡면 도형 그리기

2.1. sphere

2.2. partial sphere

2.3. 선 색깔 조정

2.4. 면 색깔 조정

3. 응용

3.1. 구와 반구

3.2. magnetocrystalline anisotropy