2015-08-03 9 views
5

Ich habe Datenpunkte, die Koordinaten für ein 2D-Array (Matrix) darstellen. Die Punkte werden regelmäßig gerastert, außer dass an einigen Gitterpositionen Datenpunkte fehlen.Machen Sie 2D Numpy Array von Koordinaten

Betrachten Sie zum Beispiel einige XYZ-Daten, die auf ein reguläres 0.1-Raster mit Form passen (3, 4). Es gibt Lücken und fehlende Punkte, so gibt es 5 Punkte, und nicht 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] 

Was Ich mag würde, um zu sehen ist in diesem Bild (Hintergründen dargestellt: schwarz = Basis-Index 0, grau = Wert koordinieren; Farbe = Matrixwert; Weiß = fehlt).

matrix

Hier ist, was ich habe, das mit einer for-Schleife ist intuitiv:

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]]  

Gibt es eine NumPythonic Weise den Ansatz der Vektorisierung ar ein 2D-Array zurück? Oder ist die for-Schleife notwendig?

Antwort

7

Sie können es mit 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

Die sparse Matrix ist die erste Lösung, die in den Sinn kam, aber da X und Y Schwimmer sind, ist es ein wenig chaotisch:

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]]) 

Es muss noch, in der einen oder anderen, diese Koordinaten übereinstimmen mit [0,1,2 ...] Indizes. Mein schneller Betrug bestand darin, die Werte linear zu skalieren. Trotzdem musste ich vorsichtig sein, wenn ich Floats in Ints umwandelte.

sparse.coo_matrix funktioniert, weil eine natürliche Art und Weise eine spärliche Matrix zu definieren, mit (i, j, data) Tupel ist, die natürlich Data Listen oder Arrays I, J, übersetzt werden.

Ich mag eher die Historgram-Lösung, obwohl ich keine Gelegenheit hatte, sie zu verwenden.

2

Sie können XY und verwenden die X-Y-Koordinaten auf einem 0.1 beabstandetes Gitter vom Z's in diese spezifischen Positionen Einfügen min to max of X und min to max of Y und dann erstreckt zu schaffen. Dies würde die Verwendung von linspace vermeiden, um Xr und Yr zu bekommen, und als solches muss ziemlich effizient sein. Hier ist die Umsetzung -

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 Tests -

In diesem Abschnitt wird die indexing-based Ansatz gegen die andere np.histogram2d based solution für die Leistung vergleichen -

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 
Verwandte Themen