Reducing 3D Interpolation size based off radial components

Question

I have the following data which represnt points in a 3D grid:

x = array([ 0, 0.08885313,  0.05077321,  0.05077321,  0.03807991,  0.03807991,
        0.03807991,  0.02538661,  0.02538661,  0.0126933 ,  0.0126933 ,
       -0.        , -0.        , -0.0126933 , -0.0126933 , -0.02538661,
       -0.02538661, -0.02538661, -0.03807991, -0.05077321, -0.05077321,
       -0.05077321, -0.06346652, -0.07615982, -0.11423973, -0.12693304,
       -0.13962634])
y = array([ 0, -0.15231964, -0.08885313, -0.17770625, -0.08885313, -0.10154643,
       -0.12693304, -0.07615982, -0.08885313, -0.07615982, -0.08885313,
       -0.07615982, -0.10154643, -0.08885313, -0.10154643, -0.07615982,
       -0.08885313, -0.17770625, -0.10154643, -0.08885313, -0.11423973,
       -0.12693304, -0.10154643, -0.08885313, -0.17770625, -0.12693304,
       -0.11423973])
z = array([ 0, 1.21839241, 0.78673339, 1.21839241, 0.70318648, 0.82850684,
       0.96078945, 0.64748854, 0.71014872, 0.63356406, 0.71711096,
       0.65445078, 0.77977115, 0.73103545, 0.77977115, 0.68926199,
       0.73799769, 1.19750569, 0.85635581, 0.84243133, 0.96078945,
       1.02344963, 0.93294048, 0.97471393, 1.24624138, 1.20446793,
       1.22535466])

I am attempting to fit a surface to these points. I want to fit the surface such that it only extends between the data points and the origin.

I have tried the following code:

data = np.array([x, y, z]).T 
  
# Define mathematical function for curve fitting 
def func(xy, a, b, c, d, e, f): 
    x, y = xy 
    return a + b*x + c*y + d*x**2 + e*y**2 + f*x*y 
  
# Perform curve fitting 
popt, pcov = curve_fit(func, (x, y), z) 
  
# Print optimized parameters 
print(popt) 
  
# Create 3D plot of the data points and the fitted curve 
fig = plt.figure() 
ax = fig.add_subplot(111, projection='3d') 
ax.scatter(x, y, z, color='blue') 
x_range = np.linspace(-0.2, 0.2, 50) 
y_range = np.linspace(-0.2, 0.2, 50) 
X, Y = np.meshgrid(x_range, y_range) 
Z = func((X, Y), *popt) 
ax.plot_surface(X, Y, Z, color='red', alpha=0.5) 

ax.set_xlim([-0.2, 0.2])  
ax.set_ylim([-0.2, 0.2])
ax.set_zlim([0,1.2])

plt.show()

This code gives me the following plot. I have drawn in blue lines where I would like the interpolation to be between. How would I add two different angles that would limit my plot in the radial direction so that the interpolation is between these lines?

Edit

Using the helpful answer provided I now am getting the following plot. However as you can see beyond the limit of y being -0.3 the interpolated surface begins to go down again. Is there a way I can apply a constraint such that the interpolation after my data only increases.

Answer 1

The limiting lines seem to lie to the left and to the right from the zero-point (0, 0, 0), across the x-axis. This means that two points with the minimum and maximum values of x could be used to define both of those lines. You could find their coefficients with numpy.linalg.lstsq, like in this answer here. The XY-projections of the lines can be used to set a limiting condition (with numpy.where) on the surface being plotted:

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
from numpy import ones, vstack
from numpy.linalg import lstsq

x = np.array([ 0, 0.08885313,  0.05077321,  0.05077321,  0.03807991,  0.03807991,
        0.03807991,  0.02538661,  0.02538661,  0.0126933 ,  0.0126933 ,
       -0.        , -0.        , -0.0126933 , -0.0126933 , -0.02538661,
       -0.02538661, -0.02538661, -0.03807991, -0.05077321, -0.05077321,
       -0.05077321, -0.06346652, -0.07615982, -0.11423973, -0.12693304,
       -0.13962634])
y = np.array([ 0, -0.15231964, -0.08885313, -0.17770625, -0.08885313, -0.10154643,
       -0.12693304, -0.07615982, -0.08885313, -0.07615982, -0.08885313,
       -0.07615982, -0.10154643, -0.08885313, -0.10154643, -0.07615982,
       -0.08885313, -0.17770625, -0.10154643, -0.08885313, -0.11423973,
       -0.12693304, -0.10154643, -0.08885313, -0.17770625, -0.12693304,
       -0.11423973])
z = np.array([ 0, 1.21839241, 0.78673339, 1.21839241, 0.70318648, 0.82850684,
       0.96078945, 0.64748854, 0.71014872, 0.63356406, 0.71711096,
       0.65445078, 0.77977115, 0.73103545, 0.77977115, 0.68926199,
       0.73799769, 1.19750569, 0.85635581, 0.84243133, 0.96078945,
       1.02344963, 0.93294048, 0.97471393, 1.24624138, 1.20446793,
       1.22535466])

# The data set seems to be within two boundaries located to the left and to the right from the (0, 0, 0) point across the x-axis
# Find the maximum and minimum x values
i_max, = np.where(np.isclose(x, np.max(x)))
i_min, = np.where(np.isclose(x, np.min(x)))
x_max = x[i_max][0]
x_min = x[i_min][0]

# Find the corresponding y values (i.e. y values of those points where x is maximum)
y_max = y[i_max][0]
y_min = y[i_min][0]

# Function finding 2D-line coefficients, based on this answer https://stackoverflow.com/a/21566184/3715182
def line_coeffs(points):
    x_coords, y_coords = zip(*points)
    A = vstack([x_coords, ones(len(x_coords))]).T
    # y = a*x + b
    a, b = lstsq(A, y_coords, rcond=None)[0]
    return (a, b)

# Find coefficients of the two lines "limiting" all points left and right across the x-axis in the XY-plane
k1_max, k2_max = line_coeffs([(x_max, y_max), (0, 0)])
k1_min, k2_min = line_coeffs([(x_min, y_min), (0, 0)])
  
# Define mathematical function for curve fitting 
def func(xy, a, b, c, d, e, f):
    x, y = xy
    return a + b*x + c*y + d*x**2 + e*y**2 + f*x*y
  
# Perform curve fitting 
popt, pcov = curve_fit(func, (x, y), z)
  
# Print optimized parameters 
print(popt)
  
# Create 3D plot of the data points and the fitted curve 
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z, color='blue')
x_range = np.linspace(-0.2, 0.2, 1000)
y_range = np.linspace(-0.2, 0.2, 1000)
X, Y = np.meshgrid(x_range, y_range)
Z = func((X, Y), *popt)
# Limit the surface with a condition, forcing its XY-projections to be within the area limited by two lines
ax.plot_surface(np.where(X >= (Y - k2_min)/k1_min, np.where(X <= (Y - k2_max)/k1_max, X, np.nan), np.nan), Y, Z, color='red', alpha=0.5)

ax.set_xlim([-0.2, 0.2])
ax.set_ylim([-0.2, 0.2])
ax.set_zlim([0,1.2])

ax.set_xlabel('x')  
ax.set_ylabel('y')
ax.set_zlabel('z')

plt.show()

The result: