{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"### Quick Note\nI originally wanted to submit this while the competition was running, \nbut only joined in the last week and got the warning that notebooks may not be made public this late.\nI hope it still provides some interesting information.\n\nMy own submissions used the recommended band08 in addition to Ash, but I did not have time left to verify if there was a substantial advantage against Ash only, I'm somewhat skeptical of that now, having not seen it used by any of the top solutions posted yet.","metadata":{}},{"cell_type":"markdown","source":"# Summary and Key Takeaway\nTaking the Ash false color image as a basis seems great, but I expect most people are doing that already.\nJust adding bands 11 to 16 to that base seems like a bad idea, as the information is almost completely contained in blue already.\nI would recommend trying to add band08 as a 4th input channel (normalized), as it is much less correlated with the Ash colors and pretty well correlated with the contrails, meaning it could add new and useful information, an implementation is here [Contrails Dataset (Ash+08)](https://www.kaggle.com/code/raki21/contrails-dataset-ash-band08-with-soft-mask). I would not add all 3 (bands 08,09 and 10) as they are highly correlated with each other and additional channels that don't contribute additional information might increase opportunity for overfitting. ","metadata":{}},{"cell_type":"markdown","source":"# Information Bands\nA lot of information about bands is already contained in other EDAs, especially [this EDA](https://www.kaggle.com/code/pranavnadimpali/comprehensive-eda-submission), a short summary:\n\nThe Advanced Baseline Imager (ABI) is a key instrument on the Geostationary Operational Environmental Satellite (GOES) series. It captures images of the Earth using 16 spectral bands, with each band focusing on a specific wavelength. This provides a wealth of information about the Earth's atmosphere, clouds, land, and water, significantly improving weather analysis and forecasting. In this dataset, you have access to 9 bands for each example.\n\nEach band provides a series of images taken at 10-minute intervals, leading to 8 images for each band spanning 80 minutes. This temporal data captures how contrails evolve over time. Two types of segmentation masks are provided: 'human_pixel_masks' that represent the consolidated ground truth, and 'human_individual_masks' that represent annotations from multiple labellers. The ground truth corresponds to the 5th image in the bands.\n\nThe image that labellers annotate is a false color image, which is not directly provided in the spectral bands but can be generated from them. The false color image uses the [Ash color scheme](https://eumetrain.org/sites/default/files/2020-05/RGB_recipes.pdf), using bands 15,14 and 11, making contrails appear darker and thus easier to detect. The ash color scheme combines red, green, and blue channels to represent different features. \n","metadata":{}},{"cell_type":"markdown","source":"## Correlation Between Bands and Contrails\nThe contribution of this notebook is trying to find bands that are useful for contrail detection.\nFor this we create a correlation matrix between the bands and the contrail groundtruth.","metadata":{}},{"cell_type":"code","source":"import os\nimport pandas as pd\nimport numpy as np\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nfrom pathlib import Path\nfrom tqdm.notebook import tqdm","metadata":{"execution":{"iopub.status.busy":"2023-08-21T11:47:17.339867Z","iopub.execute_input":"2023-08-21T11:47:17.340218Z","iopub.status.idle":"2023-08-21T11:47:17.962552Z","shell.execute_reply.started":"2023-08-21T11:47:17.340189Z","shell.execute_reply":"2023-08-21T11:47:17.961192Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data_dir = '/kaggle/input/google-research-identify-contrails-reduce-global-warming/'\n\ntrain_rs = os.listdir(data_dir + 'train')\nvalid_rs = os.listdir(data_dir + 'validation')\n\ntrain_df = pd.DataFrame(train_rs, columns=['record_id'])\nvalid_df = pd.DataFrame(valid_rs, columns=['record_id'])\n\ntrain_df['train'] = 'train'\nvalid_df['train'] = 'valid'","metadata":{"execution":{"iopub.status.busy":"2023-08-21T11:47:17.964774Z","iopub.execute_input":"2023-08-21T11:47:17.965197Z","iopub.status.idle":"2023-08-21T11:47:18.496822Z","shell.execute_reply.started":"2023-08-21T11:47:17.965167Z","shell.execute_reply":"2023-08-21T11:47:18.495538Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_correlation_matrix(correlation_matrix, band_labels, band_labels_nickname):\n    plt.figure(figsize=(10,10)) # Set the figure size\n    sns.set(font_scale=1.2)  # Set font size\n    \n    # Generate a mask for the upper triangle\n    mask = np.triu(np.ones_like(correlation_matrix, dtype=bool))\n    \n    # Generate a custom diverging colormap\n    cmap = sns.diverging_palette(230, 20, as_cmap=True)\n    \n    # Draw the heatmap with the mask and correct aspect ratio\n    sns.set(font_scale=0.7)\n    sns.heatmap(correlation_matrix, annot=True, fmt=\".3f\", mask=mask, cmap=cmap, \n                cbar_kws={\"shrink\": .5}, xticklabels=band_labels, yticklabels=band_labels_nickname)\n    plt.xlabel('Band')  # Adding label to x-axis\n    plt.ylabel('Band Nickname')  # Adding label to y-axis\n    plt.title('Correlation matrix of Band Data')\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-08-21T11:47:18.498109Z","iopub.execute_input":"2023-08-21T11:47:18.498414Z","iopub.status.idle":"2023-08-21T11:47:18.507111Z","shell.execute_reply.started":"2023-08-21T11:47:18.49839Z","shell.execute_reply":"2023-08-21T11:47:18.505421Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def correlation_create(record_id, directory):\n    record_data = {}\n    for band in range(8, 17):  # For bands 8 through 16\n        band_key = f\"band_{band:02d}\"\n        band_data = np.load(os.path.join(directory, record_id, band_key + \".npy\"))\n        record_data[band_key] = band_data[:,:,4]                                          # Only look at current time\n\n    if directory == 'train':\n        individual = np.load(os.path.join(directory, record_id, 'human_individual_masks' + \".npy\"))\n        record_data['mean_mask'] = individual.sum(axis=3) / individual.shape[3]\n    else:\n        record_data['mask'] = np.load(os.path.join(directory, record_id, 'human_pixel_masks' + \".npy\")).squeeze()\n\n    band_values = np.array([band_data.flatten() for band_data in record_data.values()])\n    correlation_matrix = np.corrcoef(band_values)\n    correlation_matrix = np.nan_to_num(correlation_matrix)\n\n    return correlation_matrix","metadata":{"execution":{"iopub.status.busy":"2023-08-21T11:47:18.508605Z","iopub.execute_input":"2023-08-21T11:47:18.50897Z","iopub.status.idle":"2023-08-21T11:47:18.520201Z","shell.execute_reply.started":"2023-08-21T11:47:18.508942Z","shell.execute_reply":"2023-08-21T11:47:18.518485Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for n,i in enumerate(train_rs):\n    if n == 0:\n        cor = correlation_create(str(i), data_dir+'train')\n        total_cor = cor\n    elif n%10 == 0:\n        cor = correlation_create(str(i), data_dir+'train')\n        total_cor += cor\n    if n == 10000 - 1:\n        break\ncorrelation_matrix = total_cor/1000","metadata":{"execution":{"iopub.status.busy":"2023-08-21T11:47:18.523867Z","iopub.execute_input":"2023-08-21T11:47:18.525496Z","iopub.status.idle":"2023-08-21T11:53:48.492594Z","shell.execute_reply.started":"2023-08-21T11:47:18.525432Z","shell.execute_reply":"2023-08-21T11:53:48.489269Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"band_labels = [\"band_08\", \"band_09\", \"band_10\", \"band_11\", \"band_12\", \"band_13\", \"band_14\", \"band_15\", \"band_16\",\"mask\"]\nband_labels_nickname = [\"Upper-Level Tropospheric Water Vapor\", \n                       \"Mid-Level Tropospheric Water Vapor\", \n                       \"Lower-level Water Vapor\", \n                       \"Cloud-Top Phase\", \n                       \"Ozone\", \n                       \"Clean IR Longwave Window\", \n                       \"IR Longwave Window\", \n                       \"Dirty Longwave Window\", \n                       \"CO2 Longwave Infrared\",\n                       \"Contrails\"]\nplot_correlation_matrix(correlation_matrix, band_labels, band_labels_nickname)","metadata":{"execution":{"iopub.status.busy":"2023-08-21T11:53:48.496263Z","iopub.execute_input":"2023-08-21T11:53:48.496775Z","iopub.status.idle":"2023-08-21T11:53:49.169889Z","shell.execute_reply.started":"2023-08-21T11:53:48.496727Z","shell.execute_reply":"2023-08-21T11:53:49.168773Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"You can see that bands 8,9,10 and bands 11 to 16 are highly correlated with each other, the bands 8 and 16 seem to be most correlated with contrails.\nI made a few small training runs with only using those bands leading to Val Dice that seem to align well with the correlation scores.","metadata":{}},{"cell_type":"markdown","source":"    Band  Dice  Val Dice\n    8 \t0.5  \t0.377\n    9 \t0.487\t0.372\n    10\t0.47 \t0.346\n    11\t0.349\t0.195\n    12\t0.323\t0.165\n    13\t0.342\t0.192\n    14\t0.454\t0.297\n    15\t0.484\t0.354\n    16\t0.518\t0.396","metadata":{}},{"cell_type":"markdown","source":"We now look at the correlation with the Ash color scheme.","metadata":{}},{"cell_type":"markdown","source":"## Correlation Ash","metadata":{}},{"cell_type":"code","source":"def normalize_range(data, bounds):\n    return (data - bounds[0]) / (bounds[1] - bounds[0])","metadata":{"execution":{"iopub.status.busy":"2023-08-21T11:53:49.171928Z","iopub.execute_input":"2023-08-21T11:53:49.172607Z","iopub.status.idle":"2023-08-21T11:53:49.181663Z","shell.execute_reply.started":"2023-08-21T11:53:49.172569Z","shell.execute_reply":"2023-08-21T11:53:49.179152Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def correlation_create_ash(record_id, directory):\n    _T11_BOUNDS = (243, 303)                   #\n    _CLOUD_TOP_TDIFF_BOUNDS = (-4, 5)\n    _TDIFF_BOUNDS = (-4, 2)\n\n    record_data = {}\n    for band in range(8, 17):  # For bands 8 through 16\n        band_key = f\"band_{band:02d}\"\n        band_data = np.load(os.path.join(directory, record_id, band_key + \".npy\"))\n        record_data[band_key] = band_data[:,:,4]                                          # Only look at current time\n\n    record_data['r'] = normalize_range(record_data[\"band_15\"] - record_data[\"band_14\"], _TDIFF_BOUNDS)\n    record_data['g'] = normalize_range(record_data[\"band_14\"] - record_data[\"band_11\"], _CLOUD_TOP_TDIFF_BOUNDS)\n    record_data['b'] = normalize_range(record_data[\"band_14\"], _T11_BOUNDS)\n    \n    if directory == 'train':\n        individual = np.load(os.path.join(directory, record_id, 'human_individual_masks' + \".npy\"))\n        record_data['mean_mask'] = individual.sum(axis=3) / individual.shape[3]\n    else:\n        record_data['mask'] = np.load(os.path.join(directory, record_id, 'human_pixel_masks' + \".npy\")).squeeze()\n\n    band_values = np.array([band_data.flatten() for band_data in record_data.values()])\n    correlation_matrix = np.corrcoef(band_values)\n    correlation_matrix = np.nan_to_num(correlation_matrix)\n\n    return correlation_matrix","metadata":{"execution":{"iopub.status.busy":"2023-08-21T11:53:49.18347Z","iopub.execute_input":"2023-08-21T11:53:49.183843Z","iopub.status.idle":"2023-08-21T11:53:49.201134Z","shell.execute_reply.started":"2023-08-21T11:53:49.183783Z","shell.execute_reply":"2023-08-21T11:53:49.198945Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for n,i in enumerate(train_rs):\n    if n == 0:\n        cor = correlation_create_ash(str(i), data_dir+'train')\n        total_cor = cor\n    elif n%10 == 0:\n        cor = correlation_create_ash(str(i), data_dir+'train')\n        total_cor += cor\n    if n == 10000 - 1:\n        break\ncorrelation_matrix = total_cor/1000","metadata":{"execution":{"iopub.status.busy":"2023-08-21T11:53:49.203366Z","iopub.execute_input":"2023-08-21T11:53:49.203828Z","iopub.status.idle":"2023-08-21T11:55:33.312291Z","shell.execute_reply.started":"2023-08-21T11:53:49.203777Z","shell.execute_reply":"2023-08-21T11:55:33.309907Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"band_labels = [\"band_08\", \"band_09\", \"band_10\", \"band_11\", \"band_12\", \"band_13\", \"band_14\", \"band_15\", \"band_16\", \"r\", \"g\", \"b\", \"mask\"]\nband_labels_nickname = [\"Upper-Level Tropospheric Water Vapor\", \n                       \"Mid-Level Tropospheric Water Vapor\", \n                       \"Lower-level Water Vapor\", \n                       \"Cloud-Top Phase\", \n                       \"Ozone\", \n                       \"Clean IR Longwave Window\", \n                       \"IR Longwave Window\", \n                       \"Dirty Longwave Window\", \n                       \"CO2 Longwave Infrared\",\n                       \"Red\", \"Green\", \"Blue\",\n                       \"Contrails\"]\nplot_correlation_matrix(correlation_matrix, band_labels, band_labels_nickname)","metadata":{"execution":{"iopub.status.busy":"2023-08-21T11:55:33.315214Z","iopub.execute_input":"2023-08-21T11:55:33.315654Z","iopub.status.idle":"2023-08-21T11:55:34.018572Z","shell.execute_reply.started":"2023-08-21T11:55:33.315618Z","shell.execute_reply":"2023-08-21T11:55:34.016Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The correlation between the red and green false colors in the ash image with contrails is very high and they are not highly correlated with any other bands. \nI except this to mean that they are really useful, which is plausible as this is what labelers saw. The exception is with Blue but this is also almost the same as band 14. ","metadata":{}},{"cell_type":"markdown","source":"Taking the Ash false color image as a basis seems great, adding bands 11 to 16 to that base seems like a bad idea, as the information is almost completely contained in blue already, I would recommend trying to add band08 as a 4th input channel, as it is much less correlated with the Ash colors and pretty well correlated with the contrails, meaning it could add new and useful information, I would not add all 3 as they are highly correlated with each other.","metadata":{}},{"cell_type":"markdown","source":"# Dataset Creation","metadata":{}},{"cell_type":"code","source":"train_df.shape, valid_df.shape","metadata":{"execution":{"iopub.status.busy":"2023-08-21T11:55:34.020939Z","iopub.execute_input":"2023-08-21T11:55:34.021411Z","iopub.status.idle":"2023-08-21T11:55:34.03076Z","shell.execute_reply.started":"2023-08-21T11:55:34.021374Z","shell.execute_reply":"2023-08-21T11:55:34.029623Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Save csvs","metadata":{}},{"cell_type":"code","source":"# Save the csvs\ntrain_df.to_csv('train_df.csv', index=False)\nvalid_df.to_csv('valid_df.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2023-08-21T11:55:34.033482Z","iopub.execute_input":"2023-08-21T11:55:34.034009Z","iopub.status.idle":"2023-08-21T11:55:34.097614Z","shell.execute_reply.started":"2023-08-21T11:55:34.033958Z","shell.execute_reply":"2023-08-21T11:55:34.095381Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Save Images as Numpy","metadata":{}},{"cell_type":"code","source":"def read_record(record_id, directory, mode):\n    record_data = {}\n    read = [\"band_08\", \"band_11\", \"band_14\", \"band_15\"]\n    if mode == 'train':\n        read.append('human_individual_masks')\n    elif mode == 'val':\n        read.append('human_pixel_masks')\n    for x in read:\n        if x == 'human_individual_masks':\n            individual = np.load(os.path.join(directory, record_id, x + \".npy\"))\n            record_data['human_pixel_masks'] = individual.sum(axis=3) / individual.shape[3]\n        else:\n            record_data[x] = np.load(os.path.join(directory, record_id, x + \".npy\"))\n    return record_data","metadata":{"execution":{"iopub.status.busy":"2023-08-21T11:55:34.099038Z","iopub.execute_input":"2023-08-21T11:55:34.099424Z","iopub.status.idle":"2023-08-21T11:55:34.108616Z","shell.execute_reply.started":"2023-08-21T11:55:34.099393Z","shell.execute_reply":"2023-08-21T11:55:34.106957Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def normalize_range(data, bounds):\n    return (data - bounds[0]) / (bounds[1] - bounds[0])\n\ndef get_false_color(record_data):\n    _T11_BOUNDS = (243, 303)\n    _CLOUD_TOP_TDIFF_BOUNDS = (-4, 5)\n    _TDIFF_BOUNDS = (-4, 2)\n    N_TIMES_BEFORE = 4\n\n    r = normalize_range(record_data[\"band_15\"] - record_data[\"band_14\"], _TDIFF_BOUNDS)\n    g = normalize_range(record_data[\"band_14\"] - record_data[\"band_11\"], _CLOUD_TOP_TDIFF_BOUNDS)\n    b = normalize_range(record_data[\"band_14\"], _T11_BOUNDS)\n    n = (record_data[\"band_08\"] - 230) / 20\n\n    false_color = np.stack([r, g, b, n], axis=2)\n    img = false_color[..., N_TIMES_BEFORE]\n    #print(img.shape, img.mean(axis=(0,1)), img.min(axis=(0,1)), img.max(axis=(0,1)))\n    \n    return img","metadata":{"execution":{"iopub.status.busy":"2023-08-21T11:55:34.112537Z","iopub.execute_input":"2023-08-21T11:55:34.113151Z","iopub.status.idle":"2023-08-21T11:55:34.124847Z","shell.execute_reply.started":"2023-08-21T11:55:34.113114Z","shell.execute_reply":"2023-08-21T11:55:34.122364Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"path = Path('contrails')\npath.mkdir(exist_ok=True, parents=True)","metadata":{"execution":{"iopub.status.busy":"2023-08-21T11:55:34.126888Z","iopub.execute_input":"2023-08-21T11:55:34.127291Z","iopub.status.idle":"2023-08-21T11:55:34.144261Z","shell.execute_reply.started":"2023-08-21T11:55:34.127247Z","shell.execute_reply":"2023-08-21T11:55:34.142974Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Val\nfor i in tqdm(valid_rs):\n    data = read_record(str(i), data_dir+'validation', mode='val')\n    img = get_false_color(data)\n    final = np.dstack([img, data['human_pixel_masks']])\n    final = final.astype(np.float16)\n    \n    pathc = path/f\"{i}.npy\"\n    np.save(str(pathc), final)","metadata":{"execution":{"iopub.status.busy":"2023-08-21T11:55:34.146026Z","iopub.execute_input":"2023-08-21T11:55:34.14662Z","iopub.status.idle":"2023-08-21T12:00:31.745408Z","shell.execute_reply.started":"2023-08-21T11:55:34.146573Z","shell.execute_reply":"2023-08-21T12:00:31.743743Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Train\nfor i in tqdm(train_rs):\n    data = read_record(str(i), data_dir+'train', mode='train')\n    img = get_false_color(data)\n    final = np.dstack([img, data['human_pixel_masks']])\n    final = final.astype(np.float16)\n\n    pathc = path/f\"{i}.npy\"\n    np.save(str(pathc), final)","metadata":{"execution":{"iopub.status.busy":"2023-08-21T12:00:31.74809Z","iopub.execute_input":"2023-08-21T12:00:31.748523Z","iopub.status.idle":"2023-08-21T12:59:46.611983Z","shell.execute_reply.started":"2023-08-21T12:00:31.748491Z","shell.execute_reply":"2023-08-21T12:59:46.60714Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}