python - 如何摆脱windrose pcolormap plot python上的不连续性

Question

我正在尝试绘制具有分级浓度值的风玫瑰图。根据这篇文章的建议和一些修改，我创建了一个情节。但是，在 0 度左右存在不连续性。任何帮助将不胜感激！

这是我的代码：

wd = list(merge_all_apr['Wind Dir (10s deg)'])
conc = list(merge_all_apr['Mean_CO2'])
ws = list(merge_all_apr['Wind Spd (km/h)'])

wd_rad = np.radians(np.array(wd))
conc = np.array(conc, dtype=np.float)

wind_speed = np.linspace(min(ws), max(ws), 16)

WD, WS = np.meshgrid(np.linspace(0, 2*np.pi, 36), wind_speed)

print (WS)

Z = interpolate.griddata((wd_rad, ws), oz, (WD, WS), method='linear')

fig, ax = plt.subplots(subplot_kw={"projection": "polar"})
cmap = plt.get_cmap('jet')
#cmap.set_under('none')
img = ax.pcolormesh(WD, WS, Z, cmap=cmap, vmin=0, vmax=40, alpha = 0.70)
ax.set_theta_zero_location('N')
ax.set_theta_direction(-1) 
plt.colorbar(img)
plt.show()

结果是：

作为一个分散，它工作正常，看起来像

我不确定如何以简洁的方式提供数据，但任何帮助将不胜感激！

标签： pythonnumpymatplotlib

Answer 1

您还可以提供wd_rad小于0和大于2*piusing的值，小值np.where相加，大值2*pi相减2*pi。np.tile(ws, 2)然后np.tile(conc, 2)将 的扩展版本wd_rad与相同的浓度值相关联。也使用这些扩展值interpolate.griddata确保浓度值环绕在0和2*pi。

顺便说一句，请注意“jet”是一个看起来不错的颜色图，但由于它会在错误的位置创建黄色高光，因此非常具有误导性。（此外，将 pandas 列转换为列表非常慢且占用内存，最好将它们保留为 numpy 数组格式。）

下面的代码假设oz问题中的数组与conc.

import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate

# wd = merge_all_apr['Wind Dir (10s deg)']
# conc = merge_all_apr['Mean_CO2']
# ws = merge_all_apr['Wind Spd (km/h)']

N = 100
wd = np.random.uniform(0, 360, N)
conc = np.random.uniform(0, 40, N)
ws = np.random.uniform(0, 45, N)

wd_rad = np.radians(np.array(wd))
conc = np.array(conc, dtype=np.float)

wd_rad_ext = np.where(wd_rad < np.pi, wd_rad + 2 * np.pi, wd_rad - 2 * np.pi)

wind_speed = np.linspace(min(ws), max(ws), 16)

WD, WS = np.meshgrid(np.linspace(0, 2 * np.pi, 37), wind_speed)

Z = interpolate.griddata((np.hstack([wd_rad, wd_rad_ext]), np.tile(ws, 2)), np.tile(conc, 2),
                         (WD, WS), method='linear', rescale=True)

fig, axes = plt.subplots(ncols=2, figsize=(10, 4), subplot_kw={"projection": "polar"})
cmap = plt.get_cmap('magma')
for ax in axes:
    if ax == axes[0]:
        img = ax.pcolormesh(WD, WS, Z, cmap=cmap, vmin=0, vmax=40)
    else:
        img = ax.scatter(wd_rad, ws, c=conc, cmap=cmap, vmin=0, vmax=40)
    ax.set_theta_zero_location('N')
    ax.set_theta_direction(-1)
    plt.colorbar(img, ax=ax, pad=0.12)
plt.show()

如果您不想要插值，但想要绘制 10 度宽的线段来表示每个区域，plt.hist2d则可以使用。

特殊参数hist2d：

• bins=(np.linspace(0, 2 * np.pi, 37), np.linspace(min(ws), max(ws), 17))：将风向划分为36个区域（37个边界）；速度将分为16个区域
• weights=conc：使用浓度而不是通常的直方图来计算每个小区域的值的数量；在同一小区域测量多个浓度时，取平均值
• cmin=0.001：当区域的浓度值小于0.001
• cmap = 'magma_r'：使用反转的“岩浆”颜色图，因此高值会变暗，而低值会变浅（请参阅有关可能更适合说明数据的其他颜色图的文档，但尽量不要使用“喷射”）

的返回值hist2d是直方图值的矩阵、bin 边界（x 和 y）和彩色补丁的集合（可用作颜色图的输入）。

import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate

N = 100
wd = np.random.uniform(0, 360, N)
conc = np.random.uniform(0, 40, N)
ws = np.random.uniform(0, 45, N)

wd_rad = np.radians(np.array(wd))
conc = np.array(conc, dtype=np.float)

fig, axes = plt.subplots(ncols=2, figsize=(10, 4), subplot_kw={"projection": "polar"})
cmap = 'magma_r'
for ax in axes:
    if ax == axes[0]:
        _, _, _, img = ax.hist2d(wd_rad, ws, bins=(np.linspace(0, 2 * np.pi, 37), np.linspace(min(ws), max(ws), 17)),
                                 weights=conc, cmin=0.001, cmap=cmap, vmin=0, vmax=40)
    else:
        img = ax.scatter(wd_rad, ws, c=conc, cmap=cmap, vmin=0, vmax=40)
    ax.set_theta_zero_location('N')
    ax.set_theta_direction(-1)
    plt.colorbar(img, ax=ax, pad=0.12)
plt.show()