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

Question

我正在尝试绘制具有分级浓度值的 windrose。遵循 this post 的建议和一些修改，我创建了一个情节。但是，在 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()

结果是:

作为散点图，效果很好，看起来像

我不确定如何以简洁的方式提供数据，但我们将不胜感激!!

Answer 1

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

顺便说一句，请注意 'jet' 是一个看起来不错的颜色图，但它是 very misleading因为它会在错误的位置创建黄色高光。 (此外，将 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':使用反转的 'magma' 颜色图，因此高值变暗，低值变浅(参见 the docs 关于可能更适合的其他颜色图说明数据，但尽量不要使用'jet')

hist2d 的返回值是直方图值矩阵、bin 边界(x 和 y)和色 block 集合(可用作颜色图的输入)。

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()