-
Notifications
You must be signed in to change notification settings - Fork 13
/
Copy pathDataGridder.py
189 lines (175 loc) · 6.7 KB
/
DataGridder.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
import numpy as np
import random
class DataGridder:
"""
Examine a list of x,y coordinates to see if they lie on a grid, and
calculate the grid order and shape if they do
"""
def __init__(self, x=None, y=None):
self.setData(x, y)
def setData(self, x, y):
"""
Set the x,y values to be tested
"""
self._x = None if x is None else x.astype(np.float64)
self._y = None if y is None else y.astype(np.float64)
self._tested = False
self._isValid = None
self._gridOrder = None
self._gridShape = None
def calcGrid(self):
"""
Calculate the grid for the current x,y values.
Returns gridshape, gridorder
gridshape is the shape (nrow,ncol) of the grid
index is the index for reordering if required, otherwise None
If the points do not lie on a grid returns gridshape as None
The grid x and y values are calculated as
xg=x[gridorder].reshape(gridshape)
xg=y[gridorder].reshape(gridshape)
if gridorder is not None, else simply reshaping x and y
"""
if self._x is None or self._y is None:
return None, None
if not self._tested:
self._tested = True
if len(self._x) > 4:
try:
self._tryGridOrdered()
valid = self._gridIsValid()
if not valid:
self._tryGridReorder()
valid = self._gridIsValid()
if not valid:
self._gridShape = None
self._gridOrder = None
self._isValid = valid
except:
pass
return self._gridShape, self._gridOrder
def _tryGridOrdered(self):
"""
Test if the data as ordered form a regular grid
"""
# Simple test assumes grid is ordered, so dot product of vector from
# first node to second with difference from n to n+1 should be positive
# except when n to n+1 jumps back to start of a new row. So negative
# dot products should identify start of rows, and be equally spaced in
# data set.
x = self._x
y = self._y
l = len(x)
xd = np.diff(x)
yd = np.diff(y)
ld = l - 1
dotprod = xd[0 : ld - 1] * xd[1:ld] + yd[0 : ld - 1] * yd[1:ld]
ends = np.flatnonzero(dotprod < 0) + 2
if len(ends) < 2:
return False
nr = ends[0]
nc = len(ends) // 2 + 1
if (nr >= 2) and (nc >= 2) and (nr * nc == l) and (all(ends % nr <= 1)):
self._gridShape = (nc, nr)
self._gridOrder = None
return True
return False
def _tryGridReorder(self):
"""
Test for grid by reordering rows and columns
"""
# More complex test. Attempts to identify grid rows by identifying
# start and end of an outside edge. Does this by finding the point
# furthest from the centre, then the point furthest from it (assumed
# to be across the diagonal, and finally the point furthest from the
# line between them, assumed to be the opposite end of an edge. Then
# uses the vector direction of the edge to sort the points into rows
# and columns. Finally checks that all the resultant grid cells are
# valid, meaning not re-entrant.
x = self._x
y = self._y
xm = np.mean(x)
ym = np.mean(y)
i0 = np.argmax(np.square(x - xm) + np.square(y - ym))
x0 = x[i0]
y0 = y[i0]
u = x - x0
v = y - y0
i3 = np.argmax(u * u + v * v)
du1 = u[i3]
dv1 = v[i3]
offset = u * dv1 - v * du1
i1 = np.argmax(offset)
i2 = np.argmin(offset)
scl = np.sqrt(du1 * du1 + dv1 * dv1)
u /= scl
v /= scl
# Transform to set the corners i0,i1,i2,i3
# to (-1,-1),(-1,1),(1,-1),(1,1). Note that
# uv(i0) = (0,0)
m = np.array(
[
[1.0, 0, 0, 0],
[1.0, u[i1], v[i1], u[i1] * v[i1]],
[1.0, u[i2], v[i2], u[i2] * v[i2]],
[1.0, u[i3], v[i3], u[i3] * v[i3]],
]
)
coefs = np.linalg.solve(m, np.array([[-1, -1], [1, -1], [-1, 1], [1, 1]]))
prms = np.vstack((np.ones(u.shape), u, v, u * v))
uv = prms.T.dot(coefs)
u = uv[:, 0]
v = uv[:, 1]
# Now guess at a spacing that will separate the
# first row from the second...
vrow = 1 + np.min(v[v + 1 > np.abs(u + 1)])
npt = u.shape[0]
# Count the elements in the first row
ncol = v[v < (vrow / 2 - 1)].shape[0]
nrow = int(npt / ncol)
# Test that row count is a divisor of the
# number of elements
if nrow * ncol != npt:
return False
# Add values to each u for each row (based on sorting of v)
# so that u values will sort into elements across each row
# in turn.
uwid = (np.max(u) - np.min(u)) * 2.0
iu = np.argsort(v)
urow = np.floor_divide(np.arange(npt), ncol) * uwid
t = np.vstack((u[iu], v[iu], urow)).T
urow1 = u
urow1[iu] += urow
# Now have ig, which is the sort order for converting into a
# grid
self._gridOrder = np.argsort(urow1)
self._gridShape = (nrow, ncol)
return True
def _gridIsValid(self):
if not self._gridShape:
return False
u = self._x.copy() if self._gridOrder is None else self._x[self._gridOrder]
v = self._y.copy() if self._gridOrder is None else self._y[self._gridOrder]
u = u.reshape(self._gridShape)
v = v.reshape(self._gridShape)
dudu = u[1:, :] - u[:-1, :]
dvdu = v[1:, :] - v[:-1, :]
dudv = u[:, 1:] - u[:, :-1]
dvdv = v[:, 1:] - v[:, :-1]
# Test the cross product at each of the internal angles in turn
# They should all have the same sign...
dotprod = dudu[:, :-1] * dvdv[:-1, :] - dvdu[:, :-1] * dudv[:-1, :]
imax = np.argmax(np.abs(dotprod))
sign = 1.0 if dotprod.flat[imax] > 0 else -1.0
dotprod *= sign
if np.any(dotprod < 0):
return False
dotprod = sign * (dudu[:, 1:] * dvdv[:-1, :] - dvdu[:, 1:] * dudv[:-1, :])
if np.any(dotprod < 0):
return False
dotprod = sign * (dudu[:, 1:] * dvdv[1:, :] - dvdu[:, 1:] * dudv[1:, :])
if np.any(dotprod < 0):
return False
dotprod = sign * (dudu[:, :-1] * dvdv[1:, :] - dvdu[:, :-1] * dudv[1:, :])
if np.any(dotprod < 0):
return False
return True