GEE_LandSat8_Volcano_Download_script.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "665380b9",
      "metadata": {
        "id": "665380b9"
      },
      "source": [
        "# Download representative images of volcanos\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "dcfacdb8",
      "metadata": {
        "id": "dcfacdb8"
      },
      "source": [
        "Credit: this script derived from [Aaron Geller's blog post.](https://sites.northwestern.edu/researchcomputing/2021/11/19/downloading-satellite-images-made-easy/), including the use of functions based on his [EarthEngineToGeoTIFF repository.](https://github.com/ageller/EarthEngineToGeoTIFF)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "74c88e1d",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "74c88e1d",
        "outputId": "109820f1-3fe7-4c2b-a96b-1bf6b394d640"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/\n",
            "Requirement already satisfied: rasterio in /usr/local/lib/python3.7/dist-packages (1.2.10)\n",
            "Requirement already satisfied: click>=4.0 in /usr/local/lib/python3.7/dist-packages (from rasterio) (7.1.2)\n",
            "Requirement already satisfied: affine in /usr/local/lib/python3.7/dist-packages (from rasterio) (2.3.1)\n",
            "Requirement already satisfied: certifi in /usr/local/lib/python3.7/dist-packages (from rasterio) (2022.6.15)\n",
            "Requirement already satisfied: cligj>=0.5 in /usr/local/lib/python3.7/dist-packages (from rasterio) (0.7.2)\n",
            "Requirement already satisfied: click-plugins in /usr/local/lib/python3.7/dist-packages (from rasterio) (1.1.1)\n",
            "Requirement already satisfied: setuptools in /usr/local/lib/python3.7/dist-packages (from rasterio) (57.4.0)\n",
            "Requirement already satisfied: snuggs>=1.4.1 in /usr/local/lib/python3.7/dist-packages (from rasterio) (1.4.7)\n",
            "Requirement already satisfied: attrs in /usr/local/lib/python3.7/dist-packages (from rasterio) (21.4.0)\n",
            "Requirement already satisfied: numpy in /usr/local/lib/python3.7/dist-packages (from rasterio) (1.21.6)\n",
            "Requirement already satisfied: pyparsing>=2.1.6 in /usr/local/lib/python3.7/dist-packages (from snuggs>=1.4.1->rasterio) (3.0.9)\n"
          ]
        }
      ],
      "source": [
        "import ee\n",
        "%pip install rasterio\n",
        "\n",
        "import rasterio\n",
        "from rasterio.plot import show as showRasterio\n",
        "\n",
        "import pandas as pd\n",
        "import re\n",
        "import matplotlib.pyplot as plt\n",
        "%matplotlib inline\n",
        "\n",
        "from datetime import date\n",
        "from datetime import timedelta\n",
        "import zipfile\n",
        "import os\n",
        "import requests\n",
        "\n",
        "import numpy as np\n"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# connect to google drive\n",
        "from google.colab import drive  \n",
        "\n",
        "# drive roots\n",
        "DRIVE_MOUNT = '/content/drive'\n",
        "drive.mount(DRIVE_MOUNT, force_remount=True)\n",
        "\n",
        "# set up project directories\n",
        "PROJECT_DIR = DRIVE_MOUNT + \"/My Drive/Volcano_datasets/\"\n",
        "OUTPUT_IMAGE_DIR = DRIVE_MOUNT +\"/My Drive/Full_landsat8_10km/\"\n",
        "PRIOR_OUTPUT_IMAGE_DIR = PROJECT_DIR +\"Full_landsat8_10km_prior_year/\"\n",
        "\n",
        "TIFS_DIR = DRIVE_MOUNT+\"/My Drive/tifs/\"\n",
        "\n",
        "lst = os.listdir(OUTPUT_IMAGE_DIR)\n",
        "print(lst)\n",
        "\n",
        "targets = pd.read_excel(PROJECT_DIR + \"Sample Eruptions.xlsx\")\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "A1iJapXtdlHe",
        "outputId": "dbe7ab0d-57c4-4787-aa82-fe65a9bca560"
      },
      "id": "A1iJapXtdlHe",
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Mounted at /content/drive\n",
            "[]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "8396951f",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "8396951f",
        "outputId": "089f8419-596a-4439-fc77-a83145aaa74e"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "To authorize access needed by Earth Engine, open the following URL in a web browser and follow the instructions. If the web browser does not start automatically, please manually browse the URL below.\n",
            "\n",
            "    https://code.earthengine.google.com/client-auth?scopes=https%3A//www.googleapis.com/auth/earthengine%20https%3A//www.googleapis.com/auth/devstorage.full_control&request_id=IhuUigsHuqywWg_QxzMmn6UpzcM1e9cF-NrgS6M3M-A&tc=4_YnE_FYZ_7unr-V6Jxelrqv9HJta7Q-6KmRuJDcEWQ&cc=575N9tAn_7fzht222xM35xKPACrNcMn3122TlN1CZBU\n",
            "\n",
            "The authorization workflow will generate a code, which you should paste in the box below. \n",
            "Enter verification code: 4/1AX4XfWjzFLLGQQ5bY2Io5SO8oAQ9omf1pv4QOszM4x7BbDiA8_S4d284d9M\n",
            "\n",
            "Successfully saved authorization token.\n"
          ]
        }
      ],
      "source": [
        "# Trigger the authentication flow.\n",
        "ee.Authenticate()\n",
        "\n",
        "# Initialize the library.\n",
        "ee.Initialize()"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "\n",
        "import time\n",
        "\n",
        "sleep_interval = 20\n",
        "wait_limit = 900\n",
        "download = True\n",
        "\n",
        "def get_landsat8_image(lon, lat, sze_m, filename_base, \n",
        "                         dateMin = '2020-04-01', dateMax = '2020-04-30', \n",
        "                         temp_dir = \"tifs\",\n",
        "                           output_dir = OUTPUT_IMAGE_DIR):\n",
        "    '''    \n",
        "    download an SARimage from the Sentinal S2 Surface Reflectance satellite, at the given coordinates\n",
        "    \n",
        "    lon : central longitude in degrees\n",
        "    lat : central latitude in degrees\n",
        "    sze_km : size of the edge of the box in km\n",
        "    dateMin : minimum date to use for image search in year-month-day (e.g., 2020-08-01)\n",
        "    dateMax : maximum date to use for image search in year-month-day (e.g., 2020-08-31)\n",
        "    base_filename : base output filename for the GeoTIFF images. Will create 6 files for each bse\n",
        "     For each of ASCENDING and DECENDING:\n",
        "      - mean of VH channel\n",
        "      - mean of VV channel\n",
        "     mean of Ascending and descending:\n",
        "      - VV channel\n",
        "      - VH channel \n",
        "    \n",
        "    https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C02_T1_TOA\n",
        "    '''\n",
        "\n",
        "    print(f\"{lon} {lat} {filename_base} {dateMin} {dateMax} {temp_dir} {output_dir}\")\n",
        "    download = True\n",
        "    if (download):\n",
        "      # Calculate Area of Interest\n",
        "      pt = ee.Geometry.Point(lon,lat)\n",
        "          \n",
        "      aoi = pt.buffer(distance=sze_m)\n",
        "\n",
        "      # Load the LandSat data\n",
        "      LS8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')\\\n",
        "                               .filterBounds(aoi)\\\n",
        "                               .filterDate(dateMin,dateMax)\\\n",
        "                               .filter(ee.Filter.lt('CLOUD_COVER', 15))\\\n",
        "                               .mean()\\\n",
        "                               .set('SENSOR_ID', 'OLI_TIRS')\n",
        "\n",
        "      # Bands to output\n",
        "      #B1 \t30 meters \t0.43 - 0.45 μm \tCoastal aerosol\n",
        "      #B2 \t30 meters \t0.45 - 0.51 μm \tBlue\n",
        "      #B3 \t30 meters \t0.53 - 0.59 μm \tGreen\n",
        "      #B4 \t30 meters \t0.64 - 0.67 μm \tRed\n",
        "      #B5 \t30 meters \t0.85 - 0.88 μm \tNear infrared\n",
        "      #B6 \t30 meters \t1.57 - 1.65 μm \tShortwave infrared 1\n",
        "      #B7 \t30 meters \t2.11 - 2.29 μm \tShortwave infrared 2\n",
        "      #B8 \t15 meters \t0.52 - 0.90 μm \tBand 8 Panchromatic\n",
        "      #B9 \t15 meters \t1.36 - 1.38 μm \tCirrus\n",
        "      #B10 \t30 meters \t10.60 - 11.19 μm Thermal infrared 1, resampled from 100m to 30m\n",
        "      #B11 \t30 meters \t11.50 - 12.51 μm Thermal infrared 2, resampled from 100m to 30m\n",
        "      bands = [\"B2\",\"B3\",\"B4\",\"B5\",\"B6\",\"B7\",\"B9\",\"B10\",\"B11\"]\n",
        "\n",
        "      task_list = []\n",
        "      for selection in bands:\n",
        "          #print(f\"Task: {filename_base} {selection}\")\n",
        "          task = ee.batch.Export.image.toDrive(image=LS8.select(selection),\n",
        "                                            description=filename_base+\"_\"+selection,\n",
        "                                            scale=30,\n",
        "                                            region=aoi,\n",
        "                                            fileNamePrefix=filename_base+\"_\"+selection,\n",
        "                                            folder=temp_dir,\n",
        "                                            crs='EPSG:4326',\n",
        "                                            fileFormat='GeoTIFF')\n",
        "\n",
        "          task.start()\n",
        "          task_list.append(task)\n",
        "      \n",
        "      loop_count= 0\n",
        "\n",
        "      for task in task_list:\n",
        "        while(task.active() and sleep_interval*loop_count < wait_limit):\n",
        "           loop_count += 1\n",
        "           time.sleep(sleep_interval)\n",
        "           #print(f'slept {sleep_interval*loop_count}')\n",
        "          \n",
        "        status = task.status()\n",
        "\n",
        "        if (status['state'] == 'COMPLETED' or status['state'] == 'READY'):\n",
        "            print(f\"Success: {task}\")#print(status['destination_uris'])\n",
        "        else:\n",
        "            print (status['state'])\n",
        "            if (status['state'] == \"FAILED\"):\n",
        "               print(f\"Failed: {task} {status['error_message']}\")\n",
        "\n",
        "      time.sleep(60) # sleep for final file to write\n",
        "\n",
        "    # Group bands by output files\n",
        "    bandsets = {\"rgb\" : [\"B4\",\"B3\",\"B2\"],\n",
        "                \"B5_B6_B7\" : [\"B5\",\"B6\",\"B7\"],\n",
        "                \"B9_B10_B11\" : [\"B9\",\"B10\",\"B11\"]}\n",
        "\n",
        "    for bandset in bandsets.keys():\n",
        "      # convert image to png\n",
        "       f,ax = plt.subplots(figsize=(15,15))\n",
        "       output = OUTPUT_IMAGE_DIR+filename_base+\"_\"+bandset+\".png\"\n",
        "       for iband, band in enumerate(bandsets[bandset]):\n",
        "         rast = rasterio.open(TIFS_DIR+filename_base+\"_\"+band+\".tif\")\n",
        "         img = rast.read()\n",
        "         if (iband == 0):\n",
        "            nz, ny, nx = img.shape\n",
        "            image = np.zeros((3,ny,nx))\n",
        "            transform=rast.transform\n",
        "         image[iband,:,:] = img[0,:,:]\n",
        "         rast.close()\n",
        "\n",
        "       showRasterio(image, ax = ax, transform=transform)    \n",
        "       f.savefig(output, bbox_inches='tight')\n",
        "       plt.close(f)\n",
        "\n",
        "       print(f\"wrote {output}\")\n",
        "\n",
        "    # clean up tifs left behind\n",
        "    lst = os.listdir(TIFS_DIR)\n",
        "   \n",
        "    regex_pattern = filename_base + \".*.tif\"\n",
        "    for f in lst:\n",
        "       if re.match(regex_pattern, f):\n",
        "         print(f\"removing {f}\")\n",
        "         os.remove(TIFS_DIR+\"/\"+f)\n"
      ],
      "metadata": {
        "id": "hKh3bVUCH_fp"
      },
      "id": "hKh3bVUCH_fp",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "target_size_m = 10000\n",
        "pull_eruption_year = True\n",
        "pull_prior_year = True\n",
        "\n",
        "for i_eruption in range(2,3):#targets.shape[0]):\n",
        "    volc = targets.loc[i_eruption]\n",
        "    v_name = re.sub(r\"[,.]\",\"\",re.sub(r\"\\s+\", \"_\",volc[\"Volcano Name\"]))\n",
        "   \n",
        "\n",
        "    # adjust eruption date by 1 year\n",
        "    erupt_start_date = date(volc[\"Start Year\"].astype(int),\n",
        "                            volc[\"Start Month\"].astype(int),\n",
        "                            volc[\"Start Day\"].astype(int))\n",
        "    erupt_end_date = date(volc[\"End Year\"].astype(int),\n",
        "                          volc[\"End Month\"].astype(int),\n",
        "                          volc[\"End Day\"].astype(int))\n",
        "\n",
        "    if (erupt_start_date >= date(2015,6,23)):\n",
        "      print(v_name,erupt_start_date)\n",
        "      pre_erupt_date = erupt_start_date - timedelta(weeks=6)\n",
        "\n",
        "      if pull_eruption_year:\n",
        "        # do prior to eruption too\n",
        "        \n",
        "        volc_file_name = f'{v_name}__{pre_erupt_date}_pre_{volc[\"Volcano Number\"]:.0f}_{volc[\"Eruption Number\"]:.0f}'\n",
        "        print(i_eruption, volc_file_name)\n",
        "          \n",
        "        try:\n",
        "            _ =get_landsat8_image(volc[\"Longitude\"], volc[\"Latitude\"], target_size_m,\n",
        "                            volc_file_name,\n",
        "                                output_dir = OUTPUT_IMAGE_DIR,\n",
        "                                dateMin = pre_erupt_date.isoformat(), dateMax=erupt_start_date.isoformat())\n",
        "        except Exception as e:\n",
        "            print(f'Failed to get data for {volc_file_name}')\n",
        "            print(e)\n",
        "      \n",
        "      if pull_prior_year:\n",
        "        prior_volc_file_name = f'{v_name}__{pre_erupt_date-timedelta(weeks=52)}_pre_{volc[\"Volcano Number\"]:.0f}_{volc[\"Eruption Number\"]:.0f}'\n",
        "  \n",
        "        print(i_eruption, prior_volc_file_name)\n",
        "          \n",
        "        try:\n",
        "            _ =get_landsat8_image(volc[\"Longitude\"], volc[\"Latitude\"], target_size_m,\n",
        "                            volc_file_name,\n",
        "                               output_dir = PRIOR_OUTPUT_IMAGE_DIR,\n",
        "                               dateMin =  (pre_erupt_date-timedelta(weeks=52)).isoformat(), dateMax=(erupt_start_date-timedelta(weeks=52)).isoformat())\n",
        "        except Exception as e:\n",
        "            print(f'Failed to get data for prior year {prior_volc_file_name}')\n",
        "            print(e)\n",
        "      \n",
        "      next_start_date = erupt_start_date\n",
        "      next_end_date = erupt_start_date + timedelta(days=7)\n",
        "      while (next_end_date <= erupt_end_date + timedelta(days=7) and next_end_date <= erupt_start_date + timedelta(days=187)):\n",
        "          volc_file_name = f'{v_name}_{next_start_date}_{volc[\"Volcano Number\"]:.0f}_{volc[\"Eruption Number\"]:.0f}'\n",
        "          prior_volc_file_name = f'{v_name}_{next_start_date-timedelta(weeks=52)}_{volc[\"Volcano Number\"]:.0f}_{volc[\"Eruption Number\"]:.0f}'\n",
        "\n",
        "          print(i_eruption, volc_file_name)\n",
        "          if pull_eruption_year:\n",
        "            try:\n",
        "                _ =get_landsat8_image(volc[\"Longitude\"], volc[\"Latitude\"], \n",
        "                            target_size_m, volc_file_name,\n",
        "                            output_dir = OUTPUT_IMAGE_DIR,\n",
        "                            dateMin = next_start_date.isoformat(), \n",
        "                            dateMax=next_end_date.isoformat())\n",
        "            except Exception as e:\n",
        "              print(f'Failed to get data for {volc_file_name}:\\n {e}')\n",
        "          if pull_prior_year:\n",
        "            try:\n",
        "                _ =get_landsat8_image(volc[\"Longitude\"], volc[\"Latitude\"], \n",
        "                            target_size_m, prior_volc_file_name,\n",
        "                            output_dir = PRIOR_OUTPUT_IMAGE_DIR,\n",
        "                            dateMin = (next_start_date-timedelta(weeks=52)).isoformat(), \n",
        "                            dateMax=(next_end_date-timedelta(weeks=52)).isoformat())\n",
        "            except Exception as e:\n",
        "              print(f'Failed to get data for {volc_file_name}:\\n {e}')\n",
        "          next_start_date = next_start_date + timedelta(days=7)\n",
        "          next_end_date = next_end_date + timedelta(days=7)\n",
        "        "
      ],
      "metadata": {
        "id": "V3-iEGY8hezd"
      },
      "id": "V3-iEGY8hezd",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "lat = 12.702;\n",
        "lon = -87.004;\n",
        "dateMin = '2020-01-01';\n",
        "dateMax = '2021-1-01';\n",
        "sze_m = 5000;\n",
        "filename_base = \"San_Cristobal\"\n",
        "\n",
        "_ =get_landsat8_image(lon,lat,sze_m,filename_base,dateMin,dateMax)"
      ],
      "metadata": {
        "id": "D--3GU2Be88R"
      },
      "id": "D--3GU2Be88R",
      "execution_count": null,
      "outputs": []
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3 (ipykernel)",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.9.7"
    },
    "colab": {
      "name": "GEE LandSat8 Volcano Download script",
      "provenance": [],
      "collapsed_sections": []
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}