From 41b8bc6ee31c68579676429ceac51f6b70243a8e Mon Sep 17 00:00:00 2001
From: tomrink <rink@ssec.wisc.edu>
Date: Wed, 3 Mar 2021 11:37:34 -0600
Subject: [PATCH] initial commit
---
modules/icing/__init__.py | 0
modules/icing/pirep_goes.py | 36 ++++++++++++
modules/icing/pireps.py | 107 ++++++++++++++++++++++++++++++++++++
3 files changed, 143 insertions(+)
create mode 100644 modules/icing/__init__.py
create mode 100644 modules/icing/pirep_goes.py
create mode 100644 modules/icing/pireps.py
diff --git a/modules/icing/__init__.py b/modules/icing/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/modules/icing/pirep_goes.py b/modules/icing/pirep_goes.py
new file mode 100644
index 00000000..2e05a235
--- /dev/null
+++ b/modules/icing/pirep_goes.py
@@ -0,0 +1,36 @@
+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)
+
diff --git a/modules/icing/pireps.py b/modules/icing/pireps.py
new file mode 100644
index 00000000..61053b00
--- /dev/null
+++ b/modules/icing/pireps.py
@@ -0,0 +1,107 @@
+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
--
GitLab