import datetime
from datetime import timezone
import numpy as np
import re
from pathlib import Path
import h5py
NO_ICE = '\s*NEG\s*|\s*NONE\s*|\s*NEGATIVE\s*|\s*NO\s*'
ICE_LVL = '\d+-\d+|FL\d+-FL\d+'
# ICE Intensity: 0:No intensity report, 1:Trace, 2:Light, 3:Light Moderate, 4:Moderate, 5:Moderate Severe, 6:Severe
# NO/NEG -1
# For Binary Observations, NO/NEG: -1, Any Icing: 1
ICE_TRACE = '\s*TRACE\s*|\s*TRC\s*|\s*TR\s*'
ICE_LGT = '\s*LGT\s*|\s*LIGHT\s*|\s*LT\s*|\s*LGHT\s*|\s*LIT\s*|\s*LITE\s*|\s*LG\s*'
ICE_MOD = '\s*MOD\s*|\s*MDT\s*|\s*MODERATE\s*|\s*HEAVY\s*'
ICE_SVR = '\s*SVR\s*|\s*SEV\s*|\s*SEVERE\s*'
FLT_LVL = '/FL\s{0,1}\d+\s*'
ICE_RPT = '/IC'
ATYPE = '/TP'
RMK = '/RM'
def time_sort_pirep_dict(the_dict):
new_dict = {}
keys = []
tidx = []
for ts in list(the_dict.keys()):
rpts = the_dict[ts]
for idx, tup in enumerate(rpts):
keys.append(ts)
tidx.append(idx)
keys = np.array(keys)
tidx = np.array(tidx)
sidxs = np.argsort(keys)
keys = keys[sidxs]
tidx = tidx[sidxs]
for idx, key in enumerate(keys):
rpts = the_dict[key]
tup = rpts[tidx[idx]]
n_rpts = new_dict.get(key)
if n_rpts is None:
n_rpts = []
new_dict[key] = n_rpts
n_rpts.append(tup)
return new_dict
# Returns icing, no_icing (no icing observed), neg-icing (icing expected, but not observed)
# from Iowa Mesonet Archived PIREP reports obtained in csv format
def pirep_icing(filename, lon_range=[-180, 180], lat_range=[-63, 63]):
# HEADER: VALID,URGENT,AIRCRAFT,REPORT,LAT,LON
ice_dict = {}
neg_ice_dict = {}
no_ice_dict = {}
cnt_ice = 0
cnt_no_ice = 0
cnt_neg_ice = 0
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
if len(toks) != 9: # Check for line format error
continue
report = toks[3]
try:
# lat = float(toks[4])
# lon = float(toks[5])
lat = float(toks[7])
lon = float(toks[8])
except Exception:
continue
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
dto = datetime.datetime.strptime(toks[0], '%Y%m%d%H%M').replace(tzinfo=timezone.utc)
timestmp = dto.timestamp()
# 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
if len(re.findall(NO_ICE, ice_str)) != 0:
rpts = neg_ice_dict.get(timestmp)
tup = (lat, lon, fl, -1, cnt_neg_ice, ice_str)
cnt_neg_ice += 1
if rpts is None:
rpts = []
rpts.append(tup)
neg_ice_dict[timestmp] = rpts
else:
rpts.append(tup)
else:
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
# Intensity
I = 0
so = re.search(ICE_TRACE, ice_str)
if so is not None:
I = 1
else:
so = re.search(ICE_LGT, ice_str)
if so is not None:
so = re.search(ICE_MOD, ice_str)
if so is not None:
I = 3
else:
I = 2
else:
so = re.search(ICE_MOD, ice_str)
if so is not None:
so = re.search(ICE_SVR, ice_str)
if so is not None:
I = 5
else:
I = 4
else:
so = re.search(ICE_SVR, ice_str)
if so is not None:
I = 6
rpts = ice_dict.get(timestmp)
tup = (lat, lon, fl, I, cnt_ice, ice_str)
cnt_ice += 1
if rpts is None:
rpts = []
rpts.append(tup)
ice_dict[timestmp] = rpts
else:
rpts.append(tup)
else:
rpts = no_ice_dict.get(timestmp)
tup = (lat, lon, fl, -1, cnt_no_ice, 'no observed icing')
cnt_no_ice += 1
if rpts is None:
rpts = []
rpts.append(tup)
no_ice_dict[timestmp] = rpts
else:
rpts.append(tup)
ice_dict = time_sort_pirep_dict(ice_dict)
no_ice_dict = time_sort_pirep_dict(no_ice_dict)
neg_ice_dict = time_sort_pirep_dict(neg_ice_dict)
return ice_dict, no_ice_dict, neg_ice_dict
# Returns icing, no_icing (no icing observed)
def pirep_icing_boeing_txt(filename, lon_range=[-180, 180], lat_range=[-63, 63]):
# year, month, day, hour, minute, sec, lat, lon, alt, degC, airspeed, icing
ice_dict = {}
no_ice_dict = {}
cnt_ice = 0
cnt_no_ice = 0
with open(filename) as file:
for idx, line in enumerate(file):
toks = line.split(' ')
if len(toks) != 12: # Check for line format error
continue
year = 2000 + int(toks[0])
month = int(toks[1])
day = int(toks[2])
hour = int(toks[3])
minute = int(toks[4])
second = int(toks[5])
try:
lat = float(toks[6])
lon = float(toks[7])
except Exception:
continue
if lon < lon_range[0] or lon > lon_range[1]:
continue
if lat < lat_range[0] or lat > lat_range[1]:
continue
flt_alt = float(toks[8]) * 0.3048 # convert feet to meters
degK = float(toks[9]) + 273.0
spd = float(toks[10])
ice_bool = int(toks[11])
dto = datetime.datetime(year=year, month=month, day=day, hour=hour, minute=minute, second=second)
dto = dto.replace(tzinfo=timezone.utc)
ts = dto.timestamp()
if ice_bool == 1:
tup = (lat, lon, flt_alt, 1, 0, 'None', degK, spd)
rpts = ice_dict.get(ts)
if rpts is None:
rpts = []
rpts.append(tup)
ice_dict[ts] = rpts
else:
print('duplicate: ', tup)
print('duplicate: ', rpts[0])
print('Skip -------------')
# rpts.append(tup)
cnt_ice += 1
else:
tup = (lat, lon, flt_alt, -1, 0, 'None', degK, spd)
rpts = no_ice_dict.get(ts)
if rpts is None:
rpts = []
rpts.append(tup)
no_ice_dict[ts] = rpts
else:
print('duplicate: ', tup)
print('duplicate: ', rpts[0])
print('Skip -------------')
# rpts.append(tup)
cnt_no_ice += 1
print('icing, no icing cnt: ', cnt_ice, cnt_no_ice)
ice_dict = time_sort_pirep_dict(ice_dict)
no_ice_dict = time_sort_pirep_dict(no_ice_dict)
return ice_dict, no_ice_dict, None
def pirep_icing_from_boeing_nc(dirname='/ships19/cloud/scratch/ICING/BOEING/icing/2020/', lon_range=[-180, 180], lat_range=[-63, 63]):
flist = []
for path in Path(dirname).rglob('*.nc'):
flist.append(path)
ice_dict = {}
no_ice_dict = {}
cnt_ice = 0
cnt_no_ice = 0
for fname in flist:
try:
h5f = h5py.File(fname, 'r')
ice_bool = h5f['ice_bool'][:]
lons = h5f['longitude'][:]
lats = h5f['latitude'][:]
years = 2000 + h5f['year'][:]
mons = h5f['month'][:]
days = h5f['day'][:]
hours = h5f['hour'][:]
mins = h5f['minute'][:]
secs = h5f['second'][:]
flt_levs = (h5f['altitude_ft'][:]).astype('float32')
flt_levs *= 0.3048 # convert feet to meters
except Exception as e:
print(e)
h5f.close()
continue
nrpts = ice_bool.shape[0]
if nrpts > 0:
for k in range(nrpts):
if lons[k] < lon_range[0] or lons[k] > lon_range[1]:
continue
if lats[k] < lat_range[0] or lats[k] > lat_range[1]:
continue
dto = datetime.datetime(year=years[k], month=mons[k], day=days[k], hour=hours[k], minute=mins[k], second=secs[k])
dto = dto.replace(tzinfo=timezone.utc)
ts = dto.timestamp()
if ice_bool[k] == 1:
tup = (lats[k], lons[k], flt_levs[k], 1, 0, 'None')
rpts = ice_dict.get(ts)
if rpts is None:
rpts = []
rpts.append(tup)
ice_dict[ts] = rpts
else:
rpts.append(tup)
cnt_ice += 1
else:
tup = (lats[k], lons[k], flt_levs[k], -1, 0, 'None')
rpts = no_ice_dict.get(ts)
if rpts is None:
rpts = []
rpts.append(tup)
no_ice_dict[ts] = rpts
else:
rpts.append(tup)
cnt_no_ice += 1
h5f.close()
print('num icing/no icing: ', cnt_ice, cnt_no_ice)
ice_dict = time_sort_pirep_dict(ice_dict)
no_ice_dict = time_sort_pirep_dict(no_ice_dict)
return ice_dict, no_ice_dict, None