使用matplotlib繪製循環數據

發表於 | 分類於

循環數據（cyclic data）的定義詳見此處：Cyclic data - Oxford Reference

引言

我在完成一個關於分析片上衍射單元的衍射場的項目時，曾需要將衍射場的光學模擬數據可視化。數據儲存於一個csv文件中，以下是數據的樣例：


數據樣例

每一列的含義如下：

x y E_real E_imag E_arg
點在平面坐標系的x值 點在平面坐標系的y值 實部 虛部 輻角（相位）

其中$-10^{-5}\leqslant x\leqslant 10^{-6},\ -5\times 10^{-6} \leqslant y\leqslant 5\times 10^{-6}$。

很顯然，數據儲存了光波在衍射場上某一點對應的的複數形式。其中光的相位是一種典型的循環數據，且由以下公式得出：

使用contourf()函數可能會遇到的問題

對於此類需求，很容易會想到使用matplotlib.pyplot中的contourf()函數繪製等高線圖。代碼如下：

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
51
52
53
54
55
56
57
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import ticker

df = pd.read_csv('/dataset/dataset.csv')
x_unique = np.sort(df['x'].unique())
y_unique = np.sort(df['y'].unique())

# 將 x 和 y 轉換為網格坐標
X, Y = np.meshgrid(x_unique, y_unique)

# 將 E_arg 轉換為網格形式
Z = np.zeros_like(X)
for i in range(len(x_unique)):
    for j in range(len(y_unique)):
        Z[j, i] = df[(df['x'] == x_unique[i]) & (df['y'] == y_unique[j])]['E_arg'].values[0]

# 繪製等高線圖
fig, ax = plt.subplots()
levels = np.linspace(-np.pi, np.pi, 129)
cs = ax.contourf(X, Y, Z, levels, cmap=plt.get_cmap('twilight'))
ax.set_title('E_arg')
ax.set_xlabel('x')
ax.set_ylabel('y')

# 設置顏色條
cbar = fig.colorbar(cs,fraction=0.1, pad=0.15, shrink=0.9, anchor=(0.0, 0.3))
tick_locator = ticker.MaxNLocator(nbins=6)  # colorbar上的刻度值個數
cbar.locator = tick_locator
cbar.ax.tick_params()
cbar.update_ticks()

# x軸設置
x_min = -1e-5
x_max = 1e-6
x_step = 2e-6
x_ticks = np.arange(x_min, x_max, x_step)

# y軸設置
y_min = -5e-6
y_max = 5e-6
y_step = 2e-6
y_ticks = np.arange(y_min, y_max, y_step)

# 設置x軸和y軸的範圍和刻度
plt.xlim(x_min, x_max)
plt.xticks(x_ticks)
plt.ylim(y_min, y_max)
plt.yticks(y_ticks)

#設置圖像像素及大小
plt.rcParams['figure.figsize']=(6.0, 4.0)
plt.rcParams['savefig.dpi'] = 600 #圖片像素
plt.rcParams['figure.dpi'] = 600 #解析度

plt.show()

因為數據是循環數據，因此我選擇了cyclic colormaps中的twilight樣式作為圖像的colormap。得到的圖像如下：


E_arg的等高線圖

仔細觀察會發現，在相位小於且接近$\pi$和相位大於且接近$-\pi$的點之間，出現的並不是平滑的顏色過渡，而是形成了一條包含所有顏色的極細的色帶。這與循環數據的特性不符，因為光波的相位週期為$2\pi$，$\pi$和$-\pi$應是等同的。

這個錯誤很可能是因為contourf()的插值方法沒有認識到$\pi$和$-\pi$的等價性，而在兩點之間進行了錯誤的插值，這在放大後的圖像中能夠清晰地觀察到：


E_arg的等高線圖中$\pi$與$-\pi$的交界處

解決方法

1. 使用scatter()函數

要解決上述問題，一個可行的替代方案是改為使用matplotlib.pyplot中的scatter()函數繪製散點圖。這裡需要注意的是，由於我們的需求仍然是類似等高線圖那樣用不同顏色填充滿整個區域的圖像，因此繪製散點圖的第一個前提是原始數據足夠大（在本文所舉的例子中，數據集有數十萬行），以達到能將散點圖在實際輸出中趨同於等高線圖的效果。這種方法的代碼如下：

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
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import ticker

df = pd.read_csv('/dataset/dataset.csv')
x_unique = np.sort(df['x'].unique())
y_unique = np.sort(df['y'].unique())

# 提取 x ， y 和 E_arg 的值
x = df['x'].values
y = df['y'].values
E_arg = df['E_arg'].values

# 繪製等高線圖
fig, ax = plt.subplots()
scatter = ax.scatter(x, y, c=E_arg, cmap='twilight')

# 設置顏色條
cbar = fig.colorbar(cs,fraction=0.1, pad=0.15, shrink=0.9, anchor=(0.0, 0.3))
tick_locator = ticker.MaxNLocator(nbins=6)  # colorbar上的刻度值個數
cbar.locator = tick_locator
cbar.ax.tick_params()
cbar.update_ticks()

# x軸設置
x_min = -1e-5
x_max = 1e-6
x_step = 2e-6
x_ticks = np.arange(x_min, x_max, x_step)

# y軸設置
y_min = -5e-6
y_max = 5e-6
y_step = 2e-6
y_ticks = np.arange(y_min, y_max, y_step)

# 設置x軸和y軸的範圍和刻度
plt.xlim(x_min, x_max)
plt.xticks(x_ticks)
plt.ylim(y_min, y_max)
plt.yticks(y_ticks)

#設置圖像像素及大小
plt.rcParams['figure.figsize']=(6.0, 4.0)
plt.rcParams['savefig.dpi'] = 600 #圖片像素
plt.rcParams['figure.dpi'] = 600 #解析度

plt.show()

這種方法得到的圖像如下：


E_arg的散點圖

可以看到，錯誤的色帶確實消失了。這是因為色帶上的數據是插值得到的，並非來源於原始數據，而散點圖只會繪製所有的原始數據，並不會做插值計算。

然而，這種方法也存在一些問題。除去上文中已經提到的原始數據需要足夠大這一要求之外，原始數據也需要足夠規律。仔細觀察得到的散點圖，會發現有些區域有錯位的情況。將圖像中錯位的區域放大觀察，能夠明顯看到不規律的數據，而錯位也是由不規律的數據所造成的視覺效果：


E_arg的散點圖中錯位的區域

實際上，本例中的原始數據相對來說已經足夠規律了。如果是更為散亂的數據，可能會得到更差的結果。

2. 使用三角函數映射

如果仍要堅持使用contourf()函數繪製等高線圖，也可以考慮對相位使用三角函數映射。以$\sin x$為例，$\sin\pi=\sin(-\pi)=0$，這樣一來，相位小於且接近$\pi$和相位大於且接近$-\pi$的點之間就是連續的平穩過渡了。代碼如下（其實只是在第17行加上了np.sin()）：

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
51
52
53
54
55
56
57
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import ticker

df = pd.read_csv('/dataset/dataset.csv')
x_unique = np.sort(df['x'].unique())
y_unique = np.sort(df['y'].unique())

# 將 x 和 y 轉換為網格坐標
X, Y = np.meshgrid(x_unique, y_unique)

# 將 E_arg 轉換為網格形式
Z = np.zeros_like(X)
for i in range(len(x_unique)):
    for j in range(len(y_unique)):
        Z[j, i] = np.sin(df[(df['x'] == x_unique[i]) & (df['y'] == y_unique[j])]['E_arg'].values[0])

# 繪製等高線圖
fig, ax = plt.subplots()
levels = np.linspace(-np.pi, np.pi, 129)
cs = ax.contourf(X, Y, Z, levels, cmap=plt.get_cmap('twilight'))
ax.set_title('E_arg')
ax.set_xlabel('x')
ax.set_ylabel('y')

# 設置顏色條
cbar = fig.colorbar(cs,fraction=0.1, pad=0.15, shrink=0.9, anchor=(0.0, 0.3))
tick_locator = ticker.MaxNLocator(nbins=6)  # colorbar上的刻度值個數
cbar.locator = tick_locator
cbar.ax.tick_params()
cbar.update_ticks()

# x軸設置
x_min = -1e-5
x_max = 1e-6
x_step = 2e-6
x_ticks = np.arange(x_min, x_max, x_step)

# y軸設置
y_min = -5e-6
y_max = 5e-6
y_step = 2e-6
y_ticks = np.arange(y_min, y_max, y_step)

# 設置x軸和y軸的範圍和刻度
plt.xlim(x_min, x_max)
plt.xticks(x_ticks)
plt.ylim(y_min, y_max)
plt.yticks(y_ticks)

#設置圖像像素及大小
plt.rcParams['figure.figsize']=(6.0, 4.0)
plt.rcParams['savefig.dpi'] = 600 #圖片像素
plt.rcParams['figure.dpi'] = 600 #解析度

plt.show()

但是，這種辦法有一個嚴重的問題：它勢必會將兩個不相等的相位映射到同一個值上。這時，此種替代方法的可行性就主要取決於本身的需求是否允許了。在本文的例子中，它最終生成了如下圖像：


sin(E_arg)的等高線圖

3. 其他的可能方案

問題的核心在於使相位$\pi$和相位$-\pi$連續，因此或許可以先使用numpy中的unwrap()函數預處理數據使之連續，再使用contourf()函數繪製等高線圖。但是這種方法會使數據不再僅限於$(-\pi,\ \pi)$區間，此時可能需要自訂colormap，將若干段長度為2的twilight拼接，以得到能夠包含所有新數據的colormap。另外，由於unwrap()函數只能用於一維數組，而數據的坐標是二維的，因此如何使用unwrap()函數也是一個難點。

最後，在StackOverflow上也有一個類似的問題可供參考：python - Handling cyclic data with matplotlib contour/contourf - Stack Overflow。

0%