snapshot...

modules/contrail/util.py

+import re
+import datetime
+from datetime import timezone
+
+import numpy as np
+import xarray as xr
+import rasterio
+from PIL import Image
+import matplotlib.pyplot as plt
+import matplotlib.image as mpimg
+import h5py
+from util.util import get_grid_values_all
+from util.gfs_reader import *
+from util.geos_nav import GEOSNavigation
+from aeolus.datasource import GFSfiles
+
+gfs_files = GFSfiles('/Users/tomrink/data/contrail/gfs/')
+# GEOSNavigation needs to be updated to support GOES-18
+# geos_nav = GEOSNavigation()
+
+
+def load_image(image_path):
+    # Extract date time string from image path
+    datetime_regex = '_\\d{8}_\\d{6}'
+    datetime_string = re.search(datetime_regex, image_path)
+    if datetime_string:
+        datetime_string = datetime_string.group()
+    dto = datetime.datetime.strptime(datetime_string, '_%Y%m%d_%H%M%S').replace(tzinfo=timezone.utc)
+    ts = dto.timestamp()
+
+    img = mpimg.imread(image_path)
+
+    return img, ts
+
+
+def get_contrail_mask_image(image, thresh=0.157):
+    image = np.where(image > thresh,1, 0)
+    return image
+
+
+def extract(mask_image, image_ts, clavrx_path):
+    gfs_file, _, _ = gfs_files.get_file(image_ts)
+    gfs_h5f = h5py.File(gfs_file, 'r')
+    xr_dataset = xr.open_dataset(gfs_file)
+
+    clvrx_h5f = h5py.File(clavrx_path, 'r')
+    cloud_top_press = get_grid_values_all(clvrx_h5f, 'cld_press_acha').flatten()
+    clvrx_lons = get_grid_values_all(clvrx_h5f, 'longitude').flatten()
+    clvrx_lats = get_grid_values_all(clvrx_h5f, 'latitude').flatten()
+
+    # Assuming GOES FD for now -------------------
+    elems, lines = np.meshgrid(np.arange(5424), np.arange(5424))
+    lines, elems = lines.flatten(), elems.flatten()
+
+    contrail_idxs = (mask_image == 1).flatten()
+    print('number of contrail pixels: ', np.sum(contrail_idxs))
+    contrail_lines, contrail_elems = lines[contrail_idxs], elems[contrail_idxs]
+
+    # See note above regarding support for GOES-18
+    # contrail_lons, contrail_lats = geos_nav.lc_to_earth(contrail_elems, contrail_lines)
+
+    contrail_press = cloud_top_press[contrail_idxs]
+    contrail_lons, contrail_lats = clvrx_lons[contrail_idxs], clvrx_lats[contrail_idxs]
+
+    keep = np.invert(np.isnan(contrail_press))
+    contrail_press = contrail_press[keep]
+    contrail_lons = contrail_lons[keep]
+    contrail_lats = contrail_lats[keep]
+
+    wind = get_point_s(xr_dataset, ['u-wind','v-wind'], contrail_lons, contrail_lats, contrail_press)
+    print(wind)
+
+
+
+
+