Построение контурной карты интерполированной функции: несоответствие результатов для разных участков данных

Question

Я определил значение z = f (x, y) для сетки точек на плоскости:

import numpy as np

n_theta, n_phi = 20, 30
thetas = np.linspace(0, 1, n_theta) * np.pi
phis = np.linspace(0, 1, n_phi) * (2 * np.pi)

res = f(thetas, phis)
# 'res' is an np.array with res.shape=(600,3), where each row is [x, y, z]

Учитывая meshgrid из x 's и y, я интерполирую функцию, используя мои предыдущие данные для z, и строю контурную карту результата:

from scipy.interpolate import griddata
import matplotlib.pyplot as plt

X, Y = np.meshgrid(res[:,0], res[:,1])
Z = griddata(res[:,0:2], res[:,2], (X, Y), method='cubic')    # Interpolate function.

fig = plt.figure()
contour = plt.contourf(X, Y, Z, levels=100)

Это отлично работает, и результат (пример данных, связанных ниже) что-то вроде этого (что на самом деле не та форма, которую я искал в моем конкретном моделировании):

Однако, когда я запускаю тот же скрипт со значениями phi (горизонтальная ось) только для половины исходного интервала (или, что то же самое, взять только эту половину моих данных), я получаю контур, подобный этому:

(Ради полноты, если я сделаю это для второй половины интервала и вручную поставлю обе цифры рядом, результат будет похож на следующий, который больше похож на то, что я ожидать получения полных данных):

(При необходимости Вот примерные данные для создания графиков , за которыми следует фрагмент, который можно использовать):

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

# Full axis:

res = np.load("stackoverflow-example-data")
X, Y = np.meshgrid(res[:,0], res[:,1])
Z = griddata(res[:,0:2], res[:,2], (X, Y), method='cubic')    # Interpolate function.
fig = plt.figure()
contour = plt.contourf(X, Y, Z, levels=100)

# Half axis:

half = np.array([res[i,:] for i in range(len(res)) if res[i,0] <= np.pi])
X, Y = np.meshgrid(res[:,0], res[:,1])
Z = griddata(res[:,0:2], res[:,2], (X, Y), method='cubic')    # Interpolate function.
fig = plt.figure()
contour = plt.contourf(X, Y, Z, levels=100)

Так чего мне здесь не хватает?

1. это правильная процедура для интерполяции и построения контура так, как я хочу?
2. Разве я не должен ожидать получения аналогичных результатов для всей оси и только для половины ее? Если нет, то почему?

Answer 1

Ваши данные уже находятся в сетке. Создание новой сетки на основе этой сетки создает что-то очень грязное.

Теперь массив res содержит 600 значений x, y, z. Использование reshape(20, 30) говорит numpy, что эти 600 записей на самом деле представляют собой 20 строк по 30 столбцов. Самая «чистая» форма для отображения данных - это график scatter, показывающий только данные.

При imshow данные могут отображаться в виде изображения. Использование interpolation='nearest' "взорвало бы" каждое значение z, чтобы заполнить прямоугольный angular пиксель. Использование interpolation='bicubic' будет плавно размывать эти пиксели.

contourf создает контуры с равным значением z. В зависимости от данных наличие нескольких уровней (например, 100) поможет получить плавное изображение, но это не будет полезным по сравнению с простым отображением сглаженного изображения. Ограниченное количество уровней (например, 20) может помочь показать общую базовую форму.

Вот код, который можно сравнить и поэкспериментировать с различными подходами:

import numpy as np
import matplotlib.pyplot as plt

# Full axis:

res = np.load("stackoverflow-example-data.npy")
X = res[:, 0].reshape(20, 30)
Y = res[:, 1].reshape(20, 30)
Z = res[:, 2].reshape(20, 30)
fig, axes = plt.subplots(ncols=3, figsize=(12, 4))
axes[0].scatter(X, Y, c=Z, s=10)
axes[1].imshow(Z, origin='lower', interpolation='bicubic', aspect='auto',
               extent=[X.min(), X.max(), Y.min(), Y.max()])
contour = axes[2].contourf(X, Y, Z, levels=20)
axes[0].set_title('scatter plot')
axes[1].set_title('imshow, bicubic interpolation')
axes[2].set_title('contour plot, 20 levels')

plt.show()

Экспериментируя с различными цветными картами , можно подчеркнуть различные свойства отображаемой функции. Вы также можете преобразовать значение Z (например, np.log(Z) или np.exp(Z)), чтобы разные регионы получали более или менее детальную информацию.

for ax, cmap in zip(axes, ('viridis', 'inferno_r', 'Purples')):
    ax.imshow(Z, origin='lower', interpolation='bicubic', aspect='auto', cmap=cmap,
               extent=[X.min(), X.max(), Y.min(), Y.max()])
    ax.set_title(cmap)