-
Notifications
You must be signed in to change notification settings - Fork 3
/
regular_gridding.py
178 lines (146 loc) · 5.92 KB
/
regular_gridding.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
import numpy as N
from matplotlib.mlab import griddata
import matplotlib
import matplotlib.cbook
import matplotlib.pyplot as plt
from numpy.random import uniform
def edges(C):
E = 0.5*(C[:-1]+C[1:])
return N.hstack((C[0]-2*(E[0]-C[0]), E, C[-1]+2*(C[-1]-E[-1])))
def regularlyGrid(xarr, yarr, zarr, xstart=None, xfinal=None, xstep=None, ystart=None, yfinal=None, ystep=None):
"Returns the regularly grided xi, yi, and zi arrays from the initial data."
# if xstart,xfinal,xstep,ystart,yfinal,ystep are NOT given, they are derived from the data
if xstart == None:
xstart = xarr.min()
if xfinal == None:
xfinal = xarr.max()
if xstep == None:
xstep = 1.0 * (xfinal - xstart) / len(xarr)
if ystart == None:
ystart = yarr.min()
if yfinal == None:
yfinal = yarr.max()
if ystep == None:
ystep = 1.0 * (yfinal - ystart) / len(yarr)
xi = N.arange(xstart, xfinal, xstep)
yi = N.arange(ystart, yfinal, ystep)
'''
Programming note:
linspace does (start, final, how many steps)
arange does (start, final, step)
'''
# grid the data.
print len(xarr), len(yarr), len(zarr), len(xi), len(yi)
xarr,yarr,zarr=remove_duplicates(xarr,yarr,zarr)
zi = griddata(xarr, yarr, zarr, xi, yi)
print "done gridding"
return xi, yi, zi
def remove_duplicates(xarr, yarr, zarr):
"""Removes duplicate points along the given axes for 2d plotting. In dataflow
only removes the points temporarily for the plotting, thus changing the
desired x and y axes will not result in a loss of data."""
uniques = [] #list of row indices to skip (ie unique pts to keep in datafile)
dups = []
numrows = len(xarr)
for index in range(numrows):
dups.append(range(numrows)) #list of lists. Each inner list is a list of every row index (0 to len)
#when a row becomes distinct from another, its index is removed from the the other's index
#for field in distinct:
for i in range(numrows):
if not uniques.__contains__(i): #if the row isn't unique
for j in range(i + 1, numrows):
if not uniques.__contains__(j) and dups[i].__contains__(j):
if xarr[i] != xarr[j] or yarr[i] != yarr[j]:
dups[i].remove(j)
dups[j].remove(i)
if len(dups[i]) == 1:
#if every row in the column is distinct from the ith row, then it is unique
#MAY NOT NEED TO KEEP TRACK OF UNIQUES...
#CAN ALWAYS GET BY len(dups[i])==1
uniques.append(i)
if len(uniques) == numrows:
#if all rows are deemed unique, return
return xarr,yarr,zarr
#ALL UNIQUE ROWS ARE INDEXED IN uniques NOW
rows_to_be_removed = []
for alist in dups:
#average the detector counts of every detector of each row in alist
#and save the resulting averages into the first row of alist.
#then delete other rows, ie remove them
if len(alist) == 1:
#if the row is unique skip this list
pass
else:
for k in range(1, len(alist)):
zarr[alist[0]] = (zarr[alist[0]] + zarr[alist[k]]) / 2.0
dups[alist[k]] = [alist[k]] #now the kth duplicate set of indices has only k, so it will be skipped
rows_to_be_removed.append(alist[k])
rows_to_be_removed.sort() #duplicate rows to be removed indices in order now
rows_to_be_removed.reverse() #duplicate rows to be removed indices in reverse order now
for k in rows_to_be_removed:
xarr = N.delete(xarr,k,0)
yarr = N.delete(yarr,k,0)
zarr = N.delete(zarr,k,0)
#update primary detector dimension
return xarr,yarr,zarr
def plotGrid(xi, yi, zi):
# contour the gridded data, plotting dots at the randomly spaced data points.
CS = plt.contour(xi, yi, zi, 15, linewidths=0.5, colors='k')
CS = plt.contourf(xi, yi, zi, 15, cmap=plt.cm.jet)
plt.colorbar() # draw colorbar
# plot data points.
plt.scatter(xarr, yarr, marker='o', c='b', s=5)
plt.xlim(xstart, xfinal) #setting xlimits
plt.ylim(ystart, yfinal) #setting ylimits
plt.title('Regular Grid')
plt.show()
def regularlyGridRandom():
"Makes a contour and contourf plot of randomly generated data pts."
# make up some randomly distributed data
npts = 1000
x = uniform(-3, 3, npts)
y = uniform(-3, 3, npts)
z = x * N.exp(-x ** 2 - y ** 2)
# define grid.
xi = N.arange(-3.1, 3.1, 0.05)
yi = N.arange(-3.1, 3.1, 0.05)
# grid the data.
zi = griddata(x, y, z, xi, yi)
# contour the gridded data, plotting dots at the randomly spaced data points.
CS = plt.contour(xi, yi, zi, 15, linewidths=0.5, colors='k')
CS = plt.contourf(xi, yi, zi, 15, cmap=plt.cm.jet)
plt.colorbar() # draw colorbar
# plot data points.
plt.scatter(x, y, marker='o', c='b', s=5)
plt.xlim(-3, 3)
plt.ylim(-3, 3)
plt.title('griddata test (%d points)' % npts)
plt.show()
#returning data
#results = [xi, yi, zi]
#return results
def pcolorRandom():
"Makes a pcolormesh plot of randomly generated data pts."
# make up some randomly distributed data
npts = 100
x = uniform(-3, 3, npts)
y = uniform(-3, 3, npts)
z = x * N.exp(-x ** 2 - y ** 2)
# define grid.
xi = N.arange(-3.1, 3.1, 0.05)
yi = N.arange(-3.1, 3.1, 0.05)
# grid the data.
zi = griddata(x, y, z, xi, yi)
# contour the gridded data, plotting dots at the randomly spaced data points.
plt.pcolormesh(xi, yi, zi)
#CS = plt.contour(xi,yi,zi,15,linewidths=0.5,colors='k')
#CS = plt.contourf(xi,yi,zi,15,cmap=plt.cm.jet)
plt.colorbar() # draw colorbar
# plot data points.
plt.scatter(x, y, marker='o', c='b', s=5)
plt.xlim(-3, 3)
plt.ylim(-3, 3)
plt.title('griddata test (%d points)' % npts)
plt.show()
if __name__ == "__main__":
regularlyGridRandom()