{"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":"# Introduction","metadata":{}},{"cell_type":"markdown","source":"### Whate are Contrails","metadata":{}},{"cell_type":"markdown","source":"Contrails, short for ‘condensation trails’, are line-shaped clouds of ice crystals that form in aircraft engine exhaust, and are created by airplanes flying through super humid regions in the atmosphere. Persistent contrails contribute as much to global warming as the fuel they burn for flights.\n\n<img style=\"float:center\" src=\"https://storage.googleapis.com/kaggle-media/competitions/Google-Contrails/waterdroplets.png\" width=\"70%\" height=\"70%\">","metadata":{}},{"cell_type":"markdown","source":"#### Contrails explained (video):\n<a href=\"https://www.youtube.com/shorts/_tWnvQfMfrM\n\" target=\"_blank\"><img src=\"https://i.ytimg.com/vi/_tWnvQfMfrM/hq720_2.jpg?sqp=-oaymwE2COgCEMoBSFXyq4qpAygIARUAAIhCGABwAcABBvABAfgBlAOAAtAFigIMCAAQARgWIEoofzAP&rs=AOn4CLDAs-1OA9Oy4gvHM6fhdvAa5gh3gQ\" \nalt=\"IMAGE ALT TEXT HERE\" width=\"750\" height=\"500\" border=\"0\" /></a>","metadata":{"execution":{"iopub.status.busy":"2023-08-03T20:18:06.699489Z","iopub.execute_input":"2023-08-03T20:18:06.699964Z","iopub.status.idle":"2023-08-03T20:18:06.712083Z","shell.execute_reply.started":"2023-08-03T20:18:06.699931Z","shell.execute_reply":"2023-08-03T20:18:06.710304Z"}}},{"cell_type":"markdown","source":"### Goal of the Competition","metadata":{"execution":{"iopub.status.busy":"2023-08-03T20:35:06.412348Z","iopub.execute_input":"2023-08-03T20:35:06.412939Z","iopub.status.idle":"2023-08-03T20:35:06.420883Z","shell.execute_reply.started":"2023-08-03T20:35:06.412892Z","shell.execute_reply":"2023-08-03T20:35:06.419543Z"}}},{"cell_type":"markdown","source":"The aim of the work is to identify aviation contrails using geostationary satellite images. The original full-disk images (<a rel=\"noreferrer nofollow\" target=\"_blank\" href=\"https://console.cloud.google.com/storage/browser/gcp-public-data-goes-16/\">Google Cloud Storage</a>) were reprojected using bilinear resampling to generate a local scene image. Because contrails are easier to identify with temporal context, a sequence of images at 10-minute intervals are provided. Each example (record_id) contains exactly one labeled frame.\n\nLearn more about the dataset from the preprint: <a rel=\"noreferrer nofollow\" target=\"_blank\" href=\"https://arxiv.org/abs/2304.02122\">OpenContrails: Benchmarking Contrail Detection on GOES-16 ABI</a>. Labeling instructions can be found at in this <a rel=\"noreferrer nofollow\" target=\"_blank\" href=\"https://storage.googleapis.com/goes_contrails_dataset/20230419/Contrail_Detection_Dataset_Instruction.pdf\">supplementary material</a>. Some key labeling guidance:\n\n- Contrails must contain at least 10 pixels\n- At some time in their life, Contrails must be at least 3x longer than they are wide\n- Contrails must either appear suddenly or enter from the sides of the image\n- Contrails should be visible in at least two image\n\nThis competition is evaluated on the global <a rel=\"noreferrer nofollow\" href=\"https://en.wikipedia.org/wiki/S%C3%B8rensen%E2%80%93Dice_coefficient\">Dice coefficient</a>. The Dice coefficient can be used to compare the pixel-wise agreement between a predicted segmentation and its corresponding ground truth. The formula is given by:\n\n$$\\frac{2 * |X \\cap Y|}{|X| + |Y|},$$\n\nwhere X is the entire set of predicted contrail pixels for all observations in the test data and Y is the ground truth set of all contrail pixels in the test data.","metadata":{}},{"cell_type":"markdown","source":"# Import dependencies","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport tensorflow as tf\nfrom matplotlib import animation\nimport matplotlib.pyplot as plt\nfrom IPython import display\nimport os","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-08-03T20:44:33.748461Z","iopub.execute_input":"2023-08-03T20:44:33.748901Z","iopub.status.idle":"2023-08-03T20:44:43.805711Z","shell.execute_reply.started":"2023-08-03T20:44:33.748869Z","shell.execute_reply":"2023-08-03T20:44:43.804466Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# import torch\n# import torch.nn as nn","metadata":{"execution":{"iopub.status.busy":"2023-08-03T20:44:43.807807Z","iopub.execute_input":"2023-08-03T20:44:43.808519Z","iopub.status.idle":"2023-08-03T20:44:43.813775Z","shell.execute_reply.started":"2023-08-03T20:44:43.808485Z","shell.execute_reply":"2023-08-03T20:44:43.812309Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Exploratory Data Analysis","metadata":{"execution":{"iopub.status.busy":"2023-05-21T08:50:15.602275Z","iopub.execute_input":"2023-05-21T08:50:15.602705Z","iopub.status.idle":"2023-05-21T08:50:15.635024Z","shell.execute_reply.started":"2023-05-21T08:50:15.602674Z","shell.execute_reply":"2023-05-21T08:50:15.634073Z"}}},{"cell_type":"markdown","source":"### Files","metadata":{}},{"cell_type":"markdown","source":"In each subdirectory named by `{record_id}`, binary files in numpy `.npy` format that corresponds to a single example are provided:\n\n* **`band_{08-16}.npy`**: array with size of `H x W x T`, where `T = n_times_before + n_times_after + 1`, representing the number of images in the sequence. There are `n_times_before` and `n_times_after` images before and after the labeled frame respectively. In our dataset all examples have  `n_times_before=4` and `n_times_after=3`. Each band represents an infrared channel at different wavelengths and is converted to brightness temperatures based on the calibration parameters. The number in the filename corresponds to the GOES-16 ABI band number. Details of the ABI bands can be found [here](https://www.goes-r.gov/mission/ABI-bands-quick-info.html).\n* **`human_individual_masks.npy`**: array with size of `H x W x 1 x R`. Each example is labeled by `R` individual human labelers. `R` is not the same for all samples. The labeled masks have value either 0 or 1 and correspond to the `(n_times_before+1)`-th image in `band_{08-16}.npy`. They are available only in the training set.\n* **`human_pixel_masks.npy`**: array with size of `H x W x 1` containing the  binary groundtruth. A pixel is regarded as contrail pixel in evaluation if it is labeled as contrail by more than half of the labelers. \n\n**`{train/validation}_metadata.json`**: contains the timestamps and the projection parameters to reproduce the satellite images.","metadata":{}},{"cell_type":"code","source":"# Files in record_id folder\n\nBASE_DIR = '/kaggle/input/google-research-identify-contrails-reduce-global-warming'\nN_TIMES_BEFORE = 4\nrecord_id = '1704010292581573769'\n\nfor dirname, _, filenames in os.walk(os.path.join(BASE_DIR, 'train', record_id)):\n    for filename in sorted(filenames):\n        print(os.path.join(dirname, filename))","metadata":{"execution":{"iopub.status.busy":"2023-08-03T20:44:43.815785Z","iopub.execute_input":"2023-08-03T20:44:43.816294Z","iopub.status.idle":"2023-08-03T20:44:43.843782Z","shell.execute_reply.started":"2023-08-03T20:44:43.81625Z","shell.execute_reply":"2023-08-03T20:44:43.842063Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Number of samples per each set\n\nfor set_ in ['train', 'validation', 'test']:\n    print(\"{0} samples: {1}\".format(set_, len(os.listdir(os.path.join(BASE_DIR, set_)))))","metadata":{"execution":{"iopub.status.busy":"2023-08-03T20:44:43.846958Z","iopub.execute_input":"2023-08-03T20:44:43.847335Z","iopub.status.idle":"2023-08-03T20:44:44.136969Z","shell.execute_reply.started":"2023-08-03T20:44:43.847305Z","shell.execute_reply":"2023-08-03T20:44:44.135572Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Visualizing","metadata":{}},{"cell_type":"markdown","source":"Contrails in GOES (Geostationary Operational Environmental Satellites) can be visualized by <b>combining bands into a false color image</b> using \"ash\" color scheme. This color scheme was originally developed for viewing volcanic ash in the atmosphere but is also useful for viewing thin cirrus, including contrails. In this color scheme, contrails appear in the image as dark blue. Better way is to use modified color scheme shown in <a rel=\"noreferrer nofollow\" target=\"_blank\" href=\"https://eumetrain.org/sites/default/files/2020-05/RGB_recipes.pdf\">Compilation of RGB Recipes</a> (page 7).","metadata":{"execution":{"iopub.status.busy":"2023-05-21T09:24:49.43397Z","iopub.execute_input":"2023-05-21T09:24:49.43447Z","iopub.status.idle":"2023-05-21T09:24:49.442875Z","shell.execute_reply.started":"2023-05-21T09:24:49.434436Z","shell.execute_reply":"2023-05-21T09:24:49.441799Z"}}},{"cell_type":"code","source":"# Load bands\n\nwith open(os.path.join(BASE_DIR, 'train', record_id, 'band_11.npy'), 'rb') as f:\n    band11 = np.load(f)\nwith open(os.path.join(BASE_DIR, 'train', record_id, 'band_14.npy'), 'rb') as f:\n    band14 = np.load(f)\nwith open(os.path.join(BASE_DIR, 'train', record_id, 'band_15.npy'), 'rb') as f:\n    band15 = np.load(f)\nwith open(os.path.join(BASE_DIR, 'train', record_id, 'human_pixel_masks.npy'), 'rb') as f:\n    human_pixel_mask = np.load(f)\nwith open(os.path.join(BASE_DIR, 'train', record_id, 'human_individual_masks.npy'), 'rb') as f:\n    human_individual_mask = np.load(f)","metadata":{"execution":{"iopub.status.busy":"2023-08-03T20:44:44.138857Z","iopub.execute_input":"2023-08-03T20:44:44.139378Z","iopub.status.idle":"2023-08-03T20:44:44.246323Z","shell.execute_reply.started":"2023-08-03T20:44:44.139337Z","shell.execute_reply":"2023-08-03T20:44:44.244984Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Combine bands to RGB\n\n_T11_BOUNDS = (243, 303)\n_CLOUD_TOP_TDIFF_BOUNDS = (-4, 5)\n_TDIFF_BOUNDS = (-4, 2)\n\ndef normalize_range(data, bounds):\n    \"\"\"Maps data to the range [0, 1].\"\"\"\n    return (data - bounds[0]) / (bounds[1] - bounds[0])\n\nr = normalize_range(band15 - band14, _TDIFF_BOUNDS)\ng = normalize_range(band14 - band11, _CLOUD_TOP_TDIFF_BOUNDS)\nb = normalize_range(band14, _T11_BOUNDS)\nfalse_color = np.clip(np.stack([r, g, b], axis=2), 0, 1)","metadata":{"execution":{"iopub.status.busy":"2023-08-03T20:44:44.24791Z","iopub.execute_input":"2023-08-03T20:44:44.248257Z","iopub.status.idle":"2023-08-03T20:44:44.283268Z","shell.execute_reply.started":"2023-08-03T20:44:44.248228Z","shell.execute_reply":"2023-08-03T20:44:44.282052Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Create false color image with contrail mask\n\nimg = false_color[..., N_TIMES_BEFORE]\n\nplt.figure(figsize=(16, 6))\nax = plt.subplot(1, 3, 1)\nax.imshow(img)\nax.set_title('False color image')\n\nax = plt.subplot(1, 3, 2)\nax.imshow(human_pixel_mask, interpolation='none', cmap='Blues')\nax.set_title('Ground truth contrail mask')\n\nax = plt.subplot(1, 3, 3)\nax.imshow(img)\nax.imshow(human_pixel_mask, cmap='Reds', alpha=.4, interpolation='none')\nax.set_title('Contrail mask on false color image');","metadata":{"execution":{"iopub.status.busy":"2023-08-03T20:44:44.285139Z","iopub.execute_input":"2023-08-03T20:44:44.285645Z","iopub.status.idle":"2023-08-03T20:44:45.280762Z","shell.execute_reply.started":"2023-08-03T20:44:44.285604Z","shell.execute_reply":"2023-08-03T20:44:45.279596Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Individual human masks\n\nn = human_individual_mask.shape[-1]\nplt.figure(figsize=(16, 4))\nplt.suptitle('Human individual mask')\nfor i in range(n):\n    plt.subplot(1, n, i+1)\n    plt.imshow(human_individual_mask[..., i], interpolation='none', cmap='Blues')","metadata":{"execution":{"iopub.status.busy":"2023-08-03T20:44:45.282458Z","iopub.execute_input":"2023-08-03T20:44:45.283097Z","iopub.status.idle":"2023-08-03T20:44:45.983703Z","shell.execute_reply.started":"2023-08-03T20:44:45.283064Z","shell.execute_reply":"2023-08-03T20:44:45.982504Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"human_individual_mask_aggregated = human_individual_mask.sum(axis=3)\n\nplt.figure(figsize=(8, 4))\nplt.suptitle('Aggregating human individual masks into groundtruth contrail')\nax = plt.subplot(1, 2, 1)\nax.imshow(human_individual_mask_aggregated, cmap='Blues')\nax.set_title('Sum of human individual masks')\n\nax = plt.subplot(1, 2, 2)\nax.imshow(human_individual_mask_aggregated > n/2, interpolation='none', cmap='Blues')\nax.set_title('Groundtruth contrail')","metadata":{"execution":{"iopub.status.busy":"2023-08-03T20:44:45.985517Z","iopub.execute_input":"2023-08-03T20:44:45.986803Z","iopub.status.idle":"2023-08-03T20:44:46.454756Z","shell.execute_reply.started":"2023-08-03T20:44:45.986763Z","shell.execute_reply":"2023-08-03T20:44:46.453577Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Animation\n\nfig = plt.figure(figsize=(6, 6))\nim = plt.imshow(false_color[..., 0])\ndef draw(i):\n    im.set_array(false_color[..., i])\n    return [im]\nanim = animation.FuncAnimation(\n    fig, draw, frames=false_color.shape[-1], interval=500, blit=True\n)\nplt.close()\ndisplay.HTML(anim.to_jshtml())","metadata":{"execution":{"iopub.status.busy":"2023-08-03T20:44:46.457965Z","iopub.execute_input":"2023-08-03T20:44:46.458333Z","iopub.status.idle":"2023-08-03T20:44:48.909814Z","shell.execute_reply.started":"2023-08-03T20:44:46.458302Z","shell.execute_reply":"2023-08-03T20:44:48.906674Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"band_list = ['08', '09', '10', '11', '12', '13', '14', '15', '16']\n\nfig, ax = plt.subplots(len(band_list), 8, figsize=(16,24))\n\nfor i, band in enumerate(band_list):\n    img = np.load(os.path.join(BASE_DIR, 'train', record_id, 'band_' + band + '.npy'))\n    for t in range(8):\n        ax[i, t].imshow(img[..., t], cmap='GnBu')\n        ax[i, t].set_title('Band {0}\\nTimestep {1}'.format(band, t+1))\n        ax[i, t].axis('off')","metadata":{"execution":{"iopub.status.busy":"2023-08-03T20:52:19.189254Z","iopub.execute_input":"2023-08-03T20:52:19.189782Z","iopub.status.idle":"2023-08-03T20:52:26.364684Z","shell.execute_reply.started":"2023-08-03T20:52:19.189727Z","shell.execute_reply":"2023-08-03T20:52:26.362931Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Check if dataset is imbalanced","metadata":{}},{"cell_type":"code","source":"# Count positive (contrails) and negative (no contrails) samples per train and validation set\n\ncounter = {'train_0': 0, 'train_1': 0, 'validation_0': 0, 'validation_1': 0}\n\nfor set_ in ['train', 'validation']:\n    for record in os.listdir(os.path.join(BASE_DIR, set_)):\n        with open(os.path.join(BASE_DIR, set_, record, 'human_pixel_masks.npy'), 'rb') as f:\n            human_pixel_mask = np.load(f)\n            if human_pixel_mask.any():\n                counter[set_ + '_1'] += 1\n            else:\n                counter[set_ + '_0'] += 1\n\n# counter = {'train_0': 11270, 'train_1': 9259, 'validation_0': 1304, 'validation_1': 552} # result of above code\n\n\n# Chart\n\nplt.figure(figsize=(10, 4))\nax = plt.subplot(1, 2, 1)\nax.bar(x=['Negative (no contrails)', 'Positive (contrails)'], height=[counter['train_0'], counter['train_1']])\nax.set_title('Train set')\n\nax = plt.subplot(1, 2, 2)\nax.bar(x=['Negative (no contrails)', 'Positive (contrails)'], height=[counter['validation_0'], counter['validation_1']])\nax.set_title('Validation set')","metadata":{"execution":{"iopub.status.busy":"2023-08-03T20:44:55.883138Z","iopub.execute_input":"2023-08-03T20:44:55.883632Z","iopub.status.idle":"2023-08-03T20:48:07.621742Z","shell.execute_reply.started":"2023-08-03T20:44:55.883601Z","shell.execute_reply":"2023-08-03T20:48:07.620332Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can try to increase positive samples in validation set using data augmentation techniques.","metadata":{}}]}