-
Notifications
You must be signed in to change notification settings - Fork 1
/
fourier_full.py
120 lines (90 loc) · 3.16 KB
/
fourier_full.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
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import blackman
maxLength = 3561
dt = 1.
periods = 1. / np.fft.rfftfreq(maxLength, d=dt)
# find the minimum length and shorten all the values to be within that range
flowData = pd.read_csv("universallyAlignedGlobalFlow_DailyQ2_column.csv")
newDataDict = {}
minColLength = 10000000
for col in flowData.columns:
mask = np.array(~flowData[col].isna())
if np.sum(mask) < minColLength:
minColLength = np.sum(mask)
newData = np.array(flowData[col])[mask]
newDataDict[col] = newData
for key in newDataDict.keys():
newDataDict[key] = newDataDict[key][:minColLength]
flowData = pd.DataFrame.from_dict(newDataDict)
maxLength = minColLength #flowData[]#11326
print(maxLength / (2 * np.pi))
dt = 1.
#periods = np.fft.rfftfreq(maxLength, d=dt) * 2 * np.pi
#print(periods)
maxLength = periods.shape[0]
#print(periods)
#quit()
#flowData = pd.read_csv("universallyAlignedGlobalFlow_DailyQ2_column.csv")
print(flowData.columns)
catchmentToFourier = {}
catchmentToFourier["period_days"] = periods
print(catchmentToFourier["period_days"])
#quit()
#print(1. / periods)
lengths = []
l = 0
for column in flowData:
if column.isnumeric():
data = flowData[column]
# get rid of everything past the first "none"
# get the first bit of data:
startIndex = 0
for i in range(len(data)):
if ~np.isnan(data[i]):
startIndex = i
# print(startIndex)
break
stopIndex = len(data)
for i in range(startIndex, len(data)):
# print(data[i])
# sdata = pd.Series(data[i])
# if pd.isna(sdata[0]):
# stopIndex = i
# break
if np.isnan(data[i]):
stopIndex = i
break
keptData = data[startIndex:stopIndex]
# print(keptData)
# keptData = pd.Series(keptData)
keptData = np.array([float(x) for x in keptData])
keptData = keptData - np.mean(keptData)
keptData = keptData / np.max(keptData)
catchment = column
w = blackman(keptData.shape[0])
result = np.fft.rfft(keptData * w)
#real = np.real(result) ** 2
#imaginary = np.imag(result) ** 2
magnitudes = np.abs(result) ** 2
#fig, axs = plt.subplots(1, 2)
#axs[0].plot(keptData)
#axs[1].plot(1. / periods[:magnitudes.shape[0]], magnitudes)
#plt.show()
lengths.append(magnitudes.shape[0])
diff = maxLength - magnitudes.shape[0]
extra = [None] * diff
magnitudes = np.concatenate((magnitudes, extra), axis=0)
catchmentToFourier[catchment] = magnitudes
#print(catchmentToFourier)
percentDone = float(l) / float(len(flowData.columns) - 5)
print(percentDone)
l = l + 1
#for key in catchmentToFourier.keys():
# print(len(catchmentToFourier[key]))
for cat in catchmentToFourier.keys():
print(cat, len(catchmentToFourier[cat]))
df = pd.DataFrame.from_dict(catchmentToFourier)
df.to_csv("glob_flow_fourier_decomposition_blackman.csv")
print(max(lengths))