插值后极坐标图中的缺失点

Question

我有 17 个测量点，我想在它们之间推断一个圆圈。这是我尝试过的：

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

# Measurement data
x1 = np.array([[-144.00, -101.80, -101.80, -75.00, -53.00, -53.00, 0.00, 0.00, 0.00, 0.00, 0.00, 53.00, 53.00, 75.00, 101.80, 101.80, 144.00]])
y1 = np.array([[0.00, 101.80, -101.80, 0.00, 53.00, -53.00, 144.00, 75.00, 0.00, -75.00, -144.00, 53.00, -53.00, 0.00, 101.80, -101.80, 0.00]])
z1 = np.array([148.3861807, 148.9051447, 148.0415147, 147.9976293, 147.98485, 147.9579673, 148.89261, 148.0217707, 147.9312247, 147.7952, 147.3225247, 148.2489567, 148.2120013, 148.3169953, 149.578092, 147.9356893, 148.556672])

# Convert to polar coords
r1 = np.sqrt(x1**2 + y1**2)
t1 = np.arctan2(y1, x1)

# Add cyclic points
for i in range(0,x1.shape[1]):
    if np.arctan2(y1[0,i], x1[0,i]) == math.pi:
        print(np.sqrt(x1[0,i]**2 + y1[0,i]**2), np.arctan2(y1[0,i], x1[0,i]))
        r1 = np.append(r1,([np.sqrt(x1[0,i]**2 + y1[0,i]**2)]))
        t1 = np.append(t1,([-np.arctan2(y1[0,i], x1[0,i])]))
        z1 = np.append(z1,([z1[i]]))

# New points
r2, t2 = np.meshgrid(np.linspace(np.min(r1), np.max(r1), num=50),np.linspace(np.min(t1), np.max(t1), num=50))

# Griddata function used to interpolate between scattered data
z2 = griddata(np.concatenate((np.array([t1]), np.array([r1])), axis=0).T, np.array([z1]).T, (t2, r2), method='linear')

# Surface plots
fig, ax = plt.subplots(figsize=(10,10), subplot_kw=dict(projection='polar'))

ax.contourf(t2, r2, np.squeeze(z2), 50, cmap='jet')

for i in range(len(t1)):
    plt.text(t1[i], r1[i], "{:.2f}".format(z1[i]), ha="center", va="center", color="k")

plt.show()

然而，生成的极坐标图在中间缺少点。我怎样才能摆脱它们？

Resultant polar plot

Answer 1

这是我的解决方法：

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

# Measurement data
x = np.array([[-144.00, -101.80, -101.80, -75.00, -53.00, -53.00, 0.00, 0.00, 0.00, 0.00, 0.00, 53.00, 53.00, 75.00, 101.80, 101.80, 144.00]])
y = np.array([[0.00, 101.80, -101.80, 0.00, 53.00, -53.00, 144.00, 75.00, 0.00, -75.00, -144.00, 53.00, -53.00, 0.00, 101.80, -101.80, 0.00]])
z1 = np.array([[148.3861807, 148.9051447, 148.0415147, 147.9976293, 147.98485, 147.9579673, 148.89261, 148.0217707, 147.9312247, 147.7952, 147.3225247, 148.2489567, 148.2120013, 148.3169953, 149.578092, 147.9356893, 148.556672]])

# Convert x & y to polar coordinates
rho   = np.sqrt(x**2 + y**2)
theta = np.arctan2(y, x)

# Round up rho values
rho = np.round(rho, 1)

# Ensure that theta axis repeats itself
tmp = z1
for i in range(0,x.shape[1]):
    if np.arctan2(y[0,i], x[0,i]) == math.pi:
        #print(np.sqrt(x[0,i]**2 + y[0,i]**2), np.arctan2(y[0,i], x[0,i]))
        rho   = np.append(rho,([np.sqrt(x[0,i]**2 + y[0,i]**2)]))
        theta = np.append(theta,([-np.arctan2(y[0,i], x[0,i])]))
        tmp    = np.append(tmp,([z1[0,i]]))
z1 = tmp
del tmp

# Create polar grid
Rho, Theta = np.meshgrid(np.unique(rho),np.unique(theta))
#plt.polar(Theta, Rho, 'ro')
#plt.show()

# Create 2D z1
Z1 = np.empty([len(np.unique(Theta)), len(np.unique(Rho))])
Z1[:] = np.NaN

for i in range(0,len(theta)):
    for m in range(0,Theta.shape[0]):
        for n in range(0,Theta.shape[1]):
            if theta[i] == Theta[m,n] and rho[i] == Rho[m,n]:
                Z1[m,n] = z1[i]
            if rho[i] == Rho[m,n] and Rho[m,n] == 0: # take care of (0,0)
                Z1[m,n] = z1[i]

#fig, ax = plt.subplots(subplot_kw=dict(projection='polar'))
#ax.contourf(Theta, Rho, Z1, 50, cmap='jet')
#plt.show()

# RectBivariateSpline requires an nxn grid
rho_interp = np.linspace(np.min(Rho[0,:]), np.max(Rho[0,:]), num=np.max(Rho.shape))
Z1_interp = np.empty([np.max(Rho.shape), np.max(Rho.shape)])

for i in range(0,Z1.shape[0]):
    f = interpolate.interp1d(Rho[i,:], Z1[i,:], kind='linear')
    Z1_interp[i,:] = f(rho_interp)
    
Rho_interp, Theta_interp = np.meshgrid(rho_interp,np.unique(Theta))

#fig, ax = plt.subplots(subplot_kw=dict(projection='polar'))
#ax.contourf(Theta_interp, Rho_interp, Z1_interp, 50, cmap='jet')
#plt.show()

# Interpolate and extrapolate onto fine grid
interp_spline = RectBivariateSpline(rho_interp, np.unique(Theta), Z1_interp, kx=1, ky=1) # kx, ky : degrees of the bivariate spline

rho_fine   = np.linspace(0, 150, num=500)
theta_fine = np.linspace(-math.pi, math.pi, num=500)
Rho_fine, Theta_fine = np.meshgrid(rho_fine, theta_fine)
Z1_fine = interp_spline(rho_fine, theta_fine)

# Generate the plot
fig, ax = plt.subplots(figsize=(10,10), subplot_kw=dict(projection='polar'))

plt.contourf(Theta_fine, Rho_fine, Z1_fine, 150, cmap='jet')

plt.colorbar()

plt.show()

脚本的某些部分看起来很荒谬，但 RectBivariateSpline 似乎可以胜任。限制是它需要一个 n×n 网格。

如果您找到解决此问题的更优雅的解决方案，请告诉我。

Final plot