Browse Source

Start converting model code to python

Petra Lamborn 5 years ago
parent
commit
178797bf0a
1 changed files with 67 additions and 0 deletions
  1. 67
    0
      py/pymodels.py

+ 67
- 0
py/pymodels.py View File

@@ -0,0 +1,67 @@
1
+from argparse import ArgumentParser, FileType
2
+import numpy as np
3
+import pandas as p
4
+import statsmodels.formula.api as smf
5
+import datetime as dt
6
+
7
+epoch = dt.datetime(2017, 1, 1)
8
+
9
+def harmonic(length, period, start = 0, harmonics = 1, prefix = ""):
10
+    x = np.arange(length) + start
11
+    y = np.ndarray((harmonics * 2, length), dtype = np.float64)
12
+    ns = []
13
+    nc = []
14
+    for i in range(harmonics):
15
+        ysi = np.sin(x/period * 2 * np.pi * (i + 1))
16
+        yci = np.cos(x/period * 2 * np.pi * (i + 1))
17
+        y[i] = ysi
18
+        ns.append("{}sin{}".format(prefix, i + 1))
19
+        y[i + harmonics] = yci
20
+        nc.append("{}cos{}".format(prefix, i + 1))
21
+    y = p.DataFrame.from_records(y.transpose(), columns = ns + nc)
22
+    return y
23
+
24
+
25
+def main():
26
+    parser = ArgumentParser(description='Harmonic models of all clusters in specified pickle')
27
+    parser.add_argument("-i", "--input",  dest="input",  
28
+                        help = "input aggregated ICP pickle path",  
29
+                        metavar="PATH", required = True,
30
+                        type = FileType('rb'))
31
+    parser.add_argument("-w", "--weather",  dest="weather",  
32
+                        help = "input weather pickle path",  
33
+                        metavar="PATH", required = True,
34
+                        type = FileType('rb'))
35
+    parser.add_argument("--weather-harmonics", dest="weather_harmonics",
36
+                        help = "number of harmonics for weather; default: 2",
37
+                        type = int, default = 2)
38
+    args = parser.parse_args()
39
+    
40
+    # print(args)
41
+
42
+    # print(harmonic(17, 8, 0, 2, "test"))
43
+    # print(epoch)
44
+
45
+    # print(p.read_pickle(args.input))
46
+    wdat = p.read_pickle(args.weather).drop(['record_no', 'station'], axis = 1)
47
+
48
+    wantedtimes = p.date_range(wdat.temp_timestamp.min(), wdat.temp_timestamp.max(), freq = "30 min")
49
+    wantedtimes = p.Series(wantedtimes.round(freq = "30 min"))
50
+    wanteddf = p.DataFrame(index = wantedtimes)
51
+    wdat = wanteddf.join(wdat.set_index(['temp_timestamp']), how = 'left').drop('temp_date', axis = 1)
52
+    wdat.index.name = "read_time"
53
+    wdat = wdat.fillna(method = 'ffill', axis = 0)
54
+    wdat['min_rolling'] = wdat['tmin_c'].rolling("1D").min() 
55
+    wdat['max_rolling'] = wdat['tmax_c'].rolling("1D").max() 
56
+    print(wdat)
57
+    wharmstart = (wdat.index.min() - epoch).total_seconds() / (60 * 30)
58
+    wharm = harmonic(wdat.shape[0], period = 365.25, start = wharmstart, 
59
+                     harmonics = args.weather_harmonics)
60
+    print(wharm)
61
+    
62
+    args.input.close()
63
+    args.weather.close()
64
+
65
+
66
+if __name__ == "__main__":
67
+    main()