|
| 1 | +#!/usr/bin/env python3 |
| 2 | +import sys |
| 3 | +#sys.path.append('/work/noaa/da/xinjin/jedi/bufr-query/build/lib/python3.10') |
| 4 | +#sys.path.append('/work/noaa/da/xinjin/jedi/ioda-bundle/build/lib/python3.10') |
| 5 | +#sys.path.append('/work2/noaa/da/nesposito/backend_20240701/ioda-bundle/build/lib/python3.10/pyioda') |
| 6 | +#sys.path.append('/work2/noaa/da/nesposito/backend_20240701/ioda-bundle/build/lib/python3.10/pyiodaconv') |
| 7 | + |
| 8 | +import bufr |
| 9 | +import netCDF4 as nc |
| 10 | +import os |
| 11 | +from pyioda.ioda.Engines.Bufr import Encoder |
| 12 | +#from wxflow import Logger |
| 13 | +import logging |
| 14 | +#from bufr2ioda_ncep_1bmusa import Bufr2IodaAmusa |
| 15 | +#from bufr2ioda_ncep_esamua import Bufr2IodaEsamusa |
| 16 | + |
| 17 | +# logger = Logger(os.path.basename(__file__), level='INFO', colored_log=True) |
| 18 | +logger = logging.getLogger(__name__) |
| 19 | +YAML_NORMAL = False # current as normal need remap for path2/1bama |
| 20 | + |
| 21 | +R1000 = 1000.0 |
| 22 | +R1000000 = 1000000.0 |
| 23 | +INVALID = R1000 |
| 24 | +# Cosmic background temperature. Taken from Mather,J.C. et. al., 1999, "Calibrator Design for the COBE |
| 25 | +# Far-Infrared Absolute Spectrophotometer (FIRAS)"Astrophysical Journal, vol 512, pp 511-520 |
| 26 | +TSPACE = 2.7253 |
| 27 | + |
| 28 | +nc_dir = '/work2/noaa/da/xinjin/gdas-validation/global-workflow/sorc/gdas.cd/ush/ioda/bufr2ioda' |
| 29 | + |
| 30 | + |
| 31 | +class ACCoeff: |
| 32 | + def __init__(self, ac_dir, sat_id='n19'): |
| 33 | + file_name = os.path.join(ac_dir, 'amsua_' + sat_id + '.ACCoeff.nc') |
| 34 | + nc_file = nc.Dataset(file_name) |
| 35 | + self.n_fovs = len(nc_file.dimensions['n_FOVs']) |
| 36 | + self.n_channels = len(nc_file.dimensions['n_Channels']) |
| 37 | + self.a_earth = nc_file.variables['A_earth'][:] |
| 38 | + self.a_platform = nc_file.variables['A_platform'][:] |
| 39 | + self.a_space = nc_file.variables['A_space'][:] |
| 40 | + self.a_ep = self.a_earth + self.a_platform |
| 41 | + self.a_sp = self.a_space * TSPACE |
| 42 | + |
| 43 | + |
| 44 | +def remove_ant_corr(i, ac, ifov, t): |
| 45 | + # AC: Structure containing the antenna correction coefficients for the sensor of interest. |
| 46 | + # iFOV: The FOV index for a scanline of the sensor of interest. |
| 47 | + # T: On input, this argument contains the brightness |
| 48 | + |
| 49 | + # print(f't before corr: {t[:100]}') |
| 50 | + t = ac.a_ep[i, ifov] * t + ac.a_sp[i, ifov] |
| 51 | + # print(f't after corr: {t[:100]}') |
| 52 | + t[(ifov < 1) | (ifov > ac.n_fovs)] = [INVALID] |
| 53 | + return t |
| 54 | + |
| 55 | + |
| 56 | +def apply_ant_corr(i, ac, ifov, t): |
| 57 | + # t: on input, this argument contains the antenna temperatures for the sensor channels. |
| 58 | + t = (t - ac.a_sp[i, ifov]) / ac.a_ep[i, ifov] |
| 59 | + t[(ifov < 1) | (ifov > ac.n_fovs)] = [INVALID] |
| 60 | + return t |
| 61 | + |
| 62 | + |
| 63 | +def get_description(yaml_path): |
| 64 | + |
| 65 | + description = bufr.encoders.Description(yaml_path) |
| 66 | + |
| 67 | + return description |
| 68 | + |
| 69 | + |
| 70 | +def apply_corr(sat_id, ta, ifov): |
| 71 | + # condition to do this? TODO check |
| 72 | + ac = ACCoeff(nc_dir) # TODO add later |
| 73 | + if sat_id not in ['n15', 'n16']: |
| 74 | + # Convert antenna temperature to brightness temperature |
| 75 | + ifov = ifov.astype(int) - 1 |
| 76 | + for i in range(ta.shape[1]): |
| 77 | + logger.info(f'inside loop for allpy ta to tb: i = {i}') |
| 78 | + x = ta[:, i] |
| 79 | + # logger.info(f'ta before correction: {x[:100]}') |
| 80 | + if YAML_NORMAL: |
| 81 | + x = apply_ant_corr(i, ac, ifov, x) |
| 82 | + else: |
| 83 | + x = remove_ant_corr(i, ac, ifov, x) |
| 84 | + # logger.info(f'ta after correction: {x[:100]}') |
| 85 | + x[x >= R1000] = R1000000 |
| 86 | + ta[:, i] = x |
| 87 | + return ta |
| 88 | + |
| 89 | + |
| 90 | +def re_map_variable(container): |
| 91 | + # read_bufrtovs.f90 |
| 92 | + # antcorr_application.f90 |
| 93 | + # search the keyword “ta2tb” for details |
| 94 | + sat_ids = container.all_sub_categories() |
| 95 | + logger.info(sat_ids) |
| 96 | + for sat_id in sat_ids: |
| 97 | + logger.info(f'Converting for {sat_id}, ...') |
| 98 | + print(f'Converting for {sat_id}, ...') |
| 99 | + ta = container.get('variables/brightnessTemperature', sat_id) |
| 100 | + if ta.shape[0]: |
| 101 | + ifov = container.get('variables/fieldOfViewNumber', sat_id) |
| 102 | + logger.info(f'ta before correction1: {ta[:100, :]}') |
| 103 | + tb = apply_corr(sat_id, ta, ifov) |
| 104 | + logger.info(f'tb after correction1: {tb[:100, :]}') |
| 105 | + container.replace('variables/brightnessTemperature', tb, sat_id) |
| 106 | + |
| 107 | + |
| 108 | +def get_one_data(input_path, yaml_path, category): |
| 109 | + cache = bufr.DataCache.has(input_path, yaml_path) |
| 110 | + if cache: |
| 111 | + logger.info(f'The cache existed get data container from it') |
| 112 | + container = bufr.DataCache.get(input_path, yaml_path) |
| 113 | + #data = Encoder(get_description(yaml_path)).encode(container)[(category,)] |
| 114 | + #bufr.DataCache.mark_finished(input_path, yaml_path, [category]) |
| 115 | + else: |
| 116 | + # If cacache does not exist, get data into cache |
| 117 | + # Get data info container first |
| 118 | + logger.info(f'The cache is not existed') |
| 119 | + container = bufr.Parser(input_path, yaml_path).parse() |
| 120 | + #logger.info(f'add original container list into a cache = {container.list()}') |
| 121 | + #bufr.DataCache.add(input_path, yaml_path, container.all_sub_categories(), container) |
| 122 | + #data = Encoder(get_description(yaml_path)).encode(container)[(category,)] |
| 123 | + #bufr.DataCache.mark_finished(input_path, yaml_path, [category]) |
| 124 | + return cache, container |
| 125 | + |
| 126 | + |
| 127 | +def mark_one_data(cache, input_path, yaml_path, category, container=None): |
| 128 | + if cache: |
| 129 | + logger.info(f'The cache existed get data container from it') |
| 130 | + bufr.DataCache.mark_finished(input_path, yaml_path, [category]) |
| 131 | + else: |
| 132 | + logger.info(f'add original container list into a cache = {container.list()}') |
| 133 | + bufr.DataCache.add(input_path, yaml_path, container.all_sub_categories(), container) |
| 134 | + bufr.DataCache.mark_finished(input_path, yaml_path, [category]) |
| 135 | + |
| 136 | + |
| 137 | +def create_obs_group(input_path1, input_path2, yaml_1b, yaml_es, category): |
| 138 | + logger.info(f'imput_path: {input_path1}, {input_path2}, and category: {category}') |
| 139 | + logger.info(f'Entering function to create obs group for {category} with yaml path {yaml_es} and {yaml_1b}') |
| 140 | + cache_1, container_1 = get_one_data(input_path1, yaml_es, category) |
| 141 | + cache_2, container_2 = get_one_data(input_path2, yaml_1b, category) |
| 142 | + |
| 143 | + re_map_variable(container_2) |
| 144 | + |
| 145 | + print(f'emily checking - Original container list = ', container_1.list()) |
| 146 | + print(f'emily checking - Original container list = ', container_2.list()) |
| 147 | + container = container_1 |
| 148 | + container.append(container_2) |
| 149 | + |
| 150 | + print('Done') |
| 151 | + data = Encoder(get_description(yaml_es)).encode(container)[(category,)] |
| 152 | + mark_one_data(cache_1, input_path1, yaml_es, category, container=container_1) |
| 153 | + mark_one_data(cache_2, input_path2, yaml_1b, category, container=container_2) |
| 154 | + return data |
| 155 | + |
0 commit comments