Dato un insieme di punti definiti nelle coordinate (X, Y, Z), interpolare il valore Z su arbitrario (X, Y)

Question

Dato un insieme di punti nelle coordinate (X, Y, Z) che sono punti su una superficie , Mi piacerebbe essere in grado di interpolare i valori Z a coordinate arbitrarie (X, Y). Ho trovato qualche successo usando mlab.griddata per interpolare i valori su una griglia, ma voglio essere in grado di chiamare una funzione di uso generale per qualsiasi coordinata (X, Y).

L'insieme di punti forma una superficie approssimativamente semisferica. Per semplificare il problema, sto provando a scrivere un metodo che interpola i valori tra i punti noti dell'emisfero definiti dalle coordinate x, yez sotto. Sebbene esista una soluzione analitica per trovare z = f (x, y) per una sfera perfetta, in modo da non dover interpolare, l'insieme effettivo di punti non sarà una sfera perfetta, quindi dovremmo assumere che abbiamo bisogno per interpolare valori a coordinate sconosciute (X, Y). Link to IPython notebook with point data

resolution = 10 
u = np.linspace(-np.pi/2, np.pi/2, resolution) 
v = np.linspace(0, np.pi, resolution) 

U, V = np.meshgrid(u, v) 

xs = np.sin(U) * np.cos(V) 
ys = np.sin(U) * np.sin(V) 
zs = np.cos(U) 

ho utilizzato scipy.interpolate.interp2d, che "restituisce una funzione la cui chiamata di metodo utilizza l'interpolazione spline per trovare il valore di nuovi punti."

def polar(xs, ys, zs, resolution=10): 
    rs = np.sqrt(np.multiply(xs, xs) + np.multiply(ys, ys)) 
    ts = np.arctan2(ys, xs) 
    func = interp2d(rs, ts, zs, kind='cubic') 
    vectorized = np.vectorize(func) 

    # Guesses 
    ri = np.linspace(0, rs.max(), resolution) 
    ti = np.linspace(0, np.pi * 2, resolution) 

    R, T = np.meshgrid(ri, ti) 
    Z = vectorized(R, T) 
    return R * np.cos(T), R * np.sin(T), Z 

Purtroppo ottengo risultati abbastanza strani, in modo simile ad un altro StackOverflow user who tried to use interp2d.

Il più successo che ho trovato finora utilizza inverse squares per stimare valori di Z a (X, Y). Ma la funzione non è perfetta per stimare i valori di Z vicino a Z = 0.

Cosa posso fare per ottenere una funzione z = f(x, y) dato un insieme di punti (x, y, z)? Mi manca qualcosa qui ... ho bisogno di più di una nuvola di punti per stimare in modo affidabile un valore su una superficie?

EDIT:

Questa è la funzione che ho finito per scrivere. La funzione accetta gli array di input di xs, ys, zs e interpola allo x, y utilizzando scipy.interpolate.griddata, che non richiede una griglia regolare. Sono sicuro che c'è un modo più intelligente per farlo e apprezzeremmo qualsiasi aggiornamento, ma funziona e non mi interessa le prestazioni. Incluso uno snippet nel caso in cui aiuti qualcuno in futuro.

def interpolate(x, y, xs, ys, zs): 
    r = np.sqrt(x*x + y*y) 
    t = np.arctan2(y, x) 

    rs = np.sqrt(np.multiply(xs, xs) + np.multiply(ys, ys)) 
    ts = np.arctan2(ys, xs) 

    rs = rs.ravel() 
    ts = ts.ravel() 
    zs = zs.ravel() 

    ts = np.concatenate((ts - np.pi * 2, ts, ts + np.pi * 2)) 
    rs = np.concatenate((rs, rs, rs)) 
    zs = np.concatenate((zs, zs, zs)) 


    Z = scipy.interpolate.griddata((rs, ts), zs, (r, t)) 
    Z = Z.ravel() 
    R, T = np.meshgrid(r, t) 
    return Z 

Answer 1

Stai dicendo che hai provato a utilizzare griddata. Allora, perché non funzionava? griddata funziona anche se i nuovi punti non sono regolarmente distanziati.Ad esempio,

# Definitions of xs, ys and zs 
nx, ny = 20, 30 
x = np.linspace(0, np.pi, nx) 
y = np.linspace(0, 2*np.pi, ny) 

X,Y = np.meshgrid(x, y) 

xs = X.reshape((nx*ny, 1)) 
ys = Y.reshape((nx*ny, 1)) 

## Arbitrary definition of zs 
zs = np.cos(3*xs/2.)+np.sin(5*ys/3.)**2 

## new points where I want the interpolations 
points = np.random.rand(1000, 2) 

import scipy.interpolate 
zs2 = scipy.interpolate.griddata(np.hstack((xs, ys)), zs, points) 

Non è questo quello che cerchi?

Answer 2

Se ho capito la tua domanda, si dispone di punti xs, ys, zs che sono definiti dai

xs = np.sin(U) * np.cos(V) 
ys = np.sin(U) * np.sin(V) 
zs = np.cos(U) 

Quello che vogliamo è essere in grado di interpolare e trovare un valore z per una data x e y? Perché hai bisogno di interpolazione? Le equazioni di cui sopra rappresentano una sfera, possono essere riscritti come xs*xs + ys*ys + zs*zs = 1, quindi non c'è una soluzione analitica facile a questo problema:

def Z(X, Y): 
    return np.sqrt(1-X**2-Y**2) 
    ## or return -np.sqrt(1-X**2-Y**2) since this equation has two solutions 

a meno che non ho capito male la domanda.