# Return Z-value Of Xy Coordinate

## 01 May 2019 - 1 answer

I have a set of xy cooridnates that generate a contour. For the code below, these cooridnates are from groups `A` and `B` in the `df`. I have also created a separate xy cooridnate that is called from `C1_X` and `C1_Y`. However this isn't used in generating the contour itself. It is a separate xy coordinate.

Question: Is it possible to return the z-value of the contour at the `C1_X``C1_Y` coordinate?

I have found a separate question that is similar: multivariate spline interpolation in python scipy?. The figure in that question displays what I'm hoping to return but I just want the z-value for one xy coordinate.

The `contour` in this question is normalised so values fall between `-1` and `1`. I'm hoping to return the z-value for `C1_X` and `C1_Y`, which is the white scatter point seen in the figure beneath the code.

I have attempted to return the z-value for this point using:

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

But I'm returning an error: ```raise TypeError('x must be strictly increasing') TypeError: x must be strictly increasing```

I have commented out this function so the code runs.

Side note: This code is written for an animation.

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

I am hoping to return the `z` value for the white scatter point. Intended Output will display the normalised `z` value `(-1,1)` for `C1_X`,`C1_Y`.

Upon visual inspection this would be between`0.6` and `0.8`

# Edit 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):

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

cov = np.array([
])

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)

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 = {
'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 = {
'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)
elif max(containing_levels) < 0:
lev = min(containing_levels)

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:

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

If you have an arbitrary cloud of (X, Y, Z) points and you want to interpolate the z-coordinate of some (x, y) point, you have a number of different options. The simplest is probably to just use `scipy.interpolate.interp2d` to get the z-value:

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

Since the grid you have appears to be regular, you may be better off using `scipy.interpolate.RectBivariateSpline`, which has a very similar interface, but is specifically made for regular grids:

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

Since you have a regular meshgrid, you can also do

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

Notice that the dimensions are flipped between the plotting arrays and the interpolation arrays. Plotting treats axis 0 as rows, i.e. Y, while the interpolation functions treat axis 0 as X. Rather than transposing, you could also switch the X and Y inputs, leaving Z intact for a similar end result, e.g.:

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

Alternatively, you could change all your plotting code to swap the inputs as well, but that would be too much work at this point. Whatever you do, pick an approach and stick with it. As long as you do it consistently, they should all work.

If you use one of the `scipy` approaches (recommended), keep the object `f` around to interpolate any further points you might want.

If you want a more manual approach, you can do something like find the three closest (X, Y, Z) points to (x, y), and find the value of the plane between them at (x, y). For example:

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

Since your data is on a regular grid, it might be easier to do something like bilinear interpolation on the quad surrounding (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))
``````