Python SciPy RectSphereBivariateSpline интерполяция, генерирующая неверные данные?

Question

У меня есть данные 3D-измерения на сфере, которая очень грубая и что я хочу интерполировать. С большой помощью от @ M4rtini и @HYRY здесь, в stackoverflow, теперь я смог создать рабочий код (на основе исходного примера из примера RectSphereBivariateSpline из SciPy).

Тестовых данных можно найти здесь: testdata

""" read csv input file, post process and plot 3D data """ 
import csv 
import numpy as np 
from mayavi import mlab 
from scipy.interpolate import RectSphereBivariateSpline 

# user input 
nElevationPoints = 17 # needs to correspond with csv file 
nAzimuthPoints = 40 # needs to correspond with csv file 
threshold = - 40 # needs to correspond with how measurement data was captured 
turnTableStepSize = 72 # needs to correspond with measurement settings 
resolution = 0.125 # needs to correspond with measurement settings 

# read data from file 
patternData = np.empty([nElevationPoints, nAzimuthPoints]) # empty buffer 
ifile = open('ttest.csv') # need the 'b' suffix to prevent blank rows being inserted 
reader = csv.reader(ifile,delimiter=',') 
reader.next() # skip first line in csv file as this is only text 
for nElevation in range (0,nElevationPoints): 
    # azimuth 
    for nAzimuth in range(0,nAzimuthPoints): 
     patternData[nElevation,nAzimuth] = reader.next()[2] 
ifile.close() 

# post process 
def r(thetaIndex,phiIndex): 
    """r(thetaIndex,phiIndex): function in 3D plotting to return positive vector length from patternData[theta,phi]""" 
    radius = -threshold + patternData[thetaIndex,phiIndex] 
    return radius 

#phi,theta = np.mgrid[0:nAzimuthPoints,0:nElevationPoints] 
theta = np.arange(0,nElevationPoints) 
phi = np.arange(0,nAzimuthPoints) 
thetaMesh, phiMesh = np.meshgrid(theta,phi) 
stepSizeRad = turnTableStepSize * resolution * np.pi/180 
theta = theta * stepSizeRad 
phi = phi * stepSizeRad 

# create new grid to interpolate on 
phiIndex = np.arange(1,361) 
phiNew = phiIndex*np.pi/180 
thetaIndex = np.arange(1,181) 
thetaNew = thetaIndex*np.pi/180 
thetaNew,phiNew = np.meshgrid(thetaNew,phiNew) 
# create interpolator object and interpolate 
data = r(thetaMesh,phiMesh) 
theta[0] += 1e-6 # zero values for theta cause program to halt; phi makes no sense at theta=0 
lut = RectSphereBivariateSpline(theta,phi,data.T) 
data_interp = lut.ev(thetaNew.ravel(),phiNew.ravel()).reshape((360,180)).T 

def rInterp(theta,phi): 
    """rInterp(theta,phi): function in 3D plotting to return positive vector length from interpolated patternData[theta,phi]""" 
    thetaIndex = theta/(np.pi/180) 
    thetaIndex = thetaIndex.astype(int) 
    phiIndex = phi/(np.pi/180) 
    phiIndex = phiIndex.astype(int) 
    radius = data_interp[thetaIndex,phiIndex] 
    return radius 
# recreate mesh minus one, needed otherwise the below gives index error, but why?? 
phiIndex = np.arange(0,360) 
phiNew = phiIndex*np.pi/180 
thetaIndex = np.arange(0,180) 
thetaNew = thetaIndex*np.pi/180 
thetaNew,phiNew = np.meshgrid(thetaNew,phiNew) 

x = (rInterp(thetaNew,phiNew)*np.cos(phiNew)*np.sin(thetaNew)) 
y = (-rInterp(thetaNew,phiNew)*np.sin(phiNew)*np.sin(thetaNew)) 
z = (rInterp(thetaNew,phiNew)*np.cos(thetaNew)) 

# plot 3D data 
obj = mlab.mesh(x, y, z, colormap='jet') 
obj.enable_contours = True 
obj.contour.filled_contours = True 
obj.contour.number_of_contours = 20 
mlab.show() 

Хотя код работает, полученный график значительно отличается от не-интерполированных данных см изображений

, как ссылка.

Кроме того, при запуске интерактивного сеанса data_interp намного больше по значению (> 3e5), чем исходные данные (это около 20 макс.).

Кто-нибудь знает, что я могу делать неправильно?

Answer 1

Я, кажется, решил это!

Для этого я попытался экстраполировать, тогда как я мог только интерполировать эти разбросанные данные. SO новая сетка интерполяции должна идти только до тета = 140 градусов или около того.

Но самое важное изменение - это добавление параметра s = 900 в вызов RectSphereBivariateSpline.

Так теперь у меня есть следующий код:

""" read csv input file, post process and plot 3D data """ 
import csv 
import numpy as np 
from mayavi import mlab 
from scipy.interpolate import RectSphereBivariateSpline 

# user input 
nElevationPoints = 17 # needs to correspond with csv file 
nAzimuthPoints = 40 # needs to correspond with csv file 
threshold = - 40 # needs to correspond with how measurement data was captured 
turnTableStepSize = 72 # needs to correspond with measurement settings 
resolution = 0.125 # needs to correspond with measurement settings 

# read data from file 
patternData = np.empty([nElevationPoints, nAzimuthPoints]) # empty buffer 
ifile = open('ttest.csv') # need the 'b' suffix to prevent blank rows being inserted 
reader = csv.reader(ifile,delimiter=',') 
reader.next() # skip first line in csv file as this is only text 
for nElevation in range (0,nElevationPoints): 
    # azimuth 
    for nAzimuth in range(0,nAzimuthPoints): 
     patternData[nElevation,nAzimuth] = reader.next()[2] 
ifile.close() 

# post process 
def r(thetaIndex,phiIndex): 
    """r(thetaIndex,phiIndex): function in 3D plotting to return positive vector length from patternData[theta,phi]""" 
    radius = -threshold + patternData[thetaIndex,phiIndex] 
    return radius 

#phi,theta = np.mgrid[0:nAzimuthPoints,0:nElevationPoints] 
theta = np.arange(0,nElevationPoints) 
phi = np.arange(0,nAzimuthPoints) 
thetaMesh, phiMesh = np.meshgrid(theta,phi) 
stepSizeRad = turnTableStepSize * resolution * np.pi/180 
theta = theta * stepSizeRad 
phi = phi * stepSizeRad 

# create new grid to interpolate on 
phiIndex = np.arange(1,361) 
phiNew = phiIndex*np.pi/180 
thetaIndex = np.arange(1,141) 
thetaNew = thetaIndex*np.pi/180 
thetaNew,phiNew = np.meshgrid(thetaNew,phiNew) 
# create interpolator object and interpolate 
data = r(thetaMesh,phiMesh) 
theta[0] += 1e-6 # zero values for theta cause program to halt; phi makes no sense at theta=0 
lut = RectSphereBivariateSpline(theta,phi,data.T,s=900) 
data_interp = lut.ev(thetaNew.ravel(),phiNew.ravel()).reshape((360,140)).T 

def rInterp(theta,phi): 
    """rInterp(theta,phi): function in 3D plotting to return positive vector length from interpolated patternData[theta,phi]""" 
    thetaIndex = theta/(np.pi/180) 
    thetaIndex = thetaIndex.astype(int) 
    phiIndex = phi/(np.pi/180) 
    phiIndex = phiIndex.astype(int) 
    radius = data_interp[thetaIndex,phiIndex] 
    return radius 
# recreate mesh minus one, needed otherwise the below gives index error, but why?? 
phiIndex = np.arange(0,360) 
phiNew = phiIndex*np.pi/180 
thetaIndex = np.arange(0,140) 
thetaNew = thetaIndex*np.pi/180 
thetaNew,phiNew = np.meshgrid(thetaNew,phiNew) 

x = (rInterp(thetaNew,phiNew)*np.cos(phiNew)*np.sin(thetaNew)) 
y = (-rInterp(thetaNew,phiNew)*np.sin(phiNew)*np.sin(thetaNew)) 
z = (rInterp(thetaNew,phiNew)*np.cos(thetaNew)) 

# plot 3D data 
intensity = rInterp(thetaNew,phiNew) 
obj = mlab.mesh(x, y, z, scalars = intensity, colormap='jet') 
obj.enable_contours = True 
obj.contour.filled_contours = True 
obj.contour.number_of_contours = 20 
mlab.show() 

Полученный график сравнивающий хорошо к оригинальному без интерполяции данных:

Я не в полной мере понять, почему ей следует устанавливается в 900, так как в документации RectSphereBivariateSpline указано, что s = 0 для интерполяции. Однако при чтении документации немного подробнее я получу некоторое представление:

Выбор оптимального значения s может быть деликатной задачей. Рекомендуемые значения s зависят от точности значений данных. Если пользователь имеет представление о статистических ошибках данных, она также может найти правильную оценку для s. Предполагая, что если она задает правильные s, интерполятор будет использовать сплайн f (u, v), который точно воспроизводит функцию, лежащую в основе данных, она может оценить сумму ((r (i, j) -s (u (i), v (j))) ** 2), чтобы найти хорошую оценку для этого s. Например, если она знает, что статистические ошибки на ее значениях r (i, j) не превышают 0,1, она может ожидать, что значение хорошего s должно иметь значение не больше, чем u.size * v.size * (0,1) ** 2. Если ничего не известно о статистической ошибке в r (i, j), s должно определяться методом проб и ошибок. Лучше всего начинать с очень большого значения s (определить полином наименьших квадратов и соответствующую верхнюю границу fp0 для s), а затем постепенно уменьшать значение s (скажем, в 10 раз в начале, т. Е. s = fp0/10, fp0/100, ... и более тщательно, поскольку приближение показывает более подробно), чтобы получить более близкие посадки.