{"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":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# NeurIPS DataVisalize [Japanese]\n\nデータの見方を理解していくためのNotebook.\n\nコンペのURL\nhttps://www.kaggle.com/competitions/ariel-data-challenge-2024/overview","metadata":{}},{"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nimport itertools\nfrom matplotlib import pyplot as plt\nimport seaborn as sns\nfrom astropy.stats import sigma_clip","metadata":{"execution":{"iopub.status.busy":"2024-08-24T13:33:33.155509Z","iopub.execute_input":"2024-08-24T13:33:33.156554Z","iopub.status.idle":"2024-08-24T13:33:33.162822Z","shell.execute_reply.started":"2024-08-24T13:33:33.156502Z","shell.execute_reply":"2024-08-24T13:33:33.1615Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#ディレクトリ\nINPUT_DIR = \"/kaggle/input/ariel-data-challenge-2024\"\nTRAIN_DIR = f\"{INPUT_DIR}/train\"","metadata":{"execution":{"iopub.status.busy":"2024-08-24T14:35:42.596507Z","iopub.execute_input":"2024-08-24T14:35:42.597061Z","iopub.status.idle":"2024-08-24T14:35:42.60368Z","shell.execute_reply.started":"2024-08-24T14:35:42.597012Z","shell.execute_reply":"2024-08-24T14:35:42.602185Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#確認するデータを指定\n\nplanet_id = 100468857","metadata":{"execution":{"iopub.status.busy":"2024-08-24T14:35:26.685929Z","iopub.execute_input":"2024-08-24T14:35:26.687614Z","iopub.status.idle":"2024-08-24T14:35:26.694096Z","shell.execute_reply.started":"2024-08-24T14:35:26.687551Z","shell.execute_reply":"2024-08-24T14:35:26.692527Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"AIRS-CH0とFGS1の2つの機器のデータがある。  \n2次元画像データに直すにはreshapeする必要がある。\n- AIRS-CH0 : 1.95 ～ 3.90 µm の感度を持つ赤外線分光計\n- FGS1 : 0.6 ~ 0.8 um(可視光)のカメラ. 星が写っているはず。\n","metadata":{}},{"cell_type":"code","source":"CH0_df = pd.read_parquet(f\"{TRAIN_DIR}/{planet_id}/AIRS-CH0_signal.parquet\")\nFGS1_df = pd.read_parquet(f\"{TRAIN_DIR}/{planet_id}/FGS1_signal.parquet\")\nCH0_signal = CH0_df.to_numpy().reshape(CH0_df.shape[0], 32, 356)   #32×356の画像\nFGS1_signal = FGS1_df.to_numpy().reshape(FGS1_df.shape[0], 32, 32) #32×32の画像\nprint(f\"AIRS-CH0 : {CH0_signal.shape}\")\nprint(f\"FGS1 : {FGS1_signal.shape}\")","metadata":{"execution":{"iopub.status.busy":"2024-08-24T13:33:35.129862Z","iopub.execute_input":"2024-08-24T13:33:35.131083Z","iopub.status.idle":"2024-08-24T13:33:37.935915Z","shell.execute_reply.started":"2024-08-24T13:33:35.13099Z","shell.execute_reply":"2024-08-24T13:33:37.934731Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"最初の画像は以下の通り。","metadata":{}},{"cell_type":"code","source":"fig,ax = plt.subplots(1, 2, figsize = (10, 4))\nsns.heatmap(CH0_signal[0, :, :], cmap=\"rocket\", ax=ax[0])\nsns.heatmap(FGS1_signal[0, :, :], cmap=\"rocket\", ax=ax[1])\nax[0].set_title(f\"planet_id = {planet_id} : AIRS-CH0\", size=12)\nax[1].set_title(f\"planet_id = {planet_id} : FGS1\", size=12)","metadata":{"execution":{"iopub.status.busy":"2024-08-24T13:33:37.93832Z","iopub.execute_input":"2024-08-24T13:33:37.938696Z","iopub.status.idle":"2024-08-24T13:33:39.315148Z","shell.execute_reply.started":"2024-08-24T13:33:37.938654Z","shell.execute_reply":"2024-08-24T13:33:39.313881Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"AIRS-CH0, FGS1はそれぞれダークフレームなどの較正用のデータを持つ。","metadata":{}},{"cell_type":"code","source":"#較正データの読み込み\nCH0_dark = pd.read_parquet(f\"{TRAIN_DIR}/{planet_id}/AIRS-CH0_calibration/dark.parquet\").to_numpy().reshape(1, 32, 356)\nCH0_dead = pd.read_parquet(f\"{TRAIN_DIR}/{planet_id}/AIRS-CH0_calibration/dead.parquet\").to_numpy().reshape(1, 32, 356)\nCH0_flat = pd.read_parquet(f\"{TRAIN_DIR}/{planet_id}/AIRS-CH0_calibration/flat.parquet\").to_numpy().reshape(1, 32, 356)\nCH0_linear_corr = pd.read_parquet(f\"{TRAIN_DIR}/{planet_id}/AIRS-CH0_calibration/linear_corr.parquet\").to_numpy().reshape(1, 32, -1) #2136 = 356 * 6\nCH0_read = pd.read_parquet(f\"{TRAIN_DIR}/{planet_id}/AIRS-CH0_calibration/read.parquet\").to_numpy().reshape(1, 32, 356)\n\nFGS1_dark = pd.read_parquet(f\"{TRAIN_DIR}/{planet_id}/FGS1_calibration/dark.parquet\").to_numpy().reshape(1, 32, 32)\nFGS1_dead = pd.read_parquet(f\"{TRAIN_DIR}/{planet_id}/FGS1_calibration/dead.parquet\").to_numpy().reshape(1, 32, 32)\nFGS1_flat = pd.read_parquet(f\"{TRAIN_DIR}/{planet_id}/FGS1_calibration/flat.parquet\").to_numpy().reshape(1, 32, 32)\nFGS1_linear_corr = pd.read_parquet(f\"{TRAIN_DIR}/{planet_id}/FGS1_calibration/linear_corr.parquet\").to_numpy().reshape(1, 32, 192) # 192 = 32*6\nFGS1_read = pd.read_parquet(f\"{TRAIN_DIR}/{planet_id}/FGS1_calibration/read.parquet\").to_numpy().reshape(1, 32, 32)\n\ncalib_data_list = [CH0_dark, CH0_dead, CH0_flat, CH0_linear_corr, CH0_read, \n                   FGS1_dark, FGS1_dead, FGS1_flat, FGS1_linear_corr, FGS1_read]\ncalib_dataname_list = [\"AIRS-CH0 : dark\", \"AIRS-CH0 : dead\", \"AIRS-CH0 : flat\", \"AIRS-CH0 : linear_corr\", \"AIRS-CH0 : read\",\n                       \"FGS1 : dark\", \"FGS1 : dead\", \"FGS1 : flat\", \"FGS1 : linear_corr\", \"FGS1 : read\"]","metadata":{"execution":{"iopub.status.busy":"2024-08-24T13:37:25.02054Z","iopub.execute_input":"2024-08-24T13:37:25.021943Z","iopub.status.idle":"2024-08-24T13:37:25.216143Z","shell.execute_reply.started":"2024-08-24T13:37:25.021885Z","shell.execute_reply":"2024-08-24T13:37:25.214691Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig,ax = plt.subplots(2, 5, figsize = (16, 4))\nax = ax.ravel()\n\nfor i in range(len(ax)):\n    data = calib_data_list[i]\n    data_name = calib_dataname_list[i]\n    sns.heatmap(data[0, :, :], cmap=\"rocket\", ax=ax[i])\n    ax[i].set_title(f\"{data_name}\", size=8)\n    ax[i].set_xticks([])\n    ax[i].set_yticks([])\n    ax[i].set_xticklabels([])","metadata":{"execution":{"iopub.status.busy":"2024-08-24T13:37:25.32486Z","iopub.execute_input":"2024-08-24T13:37:25.326207Z","iopub.status.idle":"2024-08-24T13:37:30.838678Z","shell.execute_reply.started":"2024-08-24T13:37:25.326148Z","shell.execute_reply":"2024-08-24T13:37:30.837365Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Calibration\n以下のNotebookを参考に、画像処理の部分について理解していきたい。  \nhttps://www.kaggle.com/code/gordonyip/update-calibrating-and-binning-astronomical-data","metadata":{}},{"cell_type":"markdown","source":"## Step1. Analog-to-Digital Conversion\nデータをgainで除算して、offsetを追加する必要がある.  \ngain, offsetはtrain_adc_info.csvに記載されている。","metadata":{}},{"cell_type":"code","source":"'''\nplanet_id, deviceに対応したgain, offsetを読み込む関数\n'''\ndef get_gain(planet_id, device):\n    train_adc_info_df = pd.read_csv(f\"{INPUT_DIR}/train_adc_info.csv\")\n    adc_planet = train_adc_info_df.loc[train_adc_info_df[\"planet_id\"]==int(planet_id)]\n    gain = adc_planet[f\"{device}_adc_gain\"].item()\n    \n    return gain\n\n\ndef get_offset(planet_id, device):\n    train_adc_info_df = pd.read_csv(f\"{INPUT_DIR}/train_adc_info.csv\")\n    adc_planet = train_adc_info_df.loc[train_adc_info_df[\"planet_id\"]==int(planet_id)]\n    offset = adc_planet[f\"{device}_adc_offset\"].item()\n    \n    return offset","metadata":{"execution":{"iopub.status.busy":"2024-08-24T13:34:15.421995Z","iopub.execute_input":"2024-08-24T13:34:15.423396Z","iopub.status.idle":"2024-08-24T13:34:15.431611Z","shell.execute_reply.started":"2024-08-24T13:34:15.42334Z","shell.execute_reply":"2024-08-24T13:34:15.430324Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def ADC_convert(signal, gain, offset):\n    signal = signal.astype(np.float64)\n    signal = signal / gain\n    signal = signal + offset\n    return signal","metadata":{"execution":{"iopub.status.busy":"2024-08-24T13:34:15.710868Z","iopub.execute_input":"2024-08-24T13:34:15.712426Z","iopub.status.idle":"2024-08-24T13:34:15.718666Z","shell.execute_reply.started":"2024-08-24T13:34:15.712362Z","shell.execute_reply":"2024-08-24T13:34:15.717233Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Step 2: Mask hot/dead pixel\n- hot pixel : 正しく機能していないピクセル\n- dead pixel : 光に全く反応しないピクセル","metadata":{}},{"cell_type":"code","source":"def mask_hot_dead(signal, dead, dark):\n    hot = sigma_clip(dark, sigma=5, maxiters=5).mask #ダークフレームのデータのうち、外れ値が入っている位置がhotピクセル\n    \n    hot = np.tile(hot, (signal.shape[0], 1, 1))\n    dead = np.tile(dead, (signal.shape[0], 1, 1))\n    \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-24T13:34:16.934663Z","iopub.execute_input":"2024-08-24T13:34:16.935138Z","iopub.status.idle":"2024-08-24T13:34:16.94382Z","shell.execute_reply.started":"2024-08-24T13:34:16.935094Z","shell.execute_reply":"2024-08-24T13:34:16.942248Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Step 3: linearity Correction\n[リンク先](https://www.jstage.jst.go.jp/article/itej1997/57/2/57_2_203/_pdf)で言及されている非線形現象の補正?  \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(range(clean_signal.shape[1]), range(clean_signal.shape[2])):\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-24T14:12:54.291647Z","iopub.execute_input":"2024-08-24T14:12:54.293336Z","iopub.status.idle":"2024-08-24T14:12:54.30164Z","shell.execute_reply.started":"2024-08-24T14:12:54.293271Z","shell.execute_reply":"2024-08-24T14:12:54.30018Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Step 4: dark current subtraction\nダークフレームの減算処理。  \nAIRS-CH0のdtはaxis_info.parquetに記録されている。  \nFGS1のdtは0.1s固定。","metadata":{}},{"cell_type":"code","source":"def get_dt(planet_id, device):\n    if device == \"AIRS-CH0\":\n        axis_info = pd.read_parquet(f\"{INPUT_DIR}/axis_info.parquet\")\n        dt = axis_info['AIRS-CH0-integration_time'].dropna().values\n    else:\n        dt = np.ones(135000)*0.1\n    return dt","metadata":{"execution":{"iopub.status.busy":"2024-08-24T13:34:17.880286Z","iopub.execute_input":"2024-08-24T13:34:17.880796Z","iopub.status.idle":"2024-08-24T13:34:17.887607Z","shell.execute_reply.started":"2024-08-24T13:34:17.880753Z","shell.execute_reply":"2024-08-24T13:34:17.88625Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def clean_dark(signal, dead, dark, dt):\n\n    dark = np.ma.masked_where(dead, dark) #deadピクセルはマスクする\n    dark = np.tile(dark, (signal.shape[0], 1, 1))\n\n    signal = signal - (dark * dt[:, np.newaxis, np.newaxis])\n    return signal","metadata":{"execution":{"iopub.status.busy":"2024-08-24T13:34:18.161974Z","iopub.execute_input":"2024-08-24T13:34:18.16296Z","iopub.status.idle":"2024-08-24T13:34:18.16969Z","shell.execute_reply.started":"2024-08-24T13:34:18.16291Z","shell.execute_reply":"2024-08-24T13:34:18.168406Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Step 5: Get Correlated Double Sampling (CDS)\n理解できてない。","metadata":{}},{"cell_type":"code","source":"def get_cds(signal):\n    cds = signal[:, 1::2, :, :] - signal[:, ::2, :, :]\n    return cds\n","metadata":{"execution":{"iopub.status.busy":"2024-08-24T13:34:18.703961Z","iopub.execute_input":"2024-08-24T13:34:18.705132Z","iopub.status.idle":"2024-08-24T13:34:18.710953Z","shell.execute_reply.started":"2024-08-24T13:34:18.705081Z","shell.execute_reply":"2024-08-24T13:34:18.709588Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Step 6 (Optional): Time Binning\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\n","metadata":{"execution":{"iopub.status.busy":"2024-08-24T13:34:19.26284Z","iopub.execute_input":"2024-08-24T13:34:19.263871Z","iopub.status.idle":"2024-08-24T13:34:19.27205Z","shell.execute_reply.started":"2024-08-24T13:34:19.263802Z","shell.execute_reply":"2024-08-24T13:34:19.270562Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Step 7: Flat Field Correction","metadata":{}},{"cell_type":"code","source":"def correct_flat_field(flat, dead, signal):\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-24T13:37:56.757431Z","iopub.execute_input":"2024-08-24T13:37:56.758453Z","iopub.status.idle":"2024-08-24T13:37:56.76538Z","shell.execute_reply.started":"2024-08-24T13:37:56.758398Z","shell.execute_reply":"2024-08-24T13:37:56.764086Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n元の画像から順番に処理していく。  \nlinear, CDS, Time Binningはよく理解できていないのでskip.","metadata":{}},{"cell_type":"code","source":"gain = get_gain(planet_id, \"AIRS-CH0\")\noffset  = get_offset(planet_id, \"AIRS-CH0\")\n\nCH0_dt = get_dt(planet_id, \"AIRS-CH0\")\nCH0_signal_ADC = ADC_convert(CH0_signal, gain, offset)\nCH0_signal_mask = mask_hot_dead(CH0_signal_ADC, CH0_dead, CH0_dark)\nCH0_signal_dark = clean_dark(CH0_signal_mask, CH0_dead, CH0_dark, CH0_dt)\nCH0_signal_flat = correct_flat_field(CH0_flat, CH0_dead, CH0_signal_dark)\n\nFGS1_dt = get_dt(planet_id, \"FGS1\")\nFGS1_signal_ADC = ADC_convert(FGS1_signal, gain, offset)\nFGS1_signal_mask = mask_hot_dead(FGS1_signal_ADC, FGS1_dead, FGS1_dark)\nFGS1_signal_dark = clean_dark(FGS1_signal_mask, FGS1_dead, FGS1_dark, FGS1_dt)\nFGS1_signal_flat = correct_flat_field(FGS1_flat, FGS1_dead, FGS1_signal_dark)\n","metadata":{"execution":{"iopub.status.busy":"2024-08-24T14:25:35.347377Z","iopub.execute_input":"2024-08-24T14:25:35.348859Z","iopub.status.idle":"2024-08-24T14:25:53.457803Z","shell.execute_reply.started":"2024-08-24T14:25:35.348792Z","shell.execute_reply":"2024-08-24T14:25:53.456781Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data_list = [CH0_signal, CH0_signal_ADC, CH0_signal_dark, CH0_signal_flat,\n             FGS1_signal, FGS1_signal_ADC, FGS1_signal_dark, FGS1_signal_flat]\ndata_name_list = [\"row\", \"ADC\", \"dark\", \"flat\", \n                  \"row\", \"ADC\", \"dark\", \"flat\"]\n\nfig,ax = plt.subplots(2, 4, figsize = (16, 6))\nax = ax.ravel()\n\nfor i in range(len(ax)):\n    data = data_list[i]\n    data_name = data_name_list[i]\n    sns.heatmap(data[0, :, :], cmap=\"rocket\", ax=ax[i])\n    ax[i].set_title(f\"{data_name}\", size=8)\n    ax[i].set_xticks([])\n    ax[i].set_yticks([])\n    ax[i].set_xticklabels([])","metadata":{"execution":{"iopub.status.busy":"2024-08-24T14:27:06.011757Z","iopub.execute_input":"2024-08-24T14:27:06.012204Z","iopub.status.idle":"2024-08-24T14:27:10.304269Z","shell.execute_reply.started":"2024-08-24T14:27:06.012163Z","shell.execute_reply":"2024-08-24T14:27:10.302917Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}