{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.14","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":70367,"databundleVersionId":9188054,"sourceType":"competition"}],"dockerImageVersionId":30761,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## Imports","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport itertools\nimport os\nimport glob \nfrom astropy.stats import sigma_clip\n\nfrom tqdm import tqdm","metadata":{"execution":{"iopub.status.busy":"2024-09-23T15:16:58.778876Z","iopub.execute_input":"2024-09-23T15:16:58.779933Z","iopub.status.idle":"2024-09-23T15:16:58.785475Z","shell.execute_reply.started":"2024-09-23T15:16:58.779885Z","shell.execute_reply":"2024-09-23T15:16:58.784099Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"path_folder = '/kaggle/input/ariel-data-challenge-2024/' # path to the folder containing the data\npath_out = '/kaggle/tmp/data_light_raw/' # path to the folder to store the light data\noutput_dir = '/kaggle/tmp/data_light_raw/' # path for the output directory","metadata":{"execution":{"iopub.status.busy":"2024-09-23T15:16:59.574773Z","iopub.execute_input":"2024-09-23T15:16:59.575802Z","iopub.status.idle":"2024-09-23T15:16:59.580773Z","shell.execute_reply.started":"2024-09-23T15:16:59.575753Z","shell.execute_reply":"2024-09-23T15:16:59.579433Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if not os.path.exists(path_out):\n    os.makedirs(path_out)\n    print(f\"Directory {path_out} created.\")\nelse:\n    print(f\"Directory {path_out} already exists.\")","metadata":{"execution":{"iopub.status.busy":"2024-09-23T15:17:00.746628Z","iopub.execute_input":"2024-09-23T15:17:00.747115Z","iopub.status.idle":"2024-09-23T15:17:00.75453Z","shell.execute_reply.started":"2024-09-23T15:17:00.747068Z","shell.execute_reply":"2024-09-23T15:17:00.753096Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"CHUNKS_SIZE = 1","metadata":{"execution":{"iopub.status.busy":"2024-09-23T15:17:01.450303Z","iopub.execute_input":"2024-09-23T15:17:01.45071Z","iopub.status.idle":"2024-09-23T15:17:01.456816Z","shell.execute_reply.started":"2024-09-23T15:17:01.450673Z","shell.execute_reply":"2024-09-23T15:17:01.455115Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Step 1: Analog-to-Digital Conversion","metadata":{}},{"cell_type":"markdown","source":"The Analog-to-Digital Conversion (adc) is performed by the detector to convert the pixel voltage into an integer number. We revert this operation by using the gain and offset for the calibration files 'train_adc_info.csv'.","metadata":{}},{"cell_type":"markdown","source":"1. **Division by Gain**:\n   - The original ADC process multiplied the analog signal by the gain to scale it into a range that could be captured by the detector. To revert this process, we need to divide the signal by the gain, effectively reversing the scaling.\n\n2. **Adding the Offset**:\n   - In the ADC process, an offset was subtracted from the signal to adjust its baseline. To restore the signal to its original form, we add back this offset.","metadata":{}},{"cell_type":"code","source":"def ADC_convert(signal, gain, offset):\n    signal = signal.astype(np.float64)\n    signal /= gain\n    signal += offset\n    return signal","metadata":{"execution":{"iopub.status.busy":"2024-09-22T17:59:21.250293Z","iopub.execute_input":"2024-09-22T17:59:21.250695Z","iopub.status.idle":"2024-09-22T17:59:21.264311Z","shell.execute_reply.started":"2024-09-22T17:59:21.250652Z","shell.execute_reply":"2024-09-22T17:59:21.263041Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Step 2: Mask hot/dead pixel","metadata":{}},{"cell_type":"markdown","source":"- The dead pixels map is a map of the pixels that do not respond to light and, thus, can’t be accounted for any calculation. In all these frames the dead pixels are masked using python masked arrays. The bad pixels are thus masked but left uncorrected. Some methods can be used to correct bad-pixels but this task, if needed, is left to the participants.\n\n- Sigma clipping is a method used to filter out outliers in a dataset, often applied in image processing to remove noisy pixels (like \"hot\" or \"dead\" pixels). It works by iteratively calculating the mean and standard deviation of the pixel values and masking those that deviate from the mean by a certain number of standard deviations (specified by the parameter sigma). This is particularly useful for identifying and masking pixel values that are much higher or lower than the surrounding pixels due to noise.\n\n- These masks are applied to the signal to ensure that bad pixels are ignored during further analysis. However, the pixel values themselves are not corrected or replaced, as the task of correcting bad pixels is left to the participants.","metadata":{}},{"cell_type":"markdown","source":"**Explanation of the Arguments:**\n\n1. **`signal.shape[0]:`**\n    - This is the number of frames (images) in the signal array, i.e., the size along the time or sequence dimension.\n\n2. **`(signal.shape[0], 1, 1):`**\n    - This means that the hot and dead masks are replicated along the first axis (time), while the other two dimensions (height and width) remain unchanged.\n    \n**Masking:**\nWhen masking an element, it becomes an \"invalid\" or \"ignored\" element in the masked array, which is helpful in image processing to avoid artifacts from faulty (dead or hot) pixels.\n\n3. **`np.ma.masked_where:`**\n    - This function is used to \"mask\" elements in an array based on a condition. The masked elements are essentially ignored in any future computations.\n    \n    - The first argument **`dead`** or **`hot`** provides a condition (boolean mask), where True indicates that the corresponding element in signal should be masked. The second argument **`signal`** is the array to be masked.","metadata":{}},{"cell_type":"code","source":"def mask_hot_dead(signal, dead, dark):\n    # Applies sigma clipping to the dark frame\n    hot = sigma_clip(dark, sigma=5, maxiters=5).mask\n    \n    hot = np.tile(hot, (signal.shape[0], 1, 1))\n    \n    dead = np.tile(dead, (signal.shape[0], 1, 1))\n    \n    signal = np.ma.masked_where(dead, signal)\n    \n    signal = np.ma.masked_where(hot, signal)\n    \n    return signal","metadata":{"execution":{"iopub.status.busy":"2024-09-22T17:59:21.265735Z","iopub.execute_input":"2024-09-22T17:59:21.266143Z","iopub.status.idle":"2024-09-22T17:59:21.2781Z","shell.execute_reply.started":"2024-09-22T17:59:21.266104Z","shell.execute_reply":"2024-09-22T17:59:21.276624Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Step 3: Linearity Correction","metadata":{}},{"cell_type":"markdown","source":"- The goal of the linearity correction step is to correct the non-linear response of pixels due to various effects, primarily caused by capacitive leakage in the readout electronics. This non-linearity occurs because the number of photons hitting a pixel is proportional to the number of electrons collected in the well, but the pixel's electronic response to these electrons isn't perfectly linear.\n\n- The problem arises because while the number of electrons in a pixel’s well is proportional to the number of photons (with some quantum efficiency), the output signal from the pixel doesn’t increase linearly with the number of electrons. As a result, without correction, the signal would inaccurately represent the photon count, leading to distorted readings, especially for higher photon counts.\n\n- To correct this, the calibration data provides coefficients for a polynomial function that describes this non-linearity. By applying this polynomial to the raw signal, we can adjust the signal to be more representative of the actual photon count.","metadata":{}},{"cell_type":"code","source":"linear_corr = pd.read_parquet('/kaggle/input/ariel-data-challenge-2024/train/100468857/AIRS-CH0_calibration/linear_corr.parquet')","metadata":{"execution":{"iopub.status.busy":"2024-09-22T17:59:21.279796Z","iopub.execute_input":"2024-09-22T17:59:21.280329Z","iopub.status.idle":"2024-09-22T17:59:21.481412Z","shell.execute_reply.started":"2024-09-22T17:59:21.28027Z","shell.execute_reply":"2024-09-22T17:59:21.479882Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"linear_corr.head()","metadata":{"execution":{"iopub.status.busy":"2024-09-22T17:59:21.483253Z","iopub.execute_input":"2024-09-22T17:59:21.483684Z","iopub.status.idle":"2024-09-22T17:59:21.527199Z","shell.execute_reply.started":"2024-09-22T17:59:21.483639Z","shell.execute_reply":"2024-09-22T17:59:21.52583Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"linear_corr.shape","metadata":{"execution":{"iopub.status.busy":"2024-09-22T17:59:21.529074Z","iopub.execute_input":"2024-09-22T17:59:21.529554Z","iopub.status.idle":"2024-09-22T17:59:21.538122Z","shell.execute_reply.started":"2024-09-22T17:59:21.529475Z","shell.execute_reply":"2024-09-22T17:59:21.536499Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Parameters Explanation:**\n\n1. **`linear_corr`**\n    - A 3D array (degree, row, column) representing the polynomial coefficients for each pixel. Each pixel has a set of polynomial coefficients to correct its non-linearity.\n\n2. **`(signal.shape[0], 1, 1):`**\n    - This means that the hot and dead masks are replicated along the first axis (time), while the other two dimensions (height and width) remain unchanged.\n    \n3. **`np.flip(linear_corr, axis=0):`**\n    - Polynomial coefficients are often ordered from the lowest degree (constant term) to the highest degree (e.g., quadratic, cubic terms). However, NumPy's **`poly1d`** expects the coefficients in reverse order (starting with the highest degree). \n \n4. Then this code loops through every pixel **`(x, y)`** in the signal (ignoring the time dimension).\n\n5. For each pixel **`(x, y)`**, a polynomial function is created using the linear_corr coefficients for that specific pixel.\n\n6. **`np.poly1d`** generates a polynomial from the given coefficients. For example, if the coefficients for pixel (x, y) are [a, b, c], it creates the polynomial p(z).","metadata":{}},{"cell_type":"code","source":"def apply_linear_corr(linear_corr, clean_signal):\n    # Step 1: Flip the coefficients along the first axis\n    linear_corr = np.flip(linear_corr, axis=0)\n    \n    # Step 2: Loop through all pixels (ignoring the time dimension) and apply correction\n    for x, y in itertools.product(range(clean_signal.shape[1]), range(clean_signal.shape[2])):\n        # Create a polynomial object with the correction coefficients for pixel (x, y)\n        poli = np.poly1d(linear_corr[:, x, y])\n        \n        # Apply the polynomial to each pixel's signal values (along the time axis)\n        clean_signal[:, x, y] = poli(clean_signal[:, x, y])\n    \n    # Step 3: Return the corrected signal\n    return clean_signal","metadata":{"execution":{"iopub.status.busy":"2024-09-22T17:59:21.543463Z","iopub.execute_input":"2024-09-22T17:59:21.544519Z","iopub.status.idle":"2024-09-22T17:59:21.554042Z","shell.execute_reply.started":"2024-09-22T17:59:21.544434Z","shell.execute_reply":"2024-09-22T17:59:21.552587Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Step 4: Dark Current Subtraction\n\n- The dark frame is an image that represents the detector’s response when it is not exposed to light (i.e., during a very short exposure time). This image provides the baseline of the dark current for each pixel. Since dark current varies slightly from pixel to pixel, a dark frame map is used to correct the signal from each pixel individually.\n\n- The **`dead`** pixel map (passed as the dead parameter) contains a binary array where True indicates a dead pixel and False indicates a working pixel.\n\n- This line masks the corresponding **`dead pixels`** in the dark current map (dark), ensuring that these pixels are not used in the correction process.\n\n- Why is this important? Since dead pixels don’t contribute to the signal or dark current, they need to be ignored in both the signal and the dark current map to avoid introducing errors.\n\n- **`np.tile(dark, (signal.shape[0], 1, 1))`** : Replicates the Dark Current Map for All Time Frames.\n\n- The dark current **`dark`** is multiplied by the integration time (Δt), which is provided as the **`dt`** parameter. This is necessary because the dark current accumulates over time. The result is then subtracted from the **`signal`** for each pixel in each time frame.","metadata":{}},{"cell_type":"code","source":"def clean_dark(signal, dead, dark, dt):\n    # Mask dead pixels in the dark current map\n    dark = np.ma.masked_where(dead, dark)\n    \n    # Replicate the dark current map for all time frames\n    dark = np.tile(dark, (signal.shape[0], 1, 1))\n    \n    # Subtract the dark current from the signal for each time step\n    signal -= dark * dt[:, np.newaxis, np.newaxis]\n    \n    return signal","metadata":{"execution":{"iopub.status.busy":"2024-09-22T17:59:21.55603Z","iopub.execute_input":"2024-09-22T17:59:21.556604Z","iopub.status.idle":"2024-09-22T17:59:21.567923Z","shell.execute_reply.started":"2024-09-22T17:59:21.556513Z","shell.execute_reply":"2024-09-22T17:59:21.56657Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Step 5: Get Correlated Double Sampling (CDS)","metadata":{}},{"cell_type":"markdown","source":"- Correlated Double Sampling (CDS) is a technique used in detector systems to reduce certain types of noise and improve signal accuracy, particularly in environments where weak signals, like those from distant stars or exoplanets, need to be measured. It is commonly used in imaging detectors, such as CCDs (Charge-Coupled Devices) and infrared detectors, and is very useful in astrophysical observations.\n\n- **Start-of-Exposure Reading (Reset Frame):** The pixel values are read right when the exposure starts. At this point, there is little or no photon-induced signal, and the readout gives a measure of the pixel’s initial state and baseline noise.\n\n- **End-of-Exposure Reading (Signal Frame):** After the exposure time has passed, the pixel values are read again. This measurement contains the signal induced by the incoming light as well as the noise that has accumulated during the exposure.\n\n**CDS Value = End_of_Exposure Reading  -  Start_of_Exposure Reading**","metadata":{}},{"cell_type":"code","source":"def get_cds(signal):\n    cds = signal[:,1::2,:,:] - signal[:,::2,:,:]\n    return cds","metadata":{"execution":{"iopub.status.busy":"2024-09-22T17:59:21.569429Z","iopub.execute_input":"2024-09-22T17:59:21.569943Z","iopub.status.idle":"2024-09-22T17:59:21.585501Z","shell.execute_reply.started":"2024-09-22T17:59:21.569885Z","shell.execute_reply":"2024-09-22T17:59:21.584121Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Step 6: (Optional): Time Binning","metadata":{}},{"cell_type":"markdown","source":"Time binning is a method used to reduce the size of time-series data by aggregating multiple observations into single bins over a specified time interval. This step is optional and is primarily used for:\n\n- **Data Compression:** Reducing the overall size of the data by combining several individual observations into fewer \"binned\" points.\n\n- **Noise Reduction:** By averaging or summing over multiple time points, random noise can be reduced, and the signal-to-noise ratio can improve.\n\n- **Efficient Storage:** It allows storing large datasets in a more compact form, especially for long observational periods, without losing too much essential information.","metadata":{}},{"cell_type":"markdown","source":"The **`cds_signal`** is a 4D array where:\n\n- **Dimension 0:** Time or number of exposures.\n\n- **Dimension 1:** Detector grid (first axis, usually corresponding to rows).\n\n- **Dimension 2:** Detector grid (second axis, corresponding to columns).\n\n- **Dimension 3:** The pixel values (corresponding to time).\n\n\nThe output array **`cds_binned`** is initialized with a new shape where the time dimension is divided by the binning factor (i.e., **`cds_transposed.shape[1]`**  binning). This ensures that the new binned signal has fewer time steps.","metadata":{}},{"cell_type":"markdown","source":"**Perform Binning:**\n\n- The loop iterates over the binned time intervals.\n\n- For each time bin, it selects the corresponding time steps using slicing: **`i * binning : (i + 1) * binning`**.\n\n- It then sums the signal values over this time range for each pixel.\n\n- The result is stored in the cds_binned array.","metadata":{}},{"cell_type":"code","source":"def bin_obs(cds_signal, binning):\n    # Transpose the signal to move the time dimension (axis=2) to axis=3 for easier binning\n    cds_transposed = cds_signal.transpose(0, 1, 3, 2)\n    \n    # Initialize a binned signal array with reduced size in the time dimension\n    cds_binned = np.zeros((cds_transposed.shape[0], \n                           cds_transposed.shape[1] // binning, \n                           cds_transposed.shape[2], \n                           cds_transposed.shape[3]))\n\n    # Iterate over each binning interval and sum the signal over the binning window\n    for i in range(cds_transposed.shape[1] // binning):\n        cds_binned[:, i, :, :] = np.sum(cds_transposed[:, i * binning:(i + 1) * binning, :, :], axis=1)\n\n    return cds_binned","metadata":{"execution":{"iopub.status.busy":"2024-09-22T17:59:21.587956Z","iopub.execute_input":"2024-09-22T17:59:21.588503Z","iopub.status.idle":"2024-09-22T17:59:21.599363Z","shell.execute_reply.started":"2024-09-22T17:59:21.588424Z","shell.execute_reply":"2024-09-22T17:59:21.59795Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Step 7: Flat Field Correction","metadata":{}},{"cell_type":"markdown","source":"- In an ideal detector, each pixel would respond uniformly to the same amount of light. However, real detectors have imperfections:\n\n   - **Variations in sensitivity:** Different pixels might not have the same sensitivity to light due to imperfections in the manufacturing process.\n    \n   - **Quantum efficiency differences:** The ability of each pixel to convert incoming photons into a signal can vary slightly, leading to pixel-to-pixel response differences.\n\n- Flat field correction uses a flat field map, which is an image of the detector response to a uniform light source. By applying this map, we can correct these pixel-to-pixel response variations.","metadata":{}},{"cell_type":"markdown","source":"### Example: Flat Field Correction\n\nAssume you have the following pixel responses to a flat field:\n\n| Pixel  | Flat Field Value |\n|--------|------------------|\n| (0,0)  | 0.9              |\n| (0,1)  | 1.1              |\n| (1,0)  | 0.95             |\n| (1,1)  | 1.05             |\n\nFor a given signal measurement:\n\n| Pixel  | Signal |\n|--------|--------|\n| (0,0)  | 90     |\n| (0,1)  | 110    |\n| (1,0)  | 95     |\n| (1,1)  | 105    |\n\nThe corrected signal is obtained by dividing each pixel's signal by the corresponding flat field value:\n\n| Pixel  | Corrected Signal        |\n|--------|-------------------------|\n| (0,0)  | 90 / 0.9 = 100          |\n| (0,1)  | 110 / 1.1 = 100         |\n| (1,0)  | 95 / 0.95 = 100         |\n| (1,1)  | 105 / 1.05 = 100        |\n\nAfter correction, all pixels show the same response (100), meaning the detector variations have been normalized, and the signal now reflects the true incident light intensity.\n","metadata":{}},{"cell_type":"code","source":"def correct_flat_field(flat, dead, signal):\n    # Transpose the flat field and dead pixel maps to match the signal dimensions\n    flat = flat.transpose(1, 0)\n    dead = dead.transpose(1, 0)\n\n    # Mask out the dead pixels from the flat field map\n    flat = np.ma.masked_where(dead, flat)\n\n    # Replicate the flat field map for each time slice in the signal\n    flat = np.tile(flat, (signal.shape[0], 1, 1))\n\n    # Correct the signal by dividing by the flat field map\n    signal = signal / flat\n    return signal","metadata":{"execution":{"iopub.status.busy":"2024-09-22T17:59:21.601179Z","iopub.execute_input":"2024-09-22T17:59:21.601754Z","iopub.status.idle":"2024-09-22T17:59:21.611977Z","shell.execute_reply.started":"2024-09-22T17:59:21.601683Z","shell.execute_reply":"2024-09-22T17:59:21.610473Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Step 8: Calibrating all Training Data","metadata":{}},{"cell_type":"markdown","source":"You can choose to correct the non-linearity of the pixels' response, to apply flat field, dark and dead map or to leave the data unchanged. The observations are binned in time by group of 30 frames for AIRS and 360 frames for FGS1, to obtain a lighter data-cube, easier to use. The images are cut along the wavelength axis between pixels 39 and 321, so that the 282 pixels left in the wavelength dimension match the last 282 targets' points, from AIRS. The 283rd targets' point is the one for FGS1 that will be added later on.","metadata":{}},{"cell_type":"code","source":"CHUNKS_SIZE = 1","metadata":{"execution":{"iopub.status.busy":"2024-09-22T17:59:21.61441Z","iopub.execute_input":"2024-09-22T17:59:21.615126Z","iopub.status.idle":"2024-09-22T17:59:21.623588Z","shell.execute_reply.started":"2024-09-22T17:59:21.615058Z","shell.execute_reply":"2024-09-22T17:59:21.62226Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## we will start by getting the index of the training data:\ndef get_index(files,CHUNKS_SIZE ):\n    index = []\n    for file in files :\n        file_name = file.split('/')[-1]\n        if file_name.split('_')[0] == 'AIRS-CH0' and file_name.split('_')[1] == 'signal.parquet':\n            file_index = os.path.basename(os.path.dirname(file))\n            index.append(int(file_index))\n    index = np.array(index)\n    index = np.sort(index) \n    # credit to DennisSakva\n    index=np.array_split(index, len(index)//CHUNKS_SIZE)\n    \n    return index","metadata":{"execution":{"iopub.status.busy":"2024-09-22T17:59:21.625629Z","iopub.execute_input":"2024-09-22T17:59:21.626154Z","iopub.status.idle":"2024-09-22T17:59:21.63683Z","shell.execute_reply.started":"2024-09-22T17:59:21.626094Z","shell.execute_reply":"2024-09-22T17:59:21.635591Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Let's break down and understand the code step by step, describing each part and its significance.\n\n\n","metadata":{}},{"cell_type":"markdown","source":"- **Extracting All File Paths**:\n   - This part of the code extracts all the file paths within the training data directory and uses the **`get_index`** function to organize the files.\n   \n   - The **`glob.glob`** function is used to search for all files within the specified training directory, following the '*/*' pattern to match files in subdirectories. The **`get_index`** function organizes these file paths into indexed chunks.\n\n\n- **Loading Calibration Info and Axis Info**:\n   - This loads two important calibration-related files.\n   \n   \n- **Setting Configration Flags**:\n    - These are boolean flags used to control whether certain calibration steps will be applied or skipped.\n    \n    - The **`cut_inf`** and **`cut_sup`** values are used to select a specific range of wavelengths (pixels) from the signal. These values correspond to the wavelength range of interest\n    \n- **Looping Over Data Chunks**:\n    - This loop iterates over each chunk of data indexed earlier.\n    \n    - **`AIRS_CH0_clean`** and **`FGS1_clean`** are initialized as arrays to hold the cleaned signal data for AIRS-CH0 and FGS1 respectively.\n    \n    - The dimensions of the arrays correspond to:\n\n        - **`CHUNKS_SIZE`**: The number of samples in each chunk.\n        \n        - **`11250`**: The number of time steps (for AIRS).\n        \n        - **`135000`**: The number of time steps (for FGS1).\n        \n        - **`32`**: The number of detector rows.\n        \n        - **`l`**: The length of the wavelength dimension (after cutting).\n        \n\n- **Signal Preprocessing (AIRS-CH0)**:\n    - This block processes the raw signal data for AIRS-CH0.\n    \n    - The signal data is read from the parquet file and reshaped.\n    \n    - The gain and offset from **`train_adc_info`** are used to convert the signal from raw digital values to physical units using the **`ADC_convert`** function.\n    \n- **Calibrating the data: AIRS**:\n    - The signal is cropped to only include the desired wavelength range (cut_inf:cut_sup).\n    \n    - Hot and dead pixels are masked using **`mask_hot_dead`**.\n    \n    - The dark current is subtracted using the **`clean_dark function`**, which removes the constant noise generated by the detector.\n    \n- **Calibrating the data: FGS1**:\n    - Similarly Calibration is performed for the FDS1 signals as well.\n    \n- **Time Binning to reduce space**:\n    - Time binning is used to reduce the dataset size by averaging consecutive frames, thus simplifying data handling.\n    \n    - For AIRS-CH0, frames are grouped in bins of 30, while for FGS1, frames are grouped in bins of 360.\n    \n- **Flat Field Correction (AIRS and FGS1)**:\n    - Flat field correction compensates for pixel-to-pixel variations in sensitivity, ensuring uniformity across the detector.\n    \n    - The flat field correction is applied using the **`correct_flat_field function`**, which normalizes the signal based on calibration data.\n    \n- **Saving Cleaned Data**:\n    - Saving the data in a NumPy format (.npy) is efficient for further analysis, as it allows for fast loading and processing.","metadata":{}},{"cell_type":"code","source":"#Step 1: Extracting all the file paths\n\nfiles = glob.glob(os.path.join(path_folder + 'train/', '*/*'))\n\nindex = get_index(files[:2], CHUNKS_SIZE)\n\n\n\n# Step 2: Loading Calibration Info and Axis Info\n\ntrain_adc_info = pd.read_csv(os.path.join(path_folder, 'train_adc_info.csv'))\ntrain_adc_info = train_adc_info.set_index('planet_id')\naxis_info = pd.read_parquet(os.path.join(path_folder,'axis_info.parquet'))\n\n\n# Step 3: Setting Configration Flags\n\nDO_MASK = True\nDO_THE_NL_CORR = False\nDO_DARK = True\nDO_FLAT = True\nTIME_BINNING = True\n\ncut_inf, cut_sup = 39, 321\nl = cut_sup - cut_inf\n\n\n# Step 4: Looping Over Data Chunks\n\nfor n, index_chunk in enumerate(tqdm(index)):\n    AIRS_CH0_clean = np.zeros((CHUNKS_SIZE, 11250, 32, l))\n    FGS1_clean = np.zeros((CHUNKS_SIZE, 135000, 32, 32))\n    \n    # Step 5: Signal Preprocessing (AIRS-CH0)\n    \n    for i in range (CHUNKS_SIZE): \n        df = pd.read_parquet(os.path.join(path_folder,f'train/{index_chunk[i]}/AIRS-CH0_signal.parquet'))\n        signal = df.values.astype(np.float64).reshape((df.shape[0], 32, 356))\n        gain = train_adc_info['AIRS-CH0_adc_gain'].loc[index_chunk[i]]\n        offset = train_adc_info['AIRS-CH0_adc_offset'].loc[index_chunk[i]]\n        signal = ADC_convert(signal, gain, offset)\n\n        dt_airs = axis_info['AIRS-CH0-integration_time'].dropna().values\n        dt_airs[1::2] += 0.1\n        chopped_signal = signal[:, :, cut_inf:cut_sup]\n        del signal, df\n\n        #Step 6: Calibrating the data: AIRS\n        flat = pd.read_parquet(os.path.join(path_folder,f'train/{index_chunk[i]}/AIRS-CH0_calibration/flat.parquet')).values.astype(np.float64).reshape((32, 356))[:, cut_inf:cut_sup]\n        dark = pd.read_parquet(os.path.join(path_folder,f'train/{index_chunk[i]}/AIRS-CH0_calibration/dark.parquet')).values.astype(np.float64).reshape((32, 356))[:, cut_inf:cut_sup]\n        dead_airs = pd.read_parquet(os.path.join(path_folder,f'train/{index_chunk[i]}/AIRS-CH0_calibration/dead.parquet')).values.astype(np.float64).reshape((32, 356))[:, cut_inf:cut_sup]\n        linear_corr = pd.read_parquet(os.path.join(path_folder,f'train/{index_chunk[i]}/AIRS-CH0_calibration/linear_corr.parquet')).values.astype(np.float64).reshape((6, 32, 356))[:, :, cut_inf:cut_sup]\n\n        if DO_MASK:\n            chopped_signal = mask_hot_dead(chopped_signal, dead_airs, dark)\n            AIRS_CH0_clean[i] = chopped_signal\n        else:\n            AIRS_CH0_clean[i] = chopped_signal\n\n        if DO_DARK: \n            cleaned_signal = clean_dark(AIRS_CH0_clean[i], dead_airs, dark, dt_airs)\n            AIRS_CH0_clean[i] = cleaned_signal\n        else: \n            pass\n        del dark\n\n\n\n        df = pd.read_parquet(os.path.join(path_folder,f'train/{index_chunk[i]}/FGS1_signal.parquet'))\n        fgs_signal = df.values.astype(np.float64).reshape((df.shape[0], 32, 32))\n\n        FGS1_gain = train_adc_info['FGS1_adc_gain'].loc[index_chunk[i]]\n        FGS1_offset = train_adc_info['FGS1_adc_offset'].loc[index_chunk[i]]\n\n        fgs_signal = ADC_convert(fgs_signal, FGS1_gain, FGS1_offset)\n        dt_fgs1 = np.ones(len(fgs_signal))*0.1\n        dt_fgs1[1::2] += 0.1\n        chopped_FGS1 = fgs_signal\n\n        del fgs_signal, df\n        \n        \n        #Step 7: Calibrating the data: FGS1\n        \n        flat = pd.read_parquet(os.path.join(path_folder,f'train/{index_chunk[i]}/FGS1_calibration/flat.parquet')).values.astype(np.float64).reshape((32, 32))\n        dark = pd.read_parquet(os.path.join(path_folder,f'train/{index_chunk[i]}/FGS1_calibration/dark.parquet')).values.astype(np.float64).reshape((32, 32))\n        dead_fgs1 = pd.read_parquet(os.path.join(path_folder,f'train/{index_chunk[i]}/FGS1_calibration/dead.parquet')).values.astype(np.float64).reshape((32, 32))\n        linear_corr = pd.read_parquet(os.path.join(path_folder,f'train/{index_chunk[i]}/FGS1_calibration/linear_corr.parquet')).values.astype(np.float64).reshape((6, 32, 32))\n\n\n        if DO_MASK:\n            chopped_FGS1 = mask_hot_dead(chopped_FGS1, dead_fgs1, dark)\n            FGS1_clean[i] = chopped_FGS1\n        else:\n            FGS1_clean[i] = chopped_FGS1\n\n        if DO_THE_NL_CORR: \n            linear_corr_signal = apply_linear_corr(linear_corr,FGS1_clean[i])\n            FGS1_clean[i,:, :, :] = linear_corr_signal\n        del linear_corr\n\n        if DO_DARK: \n            cleaned_signal = clean_dark(FGS1_clean[i], dead_fgs1, dark,dt_fgs1)\n            FGS1_clean[i] = cleaned_signal\n        else: \n            pass\n        del dark\n        \n        \n    # SAVE DATA AND FREE SPACE\n    AIRS_cds = get_cds(AIRS_CH0_clean)\n    FGS1_cds = get_cds(FGS1_clean)\n    \n    del AIRS_CH0_clean, FGS1_clean\n    \n    \n    ## Step 8 (Optional): Time Binning to reduce space\n    if TIME_BINNING:\n        AIRS_cds_binned = bin_obs(AIRS_cds,binning=30)\n        FGS1_cds_binned = bin_obs(FGS1_cds,binning=30*12)\n    else:\n        AIRS_cds = AIRS_cds.transpose(0,1,3,2) ## this is important to make it consistent for flat fielding, but you can always change it\n        AIRS_cds_binned = AIRS_cds\n        FGS1_cds = FGS1_cds.transpose(0,1,3,2)\n        FGS1_cds_binned = FGS1_cds\n    \n    del AIRS_cds, FGS1_cds\n    \n    \n    #Step 9: Flat Field Correction (AIRS and FGS1)\n    \n    for i in range (CHUNKS_SIZE):\n        flat_airs = pd.read_parquet(os.path.join(path_folder,f'train/{index_chunk[i]}/AIRS-CH0_calibration/flat.parquet')).values.astype(np.float64).reshape((32, 356))[:, cut_inf:cut_sup]\n        flat_fgs = pd.read_parquet(os.path.join(path_folder,f'train/{index_chunk[i]}/FGS1_calibration/flat.parquet')).values.astype(np.float64).reshape((32, 32))\n        if DO_FLAT:\n            corrected_AIRS_cds_binned = correct_flat_field(flat_airs,dead_airs, AIRS_cds_binned[i])\n            AIRS_cds_binned[i] = corrected_AIRS_cds_binned\n            corrected_FGS1_cds_binned = correct_flat_field(flat_fgs, dead_fgs1, FGS1_cds_binned[i])\n            FGS1_cds_binned[i] = corrected_FGS1_cds_binned\n        else:\n            pass\n        \n        \n    ## Step 19: Saving the data\n    \n    np.save(os.path.join(path_out, 'AIRS_clean_train_{}.npy'.format(n)), AIRS_cds_binned)\n    np.save(os.path.join(path_out, 'FGS1_train_{}.npy'.format(n)), FGS1_cds_binned)\n    del AIRS_cds_binned\n    del FGS1_cds_binned","metadata":{"execution":{"iopub.status.busy":"2024-09-22T17:59:21.638938Z","iopub.execute_input":"2024-09-22T17:59:21.639603Z","iopub.status.idle":"2024-09-22T18:01:13.56932Z","shell.execute_reply.started":"2024-09-22T17:59:21.639546Z","shell.execute_reply":"2024-09-22T18:01:13.567871Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}