{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"name":"python","version":"3.10.13","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":30746,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Calibrating and Time Binning Astronomical Data","metadata":{}},{"cell_type":"markdown","source":"このNotebookは以下のNotebookを日本語訳したものです。自分でわかる様に意訳をしているため、ところどころ原文と表現が違うところがありますがご了承ください。  \n※ 本Notebookは自分がコンペの理解を深めることを目的として作成しています。  \n　 間違っているところが多くあるかと思いますので、お気軽にコメントください！修正します。\n\nhttps://www.kaggle.com/code/gordonyip/update-calibrating-and-binning-astronomical-data","metadata":{}},{"cell_type":"markdown","source":"## Importing Package","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport polars as pl\nimport os\nimport glob\nimport itertools\nimport gc\nfrom astropy.stats import sigma_clip\n\nfrom matplotlib import pyplot as plt\n\nfrom tqdm import tqdm","metadata":{"execution":{"iopub.status.busy":"2024-08-17T14:41:32.216228Z","iopub.execute_input":"2024-08-17T14:41:32.216686Z","iopub.status.idle":"2024-08-17T14:41:33.470516Z","shell.execute_reply.started":"2024-08-17T14:41:32.216643Z","shell.execute_reply":"2024-08-17T14:41:33.469271Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Define path of data, out","metadata":{}},{"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\n\n# path_folder = 'ariel-data-challenge-2024/' # path to the local folder containing the data\n# path_out = 'tmp/data_light_raw/' # path to the folder to store the light data\n# output_dir = 'tmp/data_light_raw/' # path for the output directory","metadata":{"execution":{"iopub.status.busy":"2024-08-17T14:41:33.473753Z","iopub.execute_input":"2024-08-17T14:41:33.474454Z","iopub.status.idle":"2024-08-17T14:41:33.480921Z","shell.execute_reply.started":"2024-08-17T14:41:33.474408Z","shell.execute_reply":"2024-08-17T14:41:33.479615Z"},"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-08-17T14:41:33.482453Z","iopub.execute_input":"2024-08-17T14:41:33.482891Z","iopub.status.idle":"2024-08-17T14:41:33.494505Z","shell.execute_reply.started":"2024-08-17T14:41:33.482846Z","shell.execute_reply":"2024-08-17T14:41:33.493285Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Data Preprocessing","metadata":{}},{"cell_type":"markdown","source":"### Define chunk_size to avoid exceeding the memory capacity","metadata":{}},{"cell_type":"code","source":"CHUNKS_SIZE = 6 ","metadata":{"execution":{"iopub.status.busy":"2024-08-17T14:41:33.496153Z","iopub.execute_input":"2024-08-17T14:41:33.496601Z","iopub.status.idle":"2024-08-17T14:41:33.504543Z","shell.execute_reply.started":"2024-08-17T14:41:33.496557Z","shell.execute_reply":"2024-08-17T14:41:33.503307Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Step 1: Analog-to-Digital Conversion\nアナログーデジタル変換(ADC)は、ピクセルの電圧を整数値に変換するために検出機によって行われる。この操作を戻すために、キャリブレーションファイルである「train_adc_info.csv」に含まれるゲインとオフセットを使用する。\n\nオフセット誤差、ゲイン誤差を有効にした場合の値は以下の様に求められる。\n\n$$\n補正後の値 = (\\text{補正前の値} - \\text{offset}) \\times \\text{gain}\n$$\n今回は補正前の値が欲しいので、以下の関数で計算している。","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-08-17T14:41:33.508282Z","iopub.execute_input":"2024-08-17T14:41:33.508657Z","iopub.status.idle":"2024-08-17T14:41:33.518066Z","shell.execute_reply.started":"2024-08-17T14:41:33.508626Z","shell.execute_reply":"2024-08-17T14:41:33.516815Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Step 2: Mask hot/dead pixel\nDead pixel(pixel error) : 1つのフレームにある１個以上のプクセルが、収録した画像を正しく表示しない状態。\n→正しく表示できないので、前処理でマスクすべき。\n\n","metadata":{}},{"cell_type":"code","source":"def mask_hot_dead(signal, dead, dark):\n    hot = sigma_clip(\n        dark, sigma=5, maxiters=5\n    ).mask\n    hot = np.tile(hot, (signal.shape[0], 1, 1))\n    dead = np.tile(dead, (signal.shape[0], 1, 1))\n    signal = np.ma.masked_where(dead, signal)\n    signal = np.ma.masked_where(hot, signal)\n    return signal","metadata":{"execution":{"iopub.status.busy":"2024-08-17T14:41:33.519826Z","iopub.execute_input":"2024-08-17T14:41:33.520191Z","iopub.status.idle":"2024-08-17T14:41:33.53208Z","shell.execute_reply.started":"2024-08-17T14:41:33.52016Z","shell.execute_reply":"2024-08-17T14:41:33.530943Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Step 3: Linearity Correction\n#### ピクセル応答の非線形性について\nピクセル応答は本来線形性を持つ。以下の要因でピクセル応答が非線形となりうる。\n1. 容量漏れ（Capacity Leakage）\n    多くの電子が集まった時に、キャパシタにかかる電圧も大きくなる。電圧の容量が想定より大きくなった場合に応答が非線形になりうる。\n1. 飽和効果（Saturation Effect）\n    ウェル（光を受けて電子が集まる場所）が飽和状態となった時に、光子によって生じる電子が十分に蓄積されず、ピクセル応答が鈍くなる。→応答が非線形になる。\n他にも、量子効率の変動、温度依存性などがある。\n\nこれらの非線形効果を補正するために以下の関数を定義する。\n","metadata":{}},{"cell_type":"code","source":"def apply_linear_corr(linear_corr,clean_signal):\n    linear_corr = np.flip(linear_corr, axis=0)\n    for x, y in itertools.product(\n                range(clean_signal.shape[1]), range(clean_signal.shape[2])\n            ):\n        poli = np.poly1d(linear_corr[:, x, y])\n        clean_signal[:, x, y] = poli(clean_signal[:, x, y])\n    return clean_signal","metadata":{"execution":{"iopub.status.busy":"2024-08-17T14:41:33.533807Z","iopub.execute_input":"2024-08-17T14:41:33.534168Z","iopub.status.idle":"2024-08-17T14:41:33.546281Z","shell.execute_reply.started":"2024-08-17T14:41:33.534137Z","shell.execute_reply":"2024-08-17T14:41:33.545192Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Step ４: Dark Current Subtraction\n\nダークカレントは、入射光とは無関係に積分時間中に各ピクセルに蓄積される一定の信号を表す。積分時間に比例してダークカレントが蓄積して行く。\n\n提供されたデータには、ダークフレームやデッドピクセルマップなどのキャリブレーションファイルが含まれている。これらを使用して観測データを前処理することができます。ダークフレームは、非常に短い露光時間で撮影された検出器の応答を示すマップであり、これを用いて検出器のダークカレントを補正する。\n\n補正された画像は、次のようにして得る。\n\n$$\n\\text{image} - \\text{dark} \\times \\Delta t\n$$\n\nここで、ダークカレントマップは最初にデッドピクセルの補正が行われる。\n","metadata":{}},{"cell_type":"code","source":"def clean_dark(signal, dead, dark, dt):\n\n    dark = np.ma.masked_where(dead, dark)\n    dark = np.tile(dark, (signal.shape[0], 1, 1))\n\n    signal -= dark* dt[:, np.newaxis, np.newaxis]\n    return signal","metadata":{"execution":{"iopub.status.busy":"2024-08-17T14:41:33.547807Z","iopub.execute_input":"2024-08-17T14:41:33.548222Z","iopub.status.idle":"2024-08-17T14:41:33.559458Z","shell.execute_reply.started":"2024-08-17T14:41:33.548166Z","shell.execute_reply":"2024-08-17T14:41:33.558117Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Step ５: Get Correlated Double Sampling (CDS)\nサイエンスフレームは、露光の開始時と終了時の間で交互に取得される。読み出しの方式は、二重サンプリングを伴うランプであり、これをCorrelated Double Sampling (CDS)と呼ぶ。検出器は、露光の開始時と終了時の2回読み出す。最終的なCDSは、以下の差分として得られる。\n\n$$\n\\text{CDS} = (\\text{End of exposure}) - (\\text{Start of exposure})\n$$\n\n疑問  \nサイエンスフレームってなに？？別名があるんじゃないか？","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-08-17T14:41:33.560957Z","iopub.execute_input":"2024-08-17T14:41:33.562021Z","iopub.status.idle":"2024-08-17T14:41:33.571214Z","shell.execute_reply.started":"2024-08-17T14:41:33.561978Z","shell.execute_reply":"2024-08-17T14:41:33.570117Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Step ６ (Optional): Time Binning\nこのステップは主にメモリを節約するために実行される。  \n時系列観測データが指定された頻度でまとめてビンにされる。\n\n→メモリが足りないときはこちらを適用すると、学習に使うことができる様になるかも？？","metadata":{}},{"cell_type":"code","source":"def bin_obs(cds_signal,binning):\n    cds_transposed = cds_signal.transpose(0,1,3,2)\n    cds_binned = np.zeros((cds_transposed.shape[0], cds_transposed.shape[1]//binning, cds_transposed.shape[2], cds_transposed.shape[3]))\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    return cds_binned","metadata":{"execution":{"iopub.status.busy":"2024-08-17T14:41:33.572427Z","iopub.execute_input":"2024-08-17T14:41:33.572738Z","iopub.status.idle":"2024-08-17T14:41:33.584456Z","shell.execute_reply.started":"2024-08-17T14:41:33.572712Z","shell.execute_reply":"2024-08-17T14:41:33.583261Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Step 7: Flat Field Correction\nフラットフィールドは、検出器の均一な光子に対する応答を示すマップであり、ピクセルごとの検出器の変動、例えば各ピクセルの異なる量子効率を補正するために使用される。\n","metadata":{}},{"cell_type":"code","source":"def correct_flat_field(flat,dead, signal):\n    flat = flat.transpose(1, 0)\n    dead = dead.transpose(1, 0)\n    flat = np.ma.masked_where(dead, flat)\n    flat = np.tile(flat, (signal.shape[0], 1, 1))\n    signal = signal / flat\n    return signal","metadata":{"execution":{"iopub.status.busy":"2024-08-17T14:41:33.586014Z","iopub.execute_input":"2024-08-17T14:41:33.586541Z","iopub.status.idle":"2024-08-17T14:41:33.601861Z","shell.execute_reply.started":"2024-08-17T14:41:33.586507Z","shell.execute_reply":"2024-08-17T14:41:33.60066Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Calibrating Single training data\n### 前処理の方針\nデータの処理については、  \n- ピクセル応答の非線形性を補正する\n- フラットフィールド、ダークマップ、デッドピクセルマップを適用する\n- データをそのまま変更せずに残すか  \nを選択できる。  \n\n観測データは、AIRSでは30フレーム、FGS1では360フレームごとに時間的にビンニングされ、より軽量で使いやすいデータキューブが得られる。画像は、波長軸に沿ってピクセル39から321の間でカットされ、波長次元に残った282ピクセルが、AIRSからの最後の282ターゲットポイントに一致するように調整される。283番目のターゲットポイントはFGS1のもので、後で追加される。","metadata":{}},{"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-08-17T14:41:33.603302Z","iopub.execute_input":"2024-08-17T14:41:33.603657Z","iopub.status.idle":"2024-08-17T14:41:33.614441Z","shell.execute_reply.started":"2024-08-17T14:41:33.603626Z","shell.execute_reply":"2024-08-17T14:41:33.613154Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"files = glob.glob(os.path.join(path_folder + 'train/', '*/*'))\n# print(files)\n\nindex = get_index(files[:24],CHUNKS_SIZE)  ## 48 is hardcoded here but please feel free to remove it if you want to do it for the entire dataset\n","metadata":{"execution":{"iopub.status.busy":"2024-08-17T14:41:33.615856Z","iopub.execute_input":"2024-08-17T14:41:33.616192Z","iopub.status.idle":"2024-08-17T14:41:36.077472Z","shell.execute_reply.started":"2024-08-17T14:41:33.616153Z","shell.execute_reply":"2024-08-17T14:41:36.076089Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DO_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\ntrain_adc_info = pl.read_csv(os.path.join(path_folder, 'train_adc_info.csv'))\naxis_info = pl.read_parquet(os.path.join(path_folder,'axis_info.parquet'))\n\nAIRS_CH0_clean = np.zeros((CHUNKS_SIZE, 11250, 32, l))\nFGS1_clean = np.zeros((CHUNKS_SIZE, 135000, 32, 32))\n\ni = 0","metadata":{"execution":{"iopub.status.busy":"2024-08-17T14:41:36.081523Z","iopub.execute_input":"2024-08-17T14:41:36.081882Z","iopub.status.idle":"2024-08-17T14:41:36.332884Z","shell.execute_reply.started":"2024-08-17T14:41:36.081851Z","shell.execute_reply":"2024-08-17T14:41:36.331833Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### AIRS-CH0_signal","metadata":{}},{"cell_type":"code","source":"planet_id = index[0][0].item()\n\n\ndf = pl.read_parquet(os.path.join(path_folder,f'train/{planet_id}/AIRS-CH0_signal.parquet'))\nsignal = df.to_numpy().reshape((df.shape[0], 32, 356))\nplt.title('Natural image')\nplt.imshow(signal[0,:,:])\nplt.show()\n\ngain = train_adc_info.filter(pl.col('planet_id') == planet_id).select(['AIRS-CH0_adc_gain']).to_numpy()\noffset = train_adc_info.filter(pl.col('planet_id') == planet_id).select(['AIRS-CH0_adc_offset']).to_numpy()\n\n# AD Convert を復元している。\nsignal = ADC_convert(signal, gain, offset)\nplt.title('After ADC_convert')\nplt.imshow(signal[0,:,:])\nplt.show()\n\n# 積分時間を取ってきている。\ndt_airs = axis_info['AIRS-CH0-integration_time'].drop_nulls().to_numpy()\nchopped_signal = signal[:, :, cut_inf:cut_sup]\nplt.title('Chopped Signal')\nplt.imshow(chopped_signal[0,:,:])\nplt.show()\ndel signal, df","metadata":{"execution":{"iopub.status.busy":"2024-08-17T14:41:36.334388Z","iopub.execute_input":"2024-08-17T14:41:36.334808Z","iopub.status.idle":"2024-08-17T14:41:39.192529Z","shell.execute_reply.started":"2024-08-17T14:41:36.334767Z","shell.execute_reply":"2024-08-17T14:41:39.191259Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# CLEANING THE DATA: AIRS\nflat = pl.read_parquet(os.path.join(path_folder,f'train/{planet_id}/AIRS-CH0_calibration/flat.parquet')).to_numpy().reshape((32, 356))[:, cut_inf:cut_sup]\ndark = pl.read_parquet(os.path.join(path_folder,f'train/{planet_id}/AIRS-CH0_calibration/dark.parquet')).to_numpy().reshape((32, 356))[:, cut_inf:cut_sup]\ndead_airs = pl.read_parquet(os.path.join(path_folder,f'train/{planet_id}/AIRS-CH0_calibration/dead.parquet')).to_numpy().reshape((32, 356))[:, cut_inf:cut_sup]\nlinear_corr = pl.read_parquet(os.path.join(path_folder,f'train/{planet_id}/AIRS-CH0_calibration/linear_corr.parquet')).to_numpy().reshape((6, 32, 356))[:, :, cut_inf:cut_sup]\n","metadata":{"execution":{"iopub.status.busy":"2024-08-17T14:41:39.194108Z","iopub.execute_input":"2024-08-17T14:41:39.194592Z","iopub.status.idle":"2024-08-17T14:41:39.302886Z","shell.execute_reply.started":"2024-08-17T14:41:39.194549Z","shell.execute_reply":"2024-08-17T14:41:39.301677Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if DO_MASK:\n    chopped_signal = mask_hot_dead(chopped_signal, dead_airs, dark)\n    plt.title('Masked Chopped Signal')\n    plt.imshow(chopped_signal[0,:,:])\n    plt.show()\n    AIRS_CH0_clean[i] = chopped_signal\nelse:\n    AIRS_CH0_clean[i] = chopped_signal\n    \nif DO_THE_NL_CORR: \n    linear_corr_signal = apply_linear_corr(linear_corr,AIRS_CH0_clean[i])\n    plt.title('Linear corr chopped Signal')\n    plt.imshow(linear_corr_signal[0,:,:])\n    plt.show()\n    AIRS_CH0_clean[i,:, :, :] = linear_corr_signal\ndel linear_corr\n\nif DO_DARK: \n    cleaned_signal = clean_dark(AIRS_CH0_clean[i], dead_airs, dark, dt_airs)\n    AIRS_CH0_clean[i] = cleaned_signal\n    plt.title('Darked chopped Signal')\n    plt.imshow(cleaned_signal[0,:,:])\n    plt.show()\nelse: \n    pass\ndel dark","metadata":{"execution":{"iopub.status.busy":"2024-08-17T14:41:39.304026Z","iopub.execute_input":"2024-08-17T14:41:39.304415Z","iopub.status.idle":"2024-08-17T14:41:44.379216Z","shell.execute_reply.started":"2024-08-17T14:41:39.304383Z","shell.execute_reply":"2024-08-17T14:41:44.378098Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### FGS1_signal","metadata":{}},{"cell_type":"code","source":"df = pl.read_parquet(os.path.join(path_folder,f'train/{planet_id}/FGS1_signal.parquet'))\nfgs_signal = df.to_numpy().reshape((df.shape[0], 32, 32))\nplt.title('Natural image')\nplt.imshow(fgs_signal[0,:,:])\nplt.show()\n\nFGS1_gain = train_adc_info.filter(pl.col('planet_id') == planet_id).select(['FGS1_adc_gain']).to_numpy()\nFGS1_offset = train_adc_info.filter(pl.col('planet_id') == planet_id).select(['FGS1_adc_offset']).to_numpy()\nfgs_signal = ADC_convert(fgs_signal, FGS1_gain, FGS1_offset)\nplt.title('After ADC_convert')\nplt.imshow(fgs_signal[0,:,:])\nplt.show()\n\ndt_fgs1 = np.ones(len(fgs_signal))*0.1\nchopped_FGS1 = fgs_signal\nplt.title('Chopped Signal')\nplt.imshow(chopped_FGS1[0,:,:])\nplt.show()\ndel fgs_signal, df","metadata":{"execution":{"iopub.status.busy":"2024-08-17T14:41:44.381233Z","iopub.execute_input":"2024-08-17T14:41:44.381684Z","iopub.status.idle":"2024-08-17T14:41:46.294228Z","shell.execute_reply.started":"2024-08-17T14:41:44.381642Z","shell.execute_reply":"2024-08-17T14:41:46.293066Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# CLEANING THE DATA: FGS1\nflat = pl.read_parquet(os.path.join(path_folder,f'train/{planet_id}/FGS1_calibration/flat.parquet')).to_numpy().reshape((32, 32))\ndark = pl.read_parquet(os.path.join(path_folder,f'train/{planet_id}/FGS1_calibration/dark.parquet')).to_numpy().reshape((32, 32))\ndead_fgs1 = pl.read_parquet(os.path.join(path_folder,f'train/{planet_id}/FGS1_calibration/dead.parquet')).to_numpy().reshape((32, 32))\nlinear_corr = pl.read_parquet(os.path.join(path_folder,f'train/{planet_id}/FGS1_calibration/linear_corr.parquet')).to_numpy().reshape((6, 32, 32))","metadata":{"execution":{"iopub.status.busy":"2024-08-17T14:41:46.295909Z","iopub.execute_input":"2024-08-17T14:41:46.296386Z","iopub.status.idle":"2024-08-17T14:41:46.347741Z","shell.execute_reply.started":"2024-08-17T14:41:46.296343Z","shell.execute_reply":"2024-08-17T14:41:46.346507Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if DO_MASK:\n    chopped_FGS1 = mask_hot_dead(chopped_FGS1, dead_fgs1, dark)\n    plt.title('Masked Chopped Signal')\n    plt.imshow(chopped_FGS1[0,:,:])\n    plt.show()\n    FGS1_clean[i] = chopped_FGS1\nelse:\n    FGS1_clean[i] = chopped_FGS1\n\nif DO_THE_NL_CORR: \n    linear_corr_signal = apply_linear_corr(linear_corr,FGS1_clean[i])\n    plt.title('Linear corr chopped Signal')\n    plt.imshow(linear_corr_signal[0,:,:])\n    plt.show()\n    FGS1_clean[i,:, :, :] = linear_corr_signal\ndel linear_corr\n\nif DO_DARK: \n    cleaned_signal = clean_dark(FGS1_clean[i], dead_fgs1, dark,dt_fgs1)\n    plt.title('Darked chopped Signal')\n    plt.imshow(cleaned_signal[0,:,:])\n    plt.show()\n    FGS1_clean[i] = cleaned_signal\nelse: \n    pass\ndel dark\n\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2024-08-17T14:41:46.349369Z","iopub.execute_input":"2024-08-17T14:41:46.34973Z","iopub.status.idle":"2024-08-17T14:41:53.186631Z","shell.execute_reply.started":"2024-08-17T14:41:46.349696Z","shell.execute_reply":"2024-08-17T14:41:53.185216Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for 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    for i in range (CHUNKS_SIZE) : \n        df = pl.read_parquet(os.path.join(path_folder,f'train/{index_chunk[i]}/AIRS-CH0_signal.parquet'))\n        signal = df.to_numpy().reshape((df.shape[0], 32, 356))\n        gain = train_adc_info.filter(pl.col('planet_id') == index_chunk[i]).select(['AIRS-CH0_adc_gain']).to_numpy()\n        offset = train_adc_info.filter(pl.col('planet_id') == index_chunk[i]).select(['AIRS-CH0_adc_offset']).to_numpy()\n        signal = ADC_convert(signal, gain, offset)\n        dt_airs = axis_info['AIRS-CH0-integration_time'].drop_nulls().to_numpy()\n        chopped_signal = signal[:, :, cut_inf:cut_sup]\n        del signal, df\n        \n        # CLEANING THE DATA: AIRS\n        flat = pl.read_parquet(os.path.join(path_folder,f'train/{index_chunk[i]}/AIRS-CH0_calibration/flat.parquet')).to_numpy().reshape((32, 356))[:, cut_inf:cut_sup]\n        dark = pl.read_parquet(os.path.join(path_folder,f'train/{index_chunk[i]}/AIRS-CH0_calibration/dark.parquet')).to_numpy().reshape((32, 356))[:, cut_inf:cut_sup]\n        dead_airs = pl.read_parquet(os.path.join(path_folder,f'train/{index_chunk[i]}/AIRS-CH0_calibration/dead.parquet')).to_numpy().reshape((32, 356))[:, cut_inf:cut_sup]\n        linear_corr = pl.read_parquet(os.path.join(path_folder,f'train/{index_chunk[i]}/AIRS-CH0_calibration/linear_corr.parquet')).to_numpy().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_THE_NL_CORR: \n            linear_corr_signal = apply_linear_corr(linear_corr,AIRS_CH0_clean[i])\n            AIRS_CH0_clean[i,:, :, :] = linear_corr_signal\n        del linear_corr\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        df = pl.read_parquet(os.path.join(path_folder,f'train/{index_chunk[i]}/FGS1_signal.parquet'))\n        fgs_signal = df.to_numpy().reshape((df.shape[0], 32, 32))\n        \n        FGS1_gain = train_adc_info.filter(pl.col('planet_id') == index_chunk[i]).select(['FGS1_adc_gain']).to_numpy()\n        FGS1_offset = train_adc_info.filter(pl.col('planet_id') == index_chunk[i]).select(['FGS1_adc_offset']).to_numpy()\n        \n        fgs_signal = ADC_convert(fgs_signal, FGS1_gain, FGS1_offset)\n        dt_fgs1 = np.ones(len(fgs_signal))*0.1\n        chopped_FGS1 = fgs_signal\n        del fgs_signal, df\n        \n        # CLEANING THE DATA: FGS1\n        flat = pl.read_parquet(os.path.join(path_folder,f'train/{index_chunk[i]}/FGS1_calibration/flat.parquet')).to_numpy().reshape((32, 32))\n        dark = pl.read_parquet(os.path.join(path_folder,f'train/{index_chunk[i]}/FGS1_calibration/dark.parquet')).to_numpy().reshape((32, 32))\n        dead_fgs1 = pl.read_parquet(os.path.join(path_folder,f'train/{index_chunk[i]}/FGS1_calibration/dead.parquet')).to_numpy().reshape((32, 32))\n        linear_corr = pl.read_parquet(os.path.join(path_folder,f'train/{index_chunk[i]}/FGS1_calibration/linear_corr.parquet')).to_numpy().reshape((6, 32, 32))\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    # 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    ## (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    for i in range (CHUNKS_SIZE):\n        flat_airs = pl.read_parquet(os.path.join(path_folder,f'train/{index_chunk[i]}/AIRS-CH0_calibration/flat.parquet')).to_numpy().reshape((32, 356))[:, cut_inf:cut_sup]\n        flat_fgs = pl.read_parquet(os.path.join(path_folder,f'train/{index_chunk[i]}/FGS1_calibration/flat.parquet')).to_numpy().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    # print(gc.collect())\n    \n    ## save data\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-08-17T14:41:53.188511Z","iopub.execute_input":"2024-08-17T14:41:53.188916Z","iopub.status.idle":"2024-08-17T14:43:17.004664Z","shell.execute_reply.started":"2024-08-17T14:41:53.188882Z","shell.execute_reply":"2024-08-17T14:43:17.003495Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}