|
| 1 | +"""CoDICE Hi Sectored L1A processing functions.""" |
| 2 | + |
| 3 | +import logging |
| 4 | +from pathlib import Path |
| 5 | + |
| 6 | +import numpy as np |
| 7 | +import xarray as xr |
| 8 | + |
| 9 | +from imap_processing.cdf.imap_cdf_manager import ImapCdfAttributes |
| 10 | +from imap_processing.codice import constants |
| 11 | +from imap_processing.codice.decompress import decompress |
| 12 | +from imap_processing.codice.utils import ( |
| 13 | + CODICEAPID, |
| 14 | + ViewTabInfo, |
| 15 | + apply_replacements_to_attrs, |
| 16 | + get_codice_epoch_time, |
| 17 | + get_collapse_pattern_shape, |
| 18 | + get_energy_info, |
| 19 | + get_view_tab_info, |
| 20 | + read_sci_lut, |
| 21 | +) |
| 22 | +from imap_processing.spice.time import met_to_ttj2000ns |
| 23 | + |
| 24 | +logger = logging.getLogger(__name__) |
| 25 | + |
| 26 | + |
| 27 | +def l1a_hi_sectored(unpacked_dataset: xr.Dataset, lut_file: Path) -> xr.Dataset: |
| 28 | + """ |
| 29 | + Process CoDICE Hi Sectored L1A data. |
| 30 | +
|
| 31 | + Parameters |
| 32 | + ---------- |
| 33 | + unpacked_dataset : xarray.Dataset |
| 34 | + Unpacked dataset from L0 packet file. |
| 35 | + lut_file : Path |
| 36 | + Path to the LUT file for processing. |
| 37 | +
|
| 38 | + Returns |
| 39 | + ------- |
| 40 | + xarray.Dataset |
| 41 | + Processed L1A dataset for Hi Omni data. |
| 42 | + """ |
| 43 | + # Get these values from unpacked data. These are used to |
| 44 | + # lookup in LUT table. |
| 45 | + table_id = unpacked_dataset["table_id"].values[0] |
| 46 | + view_id = unpacked_dataset["view_id"].values[0] |
| 47 | + apid = unpacked_dataset["pkt_apid"].values[0] |
| 48 | + plan_id = unpacked_dataset["plan_id"].values[0] |
| 49 | + plan_step = unpacked_dataset["plan_step"].values[0] |
| 50 | + |
| 51 | + logger.info( |
| 52 | + f"Processing species with - APID: {apid}, View ID: {view_id}, " |
| 53 | + f"Table ID: {table_id}, Plan ID: {plan_id}, Plan Step: {plan_step}" |
| 54 | + ) |
| 55 | + |
| 56 | + # ========== Get LUT Data =========== |
| 57 | + # Read information from LUT |
| 58 | + sci_lut_data = read_sci_lut(lut_file, table_id) |
| 59 | + |
| 60 | + view_tab_info = get_view_tab_info(sci_lut_data, view_id, apid) |
| 61 | + view_tab_obj = ViewTabInfo( |
| 62 | + apid=apid, |
| 63 | + view_id=view_id, |
| 64 | + sensor=view_tab_info["sensor"], |
| 65 | + three_d_collapsed=view_tab_info["3d_collapse"], |
| 66 | + collapse_table=view_tab_info["collapse_table"], |
| 67 | + ) |
| 68 | + |
| 69 | + if view_tab_obj.sensor != 1: |
| 70 | + raise ValueError("Unsupported sensor ID for Hi Sectored processing.") |
| 71 | + |
| 72 | + # ========= Get Epoch Time Data =========== |
| 73 | + # Epoch center time and delta |
| 74 | + epoch_center, deltas = get_codice_epoch_time( |
| 75 | + unpacked_dataset["acq_start_seconds"].values, |
| 76 | + unpacked_dataset["acq_start_subseconds"].values, |
| 77 | + unpacked_dataset["spin_period"].values, |
| 78 | + view_tab_obj, |
| 79 | + ) |
| 80 | + |
| 81 | + # ========= Decompress and Calculate Reshape information =========== |
| 82 | + if view_tab_obj.apid != CODICEAPID.COD_HI_SECT_SPECIES_COUNTS: |
| 83 | + raise ValueError( |
| 84 | + f"Unknown apid {view_tab_obj.apid} in Hi Sectored species processing." |
| 85 | + ) |
| 86 | + species_data = sci_lut_data["data_product_hi_tab"]["0"]["sectored"] |
| 87 | + species_names = species_data.keys() |
| 88 | + logical_source_id = "imap_codice_l1a_hi-sectored" |
| 89 | + |
| 90 | + compression_algorithm = constants.HI_COMPRESSION_ID_LOOKUP[view_tab_obj.view_id] |
| 91 | + # Decompress data using byte count information from decommed data |
| 92 | + binary_data_list = unpacked_dataset["data"].values |
| 93 | + byte_count_list = unpacked_dataset["byte_count"].values |
| 94 | + |
| 95 | + # The decompressed data in the shape of (epoch, n). Then reshape later. |
| 96 | + decompressed_data = [ |
| 97 | + decompress( |
| 98 | + packet_data[:byte_count], |
| 99 | + compression_algorithm, |
| 100 | + ) |
| 101 | + for (packet_data, byte_count) in zip( |
| 102 | + binary_data_list, byte_count_list, strict=False |
| 103 | + ) |
| 104 | + ] |
| 105 | + |
| 106 | + num_packets = len(binary_data_list) |
| 107 | + |
| 108 | + # Use chunks of (energy_x) to put data in its energy bins as done below. |
| 109 | + # Eg. [15, 15, 15, 18, 18, 15, 18, 5, 1] |
| 110 | + # where each number is energy dimension for species 'x'. |
| 111 | + species_chunk_sizes = [ |
| 112 | + len(species_data[species]["min_energy"]) for species in species_names |
| 113 | + ] |
| 114 | + |
| 115 | + # Reshape decompressed data to in below for loop: |
| 116 | + # (num_packets, num_species, energy_bins, spin_sector, inst_az) |
| 117 | + num_species = len(species_names) |
| 118 | + energy_bins = 8 |
| 119 | + collapse_shape = get_collapse_pattern_shape( |
| 120 | + sci_lut_data, |
| 121 | + view_tab_obj.sensor, |
| 122 | + view_tab_obj.collapse_table, |
| 123 | + ) |
| 124 | + if np.unique(species_chunk_sizes) != [energy_bins]: |
| 125 | + raise ValueError("Expected energy bins to be 8 for Hi Sectored data.") |
| 126 | + |
| 127 | + # Calculate collapsed size |
| 128 | + decompressed_data = np.array(decompressed_data).reshape( |
| 129 | + num_packets, num_species, energy_bins, *collapse_shape |
| 130 | + ) |
| 131 | + |
| 132 | + # ========== Create Dataset =========== |
| 133 | + cdf_attrs = ImapCdfAttributes() |
| 134 | + cdf_attrs.add_instrument_global_attrs("codice") |
| 135 | + cdf_attrs.add_instrument_variable_attrs("codice", "l1a") |
| 136 | + |
| 137 | + l1a_dataset = xr.Dataset( |
| 138 | + coords={ |
| 139 | + "epoch": xr.DataArray( |
| 140 | + met_to_ttj2000ns(epoch_center), |
| 141 | + dims=("epoch",), |
| 142 | + attrs=cdf_attrs.get_variable_attributes("epoch", check_schema=False), |
| 143 | + ), |
| 144 | + "epoch_delta_minus": xr.DataArray( |
| 145 | + deltas, |
| 146 | + dims=("epoch",), |
| 147 | + attrs=cdf_attrs.get_variable_attributes( |
| 148 | + "epoch_delta_minus", check_schema=False |
| 149 | + ), |
| 150 | + ), |
| 151 | + "epoch_delta_plus": xr.DataArray( |
| 152 | + deltas, |
| 153 | + dims=("epoch",), |
| 154 | + attrs=cdf_attrs.get_variable_attributes( |
| 155 | + "epoch_delta_plus", check_schema=False |
| 156 | + ), |
| 157 | + ), |
| 158 | + "spin_sector": xr.DataArray( |
| 159 | + np.arange(collapse_shape[0]), |
| 160 | + dims=("spin_sector",), |
| 161 | + attrs=cdf_attrs.get_variable_attributes( |
| 162 | + "spin_sector", check_schema=False |
| 163 | + ), |
| 164 | + ), |
| 165 | + "spin_sector_label": xr.DataArray( |
| 166 | + np.arange(collapse_shape[0]).astype(str), |
| 167 | + dims=("spin_sector",), |
| 168 | + attrs=cdf_attrs.get_variable_attributes( |
| 169 | + "spin_sector_label", check_schema=False |
| 170 | + ), |
| 171 | + ), |
| 172 | + "inst_az": xr.DataArray( |
| 173 | + np.arange(collapse_shape[1]), |
| 174 | + dims=("inst_az",), |
| 175 | + attrs=cdf_attrs.get_variable_attributes("inst_az", check_schema=False), |
| 176 | + ), |
| 177 | + "inst_az_label": xr.DataArray( |
| 178 | + np.arange(collapse_shape[1]).astype(str), |
| 179 | + dims=("inst_az",), |
| 180 | + attrs=cdf_attrs.get_variable_attributes( |
| 181 | + "inst_az_label", check_schema=False |
| 182 | + ), |
| 183 | + ), |
| 184 | + }, |
| 185 | + attrs=cdf_attrs.get_global_attributes(logical_source_id), |
| 186 | + ) |
| 187 | + |
| 188 | + # Final data shape of each species is (epoch, energy_h, spin_sector, inst_az) |
| 189 | + for species_index, (species_name, data) in enumerate(species_data.items()): |
| 190 | + # Add coordinate for 'energy_{species_name}' |
| 191 | + energy_centers, energy_minus, energy_plus = get_energy_info(data) |
| 192 | + coord_attrs = cdf_attrs.get_variable_attributes( |
| 193 | + "hi-energy-attrs", check_schema=False |
| 194 | + ) |
| 195 | + coord_attrs = apply_replacements_to_attrs( |
| 196 | + coord_attrs, {"species": species_name} |
| 197 | + ) |
| 198 | + l1a_dataset = l1a_dataset.assign_coords( |
| 199 | + { |
| 200 | + f"energy_{species_name}": xr.DataArray( |
| 201 | + np.array(energy_centers), |
| 202 | + dims=(f"energy_{species_name}",), |
| 203 | + attrs=coord_attrs, |
| 204 | + ) |
| 205 | + } |
| 206 | + ) |
| 207 | + # Add energy plus and minus variables |
| 208 | + minus_attrs = cdf_attrs.get_variable_attributes("hi-energy-delta-attrs") |
| 209 | + minus_attrs = apply_replacements_to_attrs( |
| 210 | + minus_attrs, {"species": species_name, "operation": "minus"} |
| 211 | + ) |
| 212 | + l1a_dataset[f"energy_{species_name}_minus"] = xr.DataArray( |
| 213 | + energy_minus, |
| 214 | + dims=(f"energy_{species_name}",), |
| 215 | + attrs=minus_attrs, |
| 216 | + ) |
| 217 | + plus_attrs = cdf_attrs.get_variable_attributes("hi-energy-delta-attrs") |
| 218 | + plus_attrs = apply_replacements_to_attrs( |
| 219 | + plus_attrs, {"species": species_name, "operation": "plus"} |
| 220 | + ) |
| 221 | + l1a_dataset[f"energy_{species_name}_plus"] = xr.DataArray( |
| 222 | + energy_plus, |
| 223 | + dims=(f"energy_{species_name}",), |
| 224 | + attrs=plus_attrs, |
| 225 | + ) |
| 226 | + |
| 227 | + # Extract species data from decompressed data: |
| 228 | + # - (num_packets, energy_bins, spin_sector, inst_az) |
| 229 | + species_data = decompressed_data[:, species_index, :, :, :] |
| 230 | + species_attrs = cdf_attrs.get_variable_attributes("hi-species-attrs") |
| 231 | + species_attrs = apply_replacements_to_attrs( |
| 232 | + species_attrs, {"species": species_name} |
| 233 | + ) |
| 234 | + # Add DEPEND_2, DEPEND_3 |
| 235 | + species_attrs["DEPEND_2"] = "spin_sector" |
| 236 | + species_attrs["LABL_PTR_2"] = "spin_sector_label" |
| 237 | + species_attrs["DEPEND_3"] = "inst_az" |
| 238 | + species_attrs["LABL_PTR_3"] = "inst_az_label" |
| 239 | + l1a_dataset[species_name] = xr.DataArray( |
| 240 | + species_data, |
| 241 | + dims=("epoch", f"energy_{species_name}", "spin_sector", "inst_az"), |
| 242 | + attrs=species_attrs, |
| 243 | + ) |
| 244 | + # Uncertainty is sqrt of counts |
| 245 | + species_unc_attrs = cdf_attrs.get_variable_attributes("hi-species-unc-attrs") |
| 246 | + species_unc_attrs = apply_replacements_to_attrs( |
| 247 | + species_unc_attrs, {"species": species_name} |
| 248 | + ) |
| 249 | + # Add DEPEND_2, DEPEND_3 |
| 250 | + species_unc_attrs["DEPEND_2"] = "spin_sector" |
| 251 | + species_unc_attrs["LABL_PTR_2"] = "spin_sector_label" |
| 252 | + species_unc_attrs["DEPEND_3"] = "inst_az" |
| 253 | + species_unc_attrs["LABL_PTR_3"] = "inst_az_label" |
| 254 | + l1a_dataset[f"unc_{species_name}"] = xr.DataArray( |
| 255 | + np.sqrt(species_data), |
| 256 | + dims=("epoch", f"energy_{species_name}", "spin_sector", "inst_az"), |
| 257 | + attrs=species_unc_attrs, |
| 258 | + ) |
| 259 | + |
| 260 | + # ========= Add Additional Variables =========== |
| 261 | + l1a_dataset["spin_period"] = xr.DataArray( |
| 262 | + unpacked_dataset["spin_period"].values * constants.SPIN_PERIOD_CONVERSION, |
| 263 | + dims=("epoch",), |
| 264 | + attrs=cdf_attrs.get_variable_attributes("spin_period"), |
| 265 | + ) |
| 266 | + l1a_dataset["data_quality"] = xr.DataArray( |
| 267 | + unpacked_dataset["suspect"].values, |
| 268 | + dims=("epoch",), |
| 269 | + attrs=cdf_attrs.get_variable_attributes("data_quality"), |
| 270 | + ) |
| 271 | + |
| 272 | + return l1a_dataset |
0 commit comments