initial commit

modules/icing/pirep_goes.py

+from amv.intercompare import pirep_icing
+from deeplearning.amv_raob import get_images
+import numpy as np
+import pickle
+
+pirep_file = '/home/rink/data/pirep/pireps_2018040100_2018123023.csv'
+
+ice_dict, no_ice_dict = pirep_icing(pirep_file)
+print('num obs: ice, no ice', len(ice_dict), len(no_ice_dict))
+
+time_keys = list(ice_dict.keys())
+
+hist = np.zeros(20)
+
+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)
+

modules/icing/pireps.py

+import datetime
+from datetime import timezone
+import re
+
+NO_ICE = '\s*NEG\s*|\s*NONE\s*|\s*NEGATIVE\s*|\s*NO\s*'
+ICE_LVL = '\d+-\d+|FL\d+-FL\d+'
+FLT_LVL = '/FL\s{0,1}\d+\s*'
+ICE_RPT = '/IC'
+ATYPE = '/TP'
+RMK = '/RM'
+
+
+# Returns icing/non-icing
+def pirep_icing(filename, lon_range=[-180, 180], lat_range=[-55, 55]):
+    no_ice_reports = []
+    ice_reports = []
+
+    # HEADER: VALID,URGENT,AIRCRAFT,REPORT,LAT,LON
+
+    ice_dict = {}
+    no_ice_dict = {}
+
+    with open(filename) as file:
+        for idx, line in enumerate(file):
+            toks = line.split(',')
+            if toks[0] == 'VALID':  # Skip headers (concatenated files)
+                continue
+            if len(toks) != 6:  # Check for line format error
+                continue
+
+            report = toks[3]
+            lat = float(toks[4])
+            lon = float(toks[5])
+
+            if lon < lon_range[0] or lon > lon_range[1]:
+                continue
+            if lat < lat_range[0] or lat > lat_range[1]:
+                continue
+
+            # Flight level
+            so = re.search(FLT_LVL, report)
+            if so is None:
+                continue
+            else:
+                fl = so.group()
+                fl = (fl[3:].rstrip()).lstrip()
+                if fl.isnumeric():
+                    fl = float(fl) * 100 * 0.3048  # conv flight level to meters
+                else:
+                    continue
+
+            # Icing
+            ii = report.find('/IC')
+            if ii >= 0:
+                ice_str = report[ii+3:]
+                ri = ice_str.find('/RM')
+
+                if ri >= 0:
+                    ice_str = ice_str[:ri]
+                else:
+                    ri = ice_str.find(',')
+                    if ri >= 0:
+                        ice_str = ice_str[:ri]
+                ice_str = ice_str.lstrip()
+                if ice_str.find('N/A') >= 0:
+                    continue
+
+                dto = datetime.datetime.strptime(toks[0], '%Y%m%d%H%M').replace(tzinfo=timezone.utc)
+                timestmp = dto.timestamp()
+
+                if len(re.findall(NO_ICE, ice_str)) != 0:
+                    no_ice_reports.append(ice_str)
+                    rpts = no_ice_dict.get(timestmp)
+                    tup = (lat, lon, fl, ice_str)
+                    if rpts is None:
+                        rpts = []
+                        rpts.append(tup)
+                        no_ice_dict[timestmp] = rpts
+                    else:
+                        rpts.append(tup)
+                else:
+                    ice_reports.append(ice_str)
+                    so = re.search(ICE_LVL, ice_str)
+                    if so is not None:
+                        lvl_a, lvl_b = so.group().split('-')
+                        if lvl_a.find('FL') >= 0:
+                            lvl_a = float(lvl_a[2:]) * 100 * 0.3048
+                            lvl_b = float(lvl_b[2:]) * 100 * 0.3048
+                        else:
+                            lvl_a = float(lvl_a) * 100 * 0.3048
+                            lvl_b = float(lvl_b) * 100 * 0.3048
+
+                        fl = max(lvl_a, lvl_b)
+
+                    if fl > 15000.0:
+                        continue
+
+                    rpts = ice_dict.get(timestmp)
+                    tup = (lat, lon, fl, ice_str)
+                    if rpts is None:
+                        rpts = []
+                        rpts.append(tup)
+                        ice_dict[timestmp] = rpts
+                    else:
+                        rpts.append(tup)
+
+    return ice_dict, no_ice_dict