snapshot...

modules/icing/pirep_goes.py

@@ -4,7 +4,7 @@ import pickle
 import matplotlib.pyplot as plt
 import os
 from util.util import get_time_tuple_utc, GenericException, add_time_range_to_filename, is_night, is_day, \
-    get_grid_values_all, check_oblique
+    get_grid_values_all, check_oblique, make_times, find_bin_index
 from aeolus.datasource import CLAVRx, GOESL1B
 from util.geos_nav import GEOSNavigation
 import h5py
@@ -1522,6 +1522,32 @@ def time_filter(icing_dct, dt_str_0=None, dt_str_1=None, format_code='%Y-%m-%d_%
     return keep_times, keep_reports
 
 
+# dt_str_0: start datetime string in format YYYY-MM-DD_HH:MM (default)
+# dt_str_1: end datetime string in format YYYY-MM-DD_HH:MM (default)
+# format_code: Python Datetime format code, default: '%Y-%m-%d_%H:%M'
+# return a flatten list of icing reports
+def time_filter_2(times, dt_str_0=None, dt_str_1=None, format_code='%Y-%m-%d_%H:%M'):
+    ts_0 = 0
+    if dt_str_0 is not None:
+        dto_0 = datetime.datetime.strptime(dt_str_0, format_code).replace(tzinfo=timezone.utc)
+        ts_0 = dto_0.timestamp()
+
+    ts_1 = np.finfo(np.float64).max
+    if dt_str_1 is not None:
+        dto_1 = datetime.datetime.strptime(dt_str_1, format_code).replace(tzinfo=timezone.utc)
+        ts_1 = dto_1.timestamp()
+
+    keep_idxs = []
+    keep_times = []
+
+    for ts, idx in enumerate(times):
+        if ts_0 <= ts < ts_1:
+            keep_times.append(ts)
+            keep_idxs.append(idx)
+
+    return keep_times, keep_idxs
+
+
 def analyze_moon_phase(icing_dict):
     ts = api.load.timescale()
     eph = api.load('de421.bsp')
@@ -1546,6 +1572,20 @@ def analyze_moon_phase(icing_dict):
     print(len(icing_dict), cnt)
 
 
+def collect(icing_times, icing_intensity):
+    nbins = 6
+    jdays = ['071', '072', '074', '079', '080', '081', '082', '084', '085', '088', '107', '110', '113',
+             '117', '119', '129', '132', '134', '139', '150', '164', '203', '205', '232', '252', '254']
+
+    keep_idxs = [[[] for j in range(nbins)] for i in range(len(jdays))]
+
+    for jd in jdays:
+        fltr_times, fltr_idxs = time_filter_2(icing_times, '2019-'+jd+'_00:00', '2019-'+jd+'_00:00', format_code='%Y-%j_%H:%M')
+        _, ts_edges = make_times('2019-'+jd+'_00:00', None, num_steps=6, format_code='%Y-%j_%H:%M', hours=4)
+        bin_idxs = find_bin_index(np.array(ts_edges), np.array(fltr_times))
+
+
+
 # ------------ This code will not be needed when we implement a Fully Connected CNN -----------------------------------
 # Example GOES file to retrieve GEOS parameters in MetPy form (CONUS)
 exmp_file = '/Users/tomrink/data/OR_ABI-L1b-RadC-M6C14_G16_s20193140811215_e20193140813588_c20193140814070.nc'

modules/util/util.py

@@ -185,20 +185,24 @@ def value_to_index(nda, value):
     return idx
 
 
-def find_bin_index(nda, value):
+def find_bin_index(nda, value_s):
     idxs = np.arange(nda.shape[0])
+    iL_s = np.zeros(value_s.shape[0])
+    iL_s[:,] = -1
 
-    above = value >= nda
+    for k, v in enumerate(value_s):
+        above = v >= nda
         if not above.any():
-        return -1
+            continue
 
-    below = value < nda
+        below = v < nda
         if not below.any():
-        return -1
+            continue
 
         iL = idxs[above].max()
+        iL_s[k] = iL
 
-    return iL
+    return iL_s
 
 
 # array solzen must be degrees, missing values must NaN. For small roughly 50x50km regions only