snapshot...

modules/icing/pirep_goes.py

@@ -5,6 +5,8 @@ import pickle
 import os
 from util.util import get_time_tuple_utc
 from aeolus.datasource import CLAVRx
+from util.geos_nav import GEOSNavigation
+import h5py
 
 clavrx_dir = '/apollo/cloud/scratch/ICING/'
 dir_fmt = '%Y_%m_%d_%j'
@@ -31,25 +33,58 @@ def get_clavrx_datasource(timestamp):
     return ds
 
 
-for idx, time in enumerate(time_keys):
-    if (idx % 4) != 0:
-        continue
-    print(100.0*(idx/len(time_keys)))
-    reports = ice_dict[time]
-    lat_s = []
-    lon_s = []
-    for tup in reports:
-        lat, lon, fl, rpt_str = tup
-        lat_s.append(lat)
-        lon_s.append(lon)
-
-    images, _, _, idxs = get_images(lon_s, lat_s, time, ['14'], [10], [1])
-    if images is not None:
-        counts, edges = np.histogram(images[0,], range=[150, 350], bins=20)
-        hist += counts
-
-f = open('/home/rink/ice_hist.pkl', 'wb')
-pickle.dump(hist, f)
-f.close()
-print('bin edges: ', edges)
+def get_grid_values(h5f, grid_name, j_c, i_c, half_width, scale_factor_name='scale_factor', add_offset_name='add_offset'):
+    hfds = h5f[grid_name]
+    attrs = hfds.attrs
 
+    grd_vals = hfds[j_c-half_width:j_c+half_width+1, i_c-half_width:i_c+half_width+1]
+    grd_vals = np.where(grd_vals == -999, np.nan, grd_vals)
+    grd_vals = np.where(grd_vals == -32768, np.nan, grd_vals)
+
+    if attrs is None:
+        return grd_vals
+
+    if scale_factor_name is not None:
+        scale_factor = attrs.get(scale_factor_name)[0]
+        grd_vals = grd_vals * scale_factor
+
+    if add_offset_name is not None:
+        add_offset = attrs.get(add_offset_name)[0]
+        grd_vals = grd_vals + add_offset
+
+    return grd_vals
+
+
+def run():
+
+    nav = GEOSNavigation(sub_lon=-75.0, CFAC=5.6E-05, COFF=-0.101332, LFAC=-5.6E-05, LOFF=0.128212, num_elems=2500, num_lines=1500)
+
+    lon_s = np.zeros(1)
+    lat_s = np.zeros(1)
+    for idx, time in enumerate(time_keys):
+        if (idx % 4) != 0:
+            continue
+        print(100.0*(idx/len(time_keys)))
+        reports = ice_dict[time]
+        for tup in reports:
+            lat, lon, fl, rpt_str = tup
+            lat_s[0] = lat
+            lon_s[0] = lon
+
+            images, _, _, idxs = get_images(lon_s, lat_s, time, ['14'], [10], [1])
+            if images is not None:
+                counts, edges = np.histogram(images[0,], range=[150, 350], bins=20)
+
+            clvr_ds = get_clavrx_datasource(time)
+            clvr_file = clvr_ds.get_file(time)
+
+            cc, ll = nav.earth_to_lc_s(lon_s, lat_s)
+
+            try:
+                h5f = h5py.File(clvr_file, 'r')
+            except Exception:
+                h5f.close()
+                print('Problem with file: ', clvr_file)
+                continue
+
+            gvals = get_grid_values(h5f, 'temp_11_0um_nom', ll[0], cc[0], 20)