pirep_goes.py

from icing.pireps import pirep_icing
import numpy as np
import pickle
import os
from util.util import get_time_tuple_utc
from aeolus.datasource import CLAVRx, GOESL1B
from util.geos_nav import GEOSNavigation
import h5py
import datetime
from datetime import timezone
goes_date_format = '%Y%j%H'

goes16_directory = '/arcdata/goes/grb/goes16'  # /year/date/abi/L1b/RadC
clavrx_dir = '/apollo/cloud/scratch/ICING/'
dir_fmt = '%Y_%m_%d_%j'
# dir_list = [f.path for f in os.scandir('.') if f.is_dir()]
ds_dct = {}
goes_ds_dct = {}
ice_dict = None
no_ice_dict = None
time_keys = None
#pirep_file = '/home/rink/data/pireps/pireps_2019010000_2019063023.csv'
pirep_file = '/home/rink/data/pireps/pireps_20160101_20190630.csv'

ds_list = ['temp_8_5um_nom', 'temp_10_4um_nom', 'temp_11_0um_nom', 'temp_13_3um_nom', 'cld_height_acha',
           'supercooled_cloud_fraction', 'cld_temp_acha', 'cld_opd_acha', 'cloud_phase', 'solar_zenith_angle',
           'cloud_mask']
ds_grd_dct = {name: [] for name in ds_list}


def setup():
    ice_dict, no_ice_dict = pirep_icing(pirep_file)
    print('num obs: ice, no ice', len(ice_dict), len(no_ice_dict))
    return ice_dict, no_ice_dict


def get_clavrx_datasource(timestamp):
    dt_obj, time_tup = get_time_tuple_utc(timestamp)
    date_dir_str = dt_obj.strftime(dir_fmt)
    ds = ds_dct.get(date_dir_str)
    if ds is None:
        ds = CLAVRx(clavrx_dir + date_dir_str + '/')
        ds_dct[date_dir_str] = ds
    return ds


def get_goes_datasource(timestamp):
    dt_obj, time_tup = get_time_tuple_utc(timestamp)

    yr_dir = str(dt_obj.timetuple().tm_year)
    date_dir = dt_obj.strftime(dir_fmt)

    files_path = goes16_directory + '/' + yr_dir + '/' + date_dir + '/abi' + '/L1b' + '/RadC/'
    ds = goes_ds_dct.get(date_dir)
    if ds is None:
        ds = GOESL1B(files_path)
        goes_ds_dct[date_dir] = ds
    return ds


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
    ylen, xlen = hfds.shape
    j_l = j_c-half_width
    i_l = i_c-half_width
    if j_l < 0 or i_l < 0:
        return None
    j_r = j_c+half_width+1
    i_r = i_c+half_width+1
    if j_r >= ylen or i_r >= xlen:
        return None

    grd_vals = hfds[j_l:j_r, i_l:i_r]
    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(pirep_dct, outfile=None):
    time_keys = list(pirep_dct.keys())

    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)
    last_clvr_file = None
    last_h5f = None

    cnt = 0
    miss_a = 0
    miss_b = 0
    for idx, time in enumerate(time_keys):
        reports = pirep_dct[time]
        for tup in reports:
            lat, lon, fl, rpt_str = tup
            lat_s[0] = lat
            lon_s[0] = lon

            try:
                clvr_ds = get_clavrx_datasource(time)
            except Exception:
                miss_a += 1
                continue

            clvr_file = clvr_ds.get_file(time)[0]
            if clvr_file is None:
                miss_b += 1
                continue

            if clvr_file != last_clvr_file:
                try:
                    h5f = h5py.File(clvr_file, 'r')
                except Exception:
                    if h5f is not None:
                        h5f.close()
                    print('Problem with file: ', clvr_file)
                    continue
                if last_h5f is not None:
                    last_h5f.close()
                last_h5f = h5f
                last_clvr_file = clvr_file
            else:
                h5f = last_h5f

            cc, ll = nav.earth_to_lc_s(lon_s, lat_s)

            for didx, ds_name in enumerate(ds_list):
                gvals = get_grid_values(h5f, ds_name, ll[0], cc[0], 20)
                if gvals is not None:
                    ds_grd_dct[ds_name].append(gvals)

            cnt += 1

    print('num images: ', cnt, len(time_keys), miss_a, miss_b)


def analyze(ice_dct, no_ice_dct):

    last_file = None
    ice_files = []
    ice_times = []
    for ts in list(ice_dct.keys()):
        try:
            ds = get_goes_datasource(ts)
            goes_file, t_0, _ = ds.get_file(ts)
            if goes_file is not None and goes_file != last_file:
                ice_files.append(goes_file)
                ice_times.append(t_0)
                last_file = goes_file
        except Exception:
            continue

    last_file = None
    no_ice_files = []
    no_ice_times = []
    for ts in list(no_ice_dct.keys()):
        try:
            ds = get_goes_datasource(ts)
            goes_file, t_0, _ = ds.get_file(ts)
            if goes_file is not None and goes_file != last_file:
                no_ice_files.append(goes_file)
                no_ice_times.append(t_0)
                last_file = goes_file
        except Exception:
            continue

    ice_times = np.array(ice_times)
    no_ice_times = np.array(no_ice_times)

    itrsct_vals, comm1, comm2 = np.intersect1d(no_ice_times, ice_times, return_indices=True)
    print(itrsct_vals.shape, comm1.shape, comm2.shape)


def create_file(filename, ds_list, ds_types):
    pass