minor

modules/icing/pirep_goes.py

@@ -11,6 +11,8 @@ import re
 import datetime
 from datetime import timezone
 import glob
+from skyfield.api import load
+from skyfield import api, almanac
 
 goes_date_format = '%Y%j%H'
 goes16_directory = '/arcdata/goes/grb/goes16'  # /year/date/abi/L1b/RadC
@@ -20,7 +22,7 @@ dir_fmt = '%Y_%m_%d_%j'
 ds_dct = {}
 goes_ds_dct = {}
 #pirep_file = '/home/rink/data/pireps/pireps_2019010000_2019063023.csv'
-pirep_file = '/home/rink/data/pireps/pireps_20180101_20200331.csv'
+pirep_file = '/Users/tomrink/data/pirep/pireps_20180101_20200331.csv'
 
 l1b_ds_list = ['temp_10_4um_nom', 'temp_11_0um_nom', 'temp_12_0um_nom', 'temp_13_3um_nom', 'temp_3_75um_nom',
                'temp_6_2um_nom', 'temp_6_7um_nom', 'temp_7_3um_nom', 'temp_8_5um_nom', 'temp_9_7um_nom',
@@ -1354,3 +1356,27 @@ def spatial_filter(icing_dict):
             grd_bins[lat_idx, lon_idx] += 1
 
     return grd_bins
+
+
+def analyze_moon_phase(icing_dict):
+    ts = api.load.timescale()
+    eph = api.load('de421.bsp')
+
+    last_date = None
+    moon_phase = None
+    cnt = 0
+
+    for key in list(icing_dict.keys()):
+        dt_obj, dt_tup = get_time_tuple_utc(key)
+        date = datetime.date(dt_tup.tm_year, dt_tup.tm_mon, dt_tup.tm_mday)
+        if last_date != date:
+            t = ts.utc(dt_tup.tm_year, dt_tup.tm_mon, dt_tup.tm_mday)
+            moon_phase = almanac.moon_phase(eph, t)
+            if 30 < moon_phase.degrees < 330:
+                cnt += 1
+            last_date = date
+        else:
+            if 30 < moon_phase.degrees < 330:
+                cnt += 1
+
+    print(len(icing_dict), cnt)
\ No newline at end of file