Вернуть z-значение координаты xy - PullRequest
1 голос
/ 01 мая 2019

У меня есть набор координат xy, которые генерируют контур.Для приведенного ниже кода эти координаты относятся к группам A и B в df.Я также создал отдельный xy cooridnate, который вызывается из C1_X и C1_Y.Однако это не используется при создании самого контура.Это отдельная координата xy.

Вопрос: Можно ли вернуть значение z контура по координате C1_X C1_Y?

Iнашли отдельный вопрос, который похож: многомерная сплайн-интерполяция в Python Scipy? .Цифра в этом вопросе отображает то, что я надеюсь вернуть, но мне просто нужно z-значение для одной координаты xy.

contour в этом вопросе нормализовано, поэтому значения находятся в пределах от -1 до 1.Я надеюсь вернуть z-значение для C1_X и C1_Y, которое является белой точкой рассеяния, показанной на рисунке под кодом.

Я попытался вернуть z-значение для этой точки, используя:

# Attempt at returning the z-value for C1 
f = RectBivariateSpline(X, Y, normPDF)
z = f(d['C1_X'], d['C1_Y']) 
print(z)

Но я возвращаю ошибку: raise TypeError('x must be strictly increasing') TypeError: x must be strictly increasing

Я прокомментировалэта функция, поэтому код работает.

Примечание: этот код написан для анимации.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.stats as sts
import matplotlib.animation as animation
from mpl_toolkits.axes_grid1 import make_axes_locatable
from scipy.interpolate import RectBivariateSpline

DATA_LIMITS = [0, 15]

def datalimits(*data):
return DATA_LIMITS 

def mvpdf(x, y, xlim, ylim, radius=1, velocity=0, scale=0, theta=0):
    X,Y = np.meshgrid(np.linspace(*xlim), np.linspace(*ylim))
    XY = np.stack([X, Y], 2)
    PDF = sts.multivariate_normal([x, y]).pdf(XY)
    return X, Y, PDF

def mvpdfs(xs, ys, xlim, ylim, radius=None, velocity=None, scale=None, theta=None):
    PDFs = []
    for i,(x,y) in enumerate(zip(xs,ys)):
        X, Y, PDF = mvpdf(x, y, xlim, ylim)
        PDFs.append(PDF)
    return X, Y, np.sum(PDFs, axis=0)

fig, ax = plt.subplots(figsize = (10,6))
ax.set_xlim(DATA_LIMITS)
ax.set_ylim(DATA_LIMITS)

line_a, = ax.plot([], [], 'o', c='red', alpha = 0.5, markersize=5,zorder=3)
line_b, = ax.plot([], [], 'o', c='blue', alpha = 0.5, markersize=5,zorder=3)
scat = ax.scatter([], [], s=5**2,marker='o', c='white', alpha = 1,zorder=3)

lines=[line_a,line_b] 
scats=[scat] 

cfs = None

def plotmvs(tdf, xlim=datalimits(df['X']), ylim=datalimits(df['Y']), fig=fig, ax=ax):    
    global cfs  
    if cfs:
        for tp in cfs.collections:
            tp.remove()
    df = tdf[1]
    PDFs = []

    for (group, gdf), group_line in zip(df.groupby('group'), (line_a, line_b)):
        group_line.set_data(*gdf[['X','Y']].values.T)
        X, Y, PDF = mvpdfs(gdf['X'].values, gdf['Y'].values, xlim, ylim)
        PDFs.append(PDF)

    for (group, gdf), group_line in zip(df.groupby('group'), lines+scats):
            if group in ['A','B']:
                group_line.set_data(*gdf[['X','Y']].values.T)
                kwargs = {
                'xlim': xlim,
                'ylim': ylim
                }
                X, Y, PDF = mvpdfs(gdf['X'].values, gdf['Y'].values, **kwargs)
                PDFs.append(PDF)

            #plot white scatter point from C1_X, C1_Y
            elif group in ['C']:
                gdf['X'].values, gdf['Y'].values
                scat.set_offsets(gdf[['X','Y']].values)

    # normalize PDF by shifting and scaling, so that the smallest value is -1 and the largest is 1
    normPDF = (PDFs[0]-PDFs[1])/max(PDFs[0].max(),PDFs[1].max())


    ''' Attempt at returning z-value for C1_X, C1_Y '''
    ''' This is the function that I am trying to write that will '''
    ''' return the contour value '''

    #f = RectBivariateSpline(X[::-1, :], Y[::-1, :], normPDF[::-1, :]) 
    #z = f(d['C1_X'], d['C1_Y']) 
    #print(z)

    cfs = ax.contourf(X, Y, normPDF, cmap='jet', alpha = 1, levels=np.linspace(-1,1,10),zorder=1)

    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", size="5%", pad=0.1)
    cbar = fig.colorbar(cfs, ax=ax, cax=cax)
    cbar.set_ticks([-1,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1])

    return  cfs.collections + [scat] + [line_a,line_b] 

''' Sample Dataframe '''

n = 1
time = range(n)  

d = ({
    'A1_X' :    [3],
    'A1_Y' :    [6],
    'A2_X' :    [6],
    'A2_Y' :    [10],
    'B1_X' :    [12],
    'B1_Y' :    [2],
    'B2_X' :    [14],
    'B2_Y' :    [4],
    'C1_X' :    [4],
    'C1_Y' :    [6],                     
    })

# a list of tuples of the form ((time, group_id, point_id, value_label), value)
tuples = [((t, k.split('_')[0][0], int(k.split('_')[0][1:]), k.split('_')[1]), v[i])
    for k,v in d.items() for i,t in enumerate(time) ]

df = pd.Series(dict(tuples)).unstack(-1)
df.index.names = ['time', 'group', 'id']

#Code will eventually operate with multiple frames
interval_ms = 1000
delay_ms = 2000
ani = animation.FuncAnimation(fig, plotmvs, frames=df.groupby('time'), interval=interval_ms, repeat_delay=delay_ms,)

plt.show()

enter image description here

Я надеюсь вернуть zзначение для белой точки рассеяния.В предполагаемом выходе отобразится нормализованное значение z (-1,1) для C1_X, C1_Y.

При визуальном осмотре это может быть между 0.6 и 0.8

Редактировать 2:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.stats as sts
import matplotlib.animation as animation
from mpl_toolkits.axes_grid1 import make_axes_locatable
from scipy.interpolate import RectBivariateSpline
import matplotlib.transforms as transforms

DATA_LIMITS = [-85, 85]

def datalimits(*data):
    return DATA_LIMITS  # dmin - spad, dmax + spad

def rot(theta):
    theta = np.deg2rad(theta)
    return np.array([
        [np.cos(theta), -np.sin(theta)],
        [np.sin(theta), np.cos(theta)]
    ])

def getcov(radius=1, scale=1, theta=0):
    cov = np.array([
        [radius*(scale + 1), 0],
        [0, radius/(scale + 1)]
    ])

    r = rot(theta)
    return r @ cov @ r.T

def mvpdf(x, y, xlim, ylim, radius=1, velocity=0, scale=0, theta=0):

    X,Y = np.meshgrid(np.linspace(*xlim), np.linspace(*ylim))
    XY = np.stack([X, Y], 2)
    x,y = rot(theta) @ (velocity/2, 0) + (x, y)
    cov = getcov(radius=radius, scale=scale, theta=theta)

    PDF = sts.multivariate_normal([x, y], cov).pdf(XY)

    return X, Y, PDF

def mvpdfs(xs, ys, xlim, ylim, radius=None, velocity=None, scale=None, theta=None):
    PDFs = []
    for i,(x,y) in enumerate(zip(xs,ys)):
        kwargs = {
            'radius': radius[i] if radius is not None else 0.5,
            'velocity': velocity[i] if velocity is not None else 0,
            'scale': scale[i] if scale is not None else 0,
            'theta': theta[i] if theta is not None else 0,
            'xlim': xlim,
            'ylim': ylim
        }
        X, Y, PDF = mvpdf(x, y,**kwargs)
        PDFs.append(PDF)

    return X, Y, np.sum(PDFs, axis=0)

fig, ax = plt.subplots(figsize = (10,6))

ax.set_xlim(DATA_LIMITS)
ax.set_ylim(DATA_LIMITS)

line_a, = ax.plot([], [], 'o', c='red', alpha = 0.5, markersize=3,zorder=3)
line_b, = ax.plot([], [], 'o', c='blue', alpha = 0.5, markersize=3,zorder=3)
lines=[line_a,line_b] ## this is iterable!

offset = lambda p: transforms.ScaledTranslation(p/82.,0, plt.gcf().dpi_scale_trans)
trans = plt.gca().transData

scat = ax.scatter([], [], s=5,marker='o', c='white', alpha = 1,zorder=3,transform=trans+offset(+2) )
scats=[scat] 

cfs = None

def plotmvs(tdf, xlim=None, ylim=None, fig=fig, ax=ax):
    global cfs  
    if cfs:
        for tp in cfs.collections:
            tp.remove()

    df = tdf[1]

    if xlim is None: xlim = datalimits(df['X'])
    if ylim is None: ylim = datalimits(df['Y'])

    PDFs = []

    for (group, gdf), group_line in zip(df.groupby('group'), lines+scats):
        if group in ['A','B']:
            group_line.set_data(*gdf[['X','Y']].values.T)
            kwargs = {
            'radius': gdf['Radius'].values if 'Radius' in gdf else None,
            'velocity': gdf['Velocity'].values if 'Velocity' in gdf else None,
            'scale': gdf['Scaling'].values if 'Scaling' in gdf else None,
            'theta': gdf['Rotation'].values if 'Rotation' in gdf else None,
            'xlim': xlim,
            'ylim': ylim
            }
            X, Y, PDF = mvpdfs(gdf['X'].values, gdf['Y'].values, **kwargs)
            PDFs.append(PDF)
        elif group in ['C']:
            gdf['X'].values, gdf['Y'].values
            scat.set_offsets(gdf[['X','Y']].values)

    normPDF = (PDFs[0]-PDFs[1])/max(PDFs[0].max(),PDFs[1].max())


    def get_contour_value_of_point(point_x, point_y, X, Y, Z, precision=10000):

        CS = ax.contour(X, Y, Z, 100)
        containing_levels = []
        for cc, lev in zip(CS.collections, CS.levels):
            for pp in cc.get_paths():
                if pp.contains_point((point_x, point_y)):
                    containing_levels.append(lev)

        if max(containing_levels) == 0:
            return 0
        else:
            if max(containing_levels) > 0:
                lev = max(containing_levels)
                adj = 1. / precision
            elif max(containing_levels) < 0:
                lev = min(containing_levels)
                adj = -1. / precision

            is_inside = True
            while is_inside:
                CS = ax.contour(X, Y, Z, [lev])
                for pp in CS.collections[0].get_paths():
                    if not pp.contains_point((point_x, point_y)):
                       is_inside = False
                if is_inside:
                    lev += adj

            return lev - adj

    print(get_contour_value_of_point(d['C1_X'], d['C1_Y'], X, Y, normPDF))


    cfs = ax.contourf(X, Y, normPDF, cmap='viridis', alpha = 1, levels=np.linspace(-1,1,10),zorder=1)

    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", size="5%", pad=0.1)
    cbar = fig.colorbar(cfs, ax=ax, cax=cax)
    cbar.set_ticks([-1,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1])

    return  cfs.collections + [scat] + [line_a,line_b] 

''' Sample Dataframe '''

n = 10
time = range(n)  

d = ({
    'A1_X' :    [3],
    'A1_Y' :    [6],
    'A2_X' :    [6],
    'A2_Y' :    [10],
    'B1_X' :    [12],
    'B1_Y' :    [2],
    'B2_X' :    [14],
    'B2_Y' :    [4],
    'C1_X' :    [4],
    'C1_Y' :    [6],                     
    })

# a list of tuples of the form ((time, group_id, point_id, value_label), value)
tuples = [((t, k.split('_')[0][0], int(k.split('_')[0][1:]), k.split('_')[1]), v[i])
    for k,v in d.items() for i,t in enumerate(time) ]

df = pd.Series(dict(tuples)).unstack(-1)
df.index.names = ['time', 'group', 'id']

#Code will eventually operate with multiple frames
interval_ms = 1000
delay_ms = 2000
ani = animation.FuncAnimation(fig, plotmvs, frames=df.groupby('time'), interval=interval_ms, repeat_delay=delay_ms,)

plt.show()

Ответы [ 2 ]

1 голос
/ 12 мая 2019

Вот неэлегантный подход грубой силы. * Предполагая, что у нас есть значения X, Y и Z, давайте определим функцию, которая рисует пользовательские контурные линии снова и снова, пока они не пересекаются с точкой с заданным пользователем уровнем точности (в ваших данных сделайте Z = normPDF).

def get_contour_value_of_point(point_x, point_y, X, Y, Z, precision=10000):
    fig, ax = plt.subplots()
    CS = ax.contour(X, Y, Z, 100)
    containing_levels = []
    for cc, lev in zip(CS.collections, CS.levels):
        for pp in cc.get_paths():
            if pp.contains_point((point_x, point_y)):
                containing_levels.append(lev)

    if max(containing_levels) == 0:
        return 0
    else:
        if max(containing_levels) > 0:
            lev = max(containing_levels)
            adj = 1. / precision
        elif max(containing_levels) < 0:
            lev = min(containing_levels)
            adj = -1. / precision

        is_inside = True
        while is_inside:
            CS = ax.contour(X, Y, Z, [lev])
            for pp in CS.collections[0].get_paths():
                if not pp.contains_point((point_x, point_y)):
                   is_inside = False
            if is_inside:
                lev += adj

        return lev - adj

Более подробно: то, что это делает, это рисование исходной контурной карты с 100 уровнями, а затем поиск списка уровней контуров, чьи многоугольники содержат точку ввопрос.Затем мы находим самый узкий уровень (самый высокий, если уровни положительные, или самый низкий, если уровни отрицательные).Оттуда мы сжимаем уровень небольшими шагами (в соответствии с желаемым уровнем точности), проверяя, находится ли точка внутри полигонов.Когда точка больше не находится в многоугольнике контура, мы знаем, что нашли правильный уровень (последний, содержащий точку).

В качестве примера, мы можем использовать контур в библиотеке Matplotlib:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

delta = 0.025
x = np.arange(-3.0, 3.0, delta)
y = np.arange(-2.0, 2.0, delta)
X, Y = np.meshgrid(x, y)
Z1 = np.exp(-X**2 - Y**2)
Z2 = np.exp(-(X - 1)**2 - (Y - 1)**2)
Z = (Z1 - Z2) * 2

При такой настройке get_contour_value_of_point(0, -0.6) возвращает 1.338399999999998, что при визуальном осмотре кажется совпадающим.get_contour_value_of_point(0, -0.6) возвращает -1.48, что также похоже на совпадение.Графики ниже для визуальной проверки.

Contour plot with point at 0, -0.6 Contour plot with point at 1, 1.5


* Я не могу гарантировать, что это охватит все виды использованияслучаев.Это покрывало те, которые я попробовал.Я бы довольно строго проверил это, прежде чем приблизиться к любой производственной среде.Я ожидал бы, что найдутся более изящные решения, чем это (например, Ответ Безумного Физика ), но это было то, что пришло мне в голову и, казалось, работало прямым, хотя и грубым, способом.

1 голос
/ 12 мая 2019

Если у вас есть произвольное облако точек (X, Y, Z) и вы хотите интерполировать координату z некоторой точки (x, y), у вас есть несколько различных вариантов.Самое простое, вероятно, просто использовать scipy.interpolate.interp2d, чтобы получить z-значение:

f = interp2d(X.T, Y.T, Z.T)
z = f(x, y)

Поскольку сетка, которая у вас есть, выглядит регулярной, вам может быть лучше использовать scipy.interpolate.RectBivariateSpline, который имеет очень похожий интерфейс, но специально создан для обычных сеток:

f = RectBivariateSpline(X.T, Y.T, Z.T)
z = f(x, y)

Поскольку у вас есть обычная сетка, вы также можете сделать

f = RectBivariateSpline(X[0, :], Y[:, 0], Z.T)
z = f(x, y)

Обратите внимание, что размеры переворачиваются между массивами построения и интерполяционными массивами.При построении диаграммы ось 0 рассматривается как строки, т. Е. Y, в то время как функции интерполяции обрабатывают ось 0 как X. Вместо транспонирования можно также переключать входы X и Y, оставляя Z без изменений для аналогичного конечного результата, например:

f = RectBivariateSpline(Y, X, Z)
z = f(y, x)

В качестве альтернативы, вы могли бы изменить весь свой код печати, чтобы также поменять местами входные данные, но на этом этапе было бы слишком много работы.Что бы вы ни делали, выберите подход и придерживайтесь его.Пока вы делаете это последовательно, все они должны работать.

Если вы используете один из scipy подходов (рекомендуется), держите объект f рядом для интерполяции любых других точек, которые вы можете захотеть.

Если вам нужен более ручной подход, вы можете сделать что-то вроде поиска трех ближайших (X, Y, Z) точек к (x, y) и найти значение плоскости между ними в (x,у).Например:

def interp_point(x, y, X, Y, Z):
    """
    x, y: scalar coordinates to interpolate at
    X, Y, Z: arrays of coordinates corresponding to function
    """
    X = X.ravel()
    Y = Y.ravel()
    Z = Z.ravel()

    # distances from x, y to all X, Y points
    dist = np.hypot(X - x, Y - y)
    # indices of the nearest points
    nearest3 = np.argpartition(dist, 2)[:3]
    # extract the coordinates
    points = np.stack((X[nearest3], Y[nearest3], Z[nearest3]))
    # compute 2 vectors in the plane
    vecs = np.diff(points, axis=0)
    # compute normal to plane
    plane = np.cross(vecs[0], vecs[1])
    # rhs of plane equation
    d = np.dot(plane, points [:, 0])
    # The final result:
    z = (d - np.dot(plane[:2], [x, y])) / plane[-1]
    return z

print(interp_point(x, y, X.T, Y.T, Z.T))

Поскольку ваши данные находятся на регулярной сетке, может быть проще выполнить что-то вроде билинейной интерполяции на окружении четырехугольника (x, y):

def interp_grid(x, y, X, Y, Z):
    """
    x, y: scalar coordinates to interpolate at
    X, Y, Z: arrays of coordinates corresponding to function
    """
    X, Y = X[:, 0], Y[0, :]

    # find matching element
    r, c = np.searchsorted(Y, y), np.searchsorted(X, x)
    if r == 0: r += 1
    if c == 0: c += 1
    # interpolate
    z = (Z[r - 1, c - 1] * (X[c] - x) * (Y[r] - y) +
         Z[r - 1, c] * (x - X[c - 1]) * (Y[r] - y) +
         Z[r, c - 1] * (X[c] - x) * (y - Y[r - 1]) +
         Z[r, c] * (x - X[c - 1]) * (y - Y[r - 1])
    ) / ((X[c] - X[c - 1]) * (Y[r] - Y[r - 1]))
    return z

print(interpolate_grid(x, y, X.T, Y.T, Z.T))
...