Python: Как интерполировать углы на сетке? - PullRequest
0 голосов
/ 26 апреля 2019

У меня есть сетка с некоторыми данными. Эти данные задаются его углом (от 0 до π). В этой сетке у меня есть еще одна меньшая сетка.

Это может выглядеть так:

enter image description here

Теперь я хочу интерполировать углы на этой сетке.

Я попробовал это с помощью scipy.interpolate.griddata, что дает хороший результат. Но есть проблема, когда углы изменяются от почти 0 до почти π (потому что середина π/2 ...)

Вот результат, и легко понять, что идет не так.

enter image description here

Как я могу справиться с этой проблемой? Спасибо! :)

Вот код для воспроизведения:

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

ax = plt.subplot()
ax.set_aspect(1)

# Simulate some given data.
x, y = np.meshgrid(np.linspace(-10, 10, 20), np.linspace(-10, 10, 20))
data = np.arctan(y / 10) % np.pi
u = np.cos(data)
v = np.sin(data)

ax.quiver(x, y, u, v, headlength=0.01, headaxislength=0, pivot='middle', units='xy')

# Create a smaller grid within.
x1, y1 = np.meshgrid(np.linspace(-1, 5, 15), np.linspace(-6, 2, 20))
# ax.plot(x1, y1, '.', color='red', markersize=2)

# Interpolate data on grid.
interpolation = griddata((x.flatten(), y.flatten()), data.flatten(), (x1.flatten(), y1.flatten()))
u1 = np.cos(interpolation)
v1 = np.sin(interpolation)
ax.quiver(x1, y1, u1, v1, headlength=0.01, headaxislength=0, pivot='middle', units='xy',
          color='red', scale=3, width=0.03)

plt.show()

Edit:

Благодаря @bubble, есть способ отрегулировать заданные углы перед интерполяцией так, чтобы результат был желаемым. Поэтому:

  1. Определить функцию выпрямления:

    def RectifyData(data):
        for j in range(len(data)):
            step = data[j] - data[j - 1]
            if abs(step) > np.pi / 2:
                data[j] += np.pi * (2 * (step < 0) - 1)
        return data
    
  2. Интерполировать следующим образом:

    interpolation = griddata((x.flatten(), y.flatten()),
                             RectifyData(data.flatten()),
                             (x1.flatten(), y1.flatten()))
    u1 = np.cos(interpolation)
    v1 = np.sin(interpolation)
    

1 Ответ

1 голос
/ 26 апреля 2019

Я пробовал прямую интерполяцию значений cos(angle) и sin(angle), но это все равно приводило к прерываниям, которые вызывают неправильные направления линий. Основная идея заключается в сокращении разрывов, например, [2.99,3.01, 0.05,0.06] должно быть преобразовано во что-то вроде этого: [2.99, 3.01, pi+0.05, pi+0.06]. Это необходимо для правильного применения алгоритма 2D-интерполяции. Почти такая же проблема возникает в следующей записи .

def get_rectified_angles(u, v):
    angles = np.arcsin(v)
    inds = u < 0
    angles[inds] *= -1
# Direct approach of removing discontinues 
#     for j in range(len(angles[1:])):  
#         if abs(angles[j] - angles[j - 1]) > np.pi / 2:
#             sel = [abs(angles[j] + np.pi - angles[j - 1]), abs(angles[j] - np.pi - angles[j-1])]
#             if np.argmin(sel) == 0:
#                 angles[j] += np.pi
#             else:
#                 angles[j] -= np.pi
    return angles


ax.quiver(x, y, u, v, headlength=0.01, headaxislength=0, pivot='middle', units='xy')

# # Create a smaller grid within.
x1, y1 = np.meshgrid(np.linspace(-1, 5, 15), np.linspace(-6, 2, 20))

angles = get_rectified_angles(u.flatten(), v.flatten())

interpolation = griddata((x.flatten(), y.flatten()), angles, (x1.flatten(), y1.flatten()))
u1 = np.cos(interpolation)
v1 = np.sin(interpolation)
ax.quiver(x1, y1, u1, v1, headlength=0.01, headaxislength=0, pivot='middle', units='xy',
          color='red', scale=3, width=0.03)

Вероятно, функция numpy.unwrap может быть использована для устранения разрывов. В случае данных 1d, numpy.interp имеет ключевое слово period для обработки периодических данных.

enter image description here

...