From fff88000da16bafdaab0c5a0fa5f13f175ee76e1 Mon Sep 17 00:00:00 2001 From: tomrink <rink@ssec.wisc.edu> Date: Wed, 22 May 2024 11:17:52 -0500 Subject: [PATCH] snapshot... --- modules/contrail/util.py | 76 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 76 insertions(+) create mode 100644 modules/contrail/util.py diff --git a/modules/contrail/util.py b/modules/contrail/util.py new file mode 100644 index 00000000..e9517f62 --- /dev/null +++ b/modules/contrail/util.py @@ -0,0 +1,76 @@ +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) + + + + + -- GitLab