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() contrail_idxs = (mask_image == 1).flatten() print('number of contrail pixels: ', np.sum(contrail_idxs)) # Assuming GOES FD for now ------------------- # elems, lines = np.meshgrid(np.arange(5424), np.arange(5424)) # lines, elems = lines.flatten(), elems.flatten() # See note above regarding support for GOES-18 # contrail_lines, contrail_elems = lines[contrail_idxs], elems[contrail_idxs] # 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)