2015-08-03 4 views
5

Mam punkty danych, które reprezentują współrzędne dla tablicy 2D (macierzy). Punkty są regularnie siatkowane, z wyjątkiem braku punktów danych z niektórych pozycji siatki.Utwórz tablicę 2D Numpy ze współrzędnych

Weźmy na przykład niektóre dane XYZ, które pasują do zwykłej siatki 0,1 o kształcie (3, 4). Istnieją luki i brakujących punktów, więc istnieje 5 punktów, a nie 12:

import numpy as np 
X = np.array([0.4, 0.5, 0.4, 0.4, 0.7]) 
Y = np.array([1.0, 1.0, 1.1, 1.2, 1.2]) 
Z = np.array([3.3, 2.5, 3.6, 3.8, 1.8]) 
# Evaluate the regular grid dimension values 
Xr = np.linspace(X.min(), X.max(), np.round((X.max() - X.min())/np.diff(np.unique(X)).min()) + 1) 
Yr = np.linspace(Y.min(), Y.max(), np.round((Y.max() - Y.min())/np.diff(np.unique(Y)).min()) + 1) 
print('Xr={0}; Yr={1}'.format(Xr, Yr)) 
# Xr=[ 0.4 0.5 0.6 0.7]; Yr=[ 1. 1.1 1.2] 

Co chciałbym zobaczyć jest pokazane w ten obraz (tła: czarny = base-0 Index; szary = wartość współrzędnych; kolor = wartość matrycy, biały = brakujący).

matrix

Oto co mam, który jest intuicyjny z pętli for:

ar = np.ma.array(np.zeros((len(Yr), len(Xr)), dtype=Z.dtype), mask=True) 
for x, y, z in zip(X, Y, Z): 
    j = (np.abs(Xr - x)).argmin() 
    i = (np.abs(Yr - y)).argmin() 
    ar[i, j] = z 
print(ar) 
# [[3.3 2.5 -- --] 
# [3.6 -- -- --] 
# [3.8 -- -- 1.8]]  

Czy istnieje bardziej NumPythonic sposób vectorising podejście do powrotu tablicę 2D ar? Lub czy pętla for jest konieczna?

Odpowiedz

7

Można to zrobić w jednej linii z np.histogram2d

data = np.histogram2d(Y, X, bins=[len(Yr),len(Xr)], weights=Z) 
print(data[0]) 
[[ 3.3 2.5 0. 0. ] 
[ 3.6 0. 0. 0. ] 
[ 3.8 0. 0. 1.8]] 
1

sparse matryca jest pierwszym rozwiązaniem, które przyszło mi do głowy, ale od X i Y są pływaki, to trochę niechlujny:

In [624]: I=((X-.4)*10).round().astype(int) 
In [625]: J=((Y-1)*10).round().astype(int) 
In [626]: I,J 
Out[626]: (array([0, 1, 0, 0, 3]), array([0, 0, 1, 2, 2])) 

In [627]: sparse.coo_matrix((Z,(J,I))).A 
Out[627]: 
array([[ 3.3, 2.5, 0. , 0. ], 
     [ 3.6, 0. , 0. , 0. ], 
     [ 3.8, 0. , 0. , 1.8]]) 

nadal wymaga on, w taki czy inny sposób, aby dopasować te współrzędne z indeksami [0,1,2 ...]. Moje szybkie oszukiwanie polegało na liniowym skalowaniu wartości. Mimo to musiałem zachować ostrożność podczas konwersji pływaków na ints.

sparse.coo_matrix prace bo naturalnym sposobem zdefiniowania macierzy jest rzadki z (i, j, data) krotek, które oczywiście mogą być tłumaczone na I, J, Data list lub tablic.

Raczej podoba mi się rozwiązanie historgramowe, mimo że nie miałem okazji go użyć.

2

można wykorzystywać X i Y tworzyć X Y współrzędne na 0.1 odstępie siatki rozciągające się od min to max of X i min to max of Y a następnie wstawianie Z's do tych konkretnych pozycjach. Unikałoby to używania linspace do uzyskania Xr i Yr i jako takie musi być dość wydajne. Oto realizacja -

def indexing_based(X,Y,Z): 
    # Convert X's and Y's to indices on a 0.1 spaced grid 
    X_int = np.round((X*10)).astype(int) 
    Y_int = np.round((Y*10)).astype(int) 
    X_idx = X_int - X_int.min() 
    Y_idx = Y_int - Y_int.min() 

    # Setup output array and index it with X_idx & Y_idx to set those as Z 
    out = np.zeros((Y_idx.max()+1,X_idx.max()+1)) 
    out[Y_idx,X_idx] = Z 

    return out 

Runtime testy -

Ta sekcja porównać podejście indexing-based przeciwko drugiej np.histogram2d based solution dla Performance -

In [132]: # Create unique couples X-Y (as needed to work with histogram2d) 
    ...: data = np.random.randint(0,1000,(5000,2)) 
    ...: data1 = data[np.lexsort(data.T),:] 
    ...: mask = ~np.all(np.diff(data1,axis=0)==0,axis=1) 
    ...: data2 = data1[np.append([True],mask)] 
    ...: 
    ...: X = (data2[:,0]).astype(float)/10 
    ...: Y = (data2[:,1]).astype(float)/10 
    ...: Z = np.random.randint(0,1000,(X.size)) 
    ...: 

In [133]: def histogram_based(X,Y,Z): # From other np.histogram2d based solution 
    ...: Xr = np.linspace(X.min(), X.max(), np.round((X.max() - X.min())/np.diff(np.unique(X)).min()) + 1) 
    ...: Yr = np.linspace(Y.min(), Y.max(), np.round((Y.max() - Y.min())/np.diff(np.unique(Y)).min()) + 1) 
    ...: data = np.histogram2d(Y, X, bins=[len(Yr),len(Xr)], weights=Z) 
    ...: return data[0] 
    ...: 

In [134]: %timeit histogram_based(X,Y,Z) 
10 loops, best of 3: 22.8 ms per loop 

In [135]: %timeit indexing_based(X,Y,Z) 
100 loops, best of 3: 2.11 ms per loop