mask2seed.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "machine_shape": "hm",
      "authorship_tag": "ABX9TyPOAu9ZOiwbieu3WQ8S3Yd1",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    },
    "gpuClass": "standard"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/ericslevenson/arctic-surface-water/blob/main/mask2seed.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Ff50g3oRrrFD"
      },
      "outputs": [],
      "source": [
        "'''\n",
        "author: @ericslevenson\n",
        "date: 11/28/2022\n",
        "description: Convert lake occurrence raster to clean buffered lake shapefile to\n",
        "seed the time series in GEE\n",
        "'''"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Authenticate private account (only required for exporting to drive/gee/gcp)\n",
        "from google.colab import auth\n",
        "auth.authenticate_user()\n",
        "# Google Drive setup\n",
        "from google.colab import drive\n",
        "drive.mount('/content/drive')\n",
        "# Complete the environment\n",
        "!pip install rasterio\n",
        "!pip install rioxarray\n",
        "! pip install geopandas"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "_EUA8SYxr_7N",
        "outputId": "6759d097-b046-4141-b46c-69aa6f537e2c"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Mounted at /content/drive\n",
            "Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/\n",
            "Collecting rasterio\n",
            "  Downloading rasterio-1.3.4-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (20.9 MB)\n",
            "\u001b[K     |████████████████████████████████| 20.9 MB 4.3 MB/s \n",
            "\u001b[?25hCollecting affine\n",
            "  Downloading affine-2.3.1-py2.py3-none-any.whl (16 kB)\n",
            "Requirement already satisfied: certifi in /usr/local/lib/python3.8/dist-packages (from rasterio) (2022.9.24)\n",
            "Collecting snuggs>=1.4.1\n",
            "  Downloading snuggs-1.4.7-py3-none-any.whl (5.4 kB)\n",
            "Requirement already satisfied: click>=4.0 in /usr/local/lib/python3.8/dist-packages (from rasterio) (7.1.2)\n",
            "Requirement already satisfied: setuptools in /usr/local/lib/python3.8/dist-packages (from rasterio) (57.4.0)\n",
            "Collecting cligj>=0.5\n",
            "  Downloading cligj-0.7.2-py3-none-any.whl (7.1 kB)\n",
            "Requirement already satisfied: attrs in /usr/local/lib/python3.8/dist-packages (from rasterio) (22.1.0)\n",
            "Collecting click-plugins\n",
            "  Downloading click_plugins-1.1.1-py2.py3-none-any.whl (7.5 kB)\n",
            "Requirement already satisfied: numpy>=1.18 in /usr/local/lib/python3.8/dist-packages (from rasterio) (1.21.6)\n",
            "Requirement already satisfied: pyparsing>=2.1.6 in /usr/local/lib/python3.8/dist-packages (from snuggs>=1.4.1->rasterio) (3.0.9)\n",
            "Installing collected packages: snuggs, cligj, click-plugins, affine, rasterio\n",
            "Successfully installed affine-2.3.1 click-plugins-1.1.1 cligj-0.7.2 rasterio-1.3.4 snuggs-1.4.7\n",
            "Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/\n",
            "Collecting rioxarray\n",
            "  Downloading rioxarray-0.13.1-py3-none-any.whl (53 kB)\n",
            "\u001b[K     |████████████████████████████████| 53 kB 1.2 MB/s \n",
            "\u001b[?25hCollecting pyproj>=2.2\n",
            "  Downloading pyproj-3.4.0-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (7.8 MB)\n",
            "\u001b[K     |████████████████████████████████| 7.8 MB 8.4 MB/s \n",
            "\u001b[?25hRequirement already satisfied: xarray>=0.17 in /usr/local/lib/python3.8/dist-packages (from rioxarray) (0.20.2)\n",
            "Requirement already satisfied: packaging in /usr/local/lib/python3.8/dist-packages (from rioxarray) (21.3)\n",
            "Requirement already satisfied: rasterio>=1.1.1 in /usr/local/lib/python3.8/dist-packages (from rioxarray) (1.3.4)\n",
            "Requirement already satisfied: certifi in /usr/local/lib/python3.8/dist-packages (from pyproj>=2.2->rioxarray) (2022.9.24)\n",
            "Requirement already satisfied: click>=4.0 in /usr/local/lib/python3.8/dist-packages (from rasterio>=1.1.1->rioxarray) (7.1.2)\n",
            "Requirement already satisfied: attrs in /usr/local/lib/python3.8/dist-packages (from rasterio>=1.1.1->rioxarray) (22.1.0)\n",
            "Requirement already satisfied: affine in /usr/local/lib/python3.8/dist-packages (from rasterio>=1.1.1->rioxarray) (2.3.1)\n",
            "Requirement already satisfied: click-plugins in /usr/local/lib/python3.8/dist-packages (from rasterio>=1.1.1->rioxarray) (1.1.1)\n",
            "Requirement already satisfied: numpy>=1.18 in /usr/local/lib/python3.8/dist-packages (from rasterio>=1.1.1->rioxarray) (1.21.6)\n",
            "Requirement already satisfied: setuptools in /usr/local/lib/python3.8/dist-packages (from rasterio>=1.1.1->rioxarray) (57.4.0)\n",
            "Requirement already satisfied: snuggs>=1.4.1 in /usr/local/lib/python3.8/dist-packages (from rasterio>=1.1.1->rioxarray) (1.4.7)\n",
            "Requirement already satisfied: cligj>=0.5 in /usr/local/lib/python3.8/dist-packages (from rasterio>=1.1.1->rioxarray) (0.7.2)\n",
            "Requirement already satisfied: pyparsing>=2.1.6 in /usr/local/lib/python3.8/dist-packages (from snuggs>=1.4.1->rasterio>=1.1.1->rioxarray) (3.0.9)\n",
            "Requirement already satisfied: pandas>=1.1 in /usr/local/lib/python3.8/dist-packages (from xarray>=0.17->rioxarray) (1.3.5)\n",
            "Requirement already satisfied: python-dateutil>=2.7.3 in /usr/local/lib/python3.8/dist-packages (from pandas>=1.1->xarray>=0.17->rioxarray) (2.8.2)\n",
            "Requirement already satisfied: pytz>=2017.3 in /usr/local/lib/python3.8/dist-packages (from pandas>=1.1->xarray>=0.17->rioxarray) (2022.6)\n",
            "Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.8/dist-packages (from python-dateutil>=2.7.3->pandas>=1.1->xarray>=0.17->rioxarray) (1.15.0)\n",
            "Installing collected packages: pyproj, rioxarray\n",
            "Successfully installed pyproj-3.4.0 rioxarray-0.13.1\n",
            "Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/\n",
            "Collecting geopandas\n",
            "  Downloading geopandas-0.12.1-py3-none-any.whl (1.1 MB)\n",
            "\u001b[K     |████████████████████████████████| 1.1 MB 3.9 MB/s \n",
            "\u001b[?25hRequirement already satisfied: pandas>=1.0.0 in /usr/local/lib/python3.8/dist-packages (from geopandas) (1.3.5)\n",
            "Requirement already satisfied: shapely>=1.7 in /usr/local/lib/python3.8/dist-packages (from geopandas) (1.8.5.post1)\n",
            "Collecting fiona>=1.8\n",
            "  Downloading Fiona-1.8.22-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (16.6 MB)\n",
            "\u001b[K     |████████████████████████████████| 16.6 MB 60.3 MB/s \n",
            "\u001b[?25hRequirement already satisfied: packaging in /usr/local/lib/python3.8/dist-packages (from geopandas) (21.3)\n",
            "Requirement already satisfied: pyproj>=2.6.1.post1 in /usr/local/lib/python3.8/dist-packages (from geopandas) (3.4.0)\n",
            "Requirement already satisfied: attrs>=17 in /usr/local/lib/python3.8/dist-packages (from fiona>=1.8->geopandas) (22.1.0)\n",
            "Requirement already satisfied: six>=1.7 in /usr/local/lib/python3.8/dist-packages (from fiona>=1.8->geopandas) (1.15.0)\n",
            "Requirement already satisfied: click>=4.0 in /usr/local/lib/python3.8/dist-packages (from fiona>=1.8->geopandas) (7.1.2)\n",
            "Requirement already satisfied: setuptools in /usr/local/lib/python3.8/dist-packages (from fiona>=1.8->geopandas) (57.4.0)\n",
            "Requirement already satisfied: cligj>=0.5 in /usr/local/lib/python3.8/dist-packages (from fiona>=1.8->geopandas) (0.7.2)\n",
            "Requirement already satisfied: certifi in /usr/local/lib/python3.8/dist-packages (from fiona>=1.8->geopandas) (2022.9.24)\n",
            "Requirement already satisfied: click-plugins>=1.0 in /usr/local/lib/python3.8/dist-packages (from fiona>=1.8->geopandas) (1.1.1)\n",
            "Collecting munch\n",
            "  Downloading munch-2.5.0-py2.py3-none-any.whl (10 kB)\n",
            "Requirement already satisfied: numpy>=1.17.3 in /usr/local/lib/python3.8/dist-packages (from pandas>=1.0.0->geopandas) (1.21.6)\n",
            "Requirement already satisfied: pytz>=2017.3 in /usr/local/lib/python3.8/dist-packages (from pandas>=1.0.0->geopandas) (2022.6)\n",
            "Requirement already satisfied: python-dateutil>=2.7.3 in /usr/local/lib/python3.8/dist-packages (from pandas>=1.0.0->geopandas) (2.8.2)\n",
            "Requirement already satisfied: pyparsing!=3.0.5,>=2.0.2 in /usr/local/lib/python3.8/dist-packages (from packaging->geopandas) (3.0.9)\n",
            "Installing collected packages: munch, fiona, geopandas\n",
            "Successfully installed fiona-1.8.22 geopandas-0.12.1 munch-2.5.0\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import os\n",
        "import rasterio\n",
        "import geopandas\n",
        "import scipy\n",
        "import pandas as pd\n",
        "import numpy as np\n",
        "from rasterio.plot import show\n",
        "from geopandas import GeoDataFrame\n",
        "from scipy import ndimage\n",
        "from pandas import DataFrame\n",
        "from shapely.geometry import shape\n",
        "from rasterio.features import shapes\n",
        "from rasterio.merge import merge\n",
        "import geopandas as gpd"
      ],
      "metadata": {
        "id": "VylB8lkxtR-U"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# input directory\n",
        "LOdirectory = '/content/drive/MyDrive/lakeTransectResults/'\n",
        "outDirectory = '/content/drive/MyDrive/AGU/Lakes75/'\n",
        "\n",
        "\n",
        "IDs = []\n",
        "files = [i for i in os.listdir(LOdirectory) if i.split('.')[1] == 'tif']\n",
        "for f in files:\n",
        "  IDs.append(f.split('_')[0])\n",
        "IDs = list(set(IDs))"
      ],
      "metadata": {
        "id": "ZRIBcyOPxvwj"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "IDs"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "7Rf0rZ4q4m7I",
        "outputId": "fd2eec50-4b12-4e2e-c9a6-e9f205b1af3d"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "['04WEA',\n",
              " '04WEB',\n",
              " '04WES',\n",
              " '04WEU',\n",
              " '04WED',\n",
              " '04WEV',\n",
              " '04VEN',\n",
              " '04WEE',\n",
              " '04VEP',\n",
              " '04VEM',\n",
              " '04WEC',\n",
              " '04VEQ',\n",
              " '04VER',\n",
              " '04WET']"
            ]
          },
          "metadata": {},
          "execution_count": 8
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "##BINARIZE | DILATE | VECTORIZE | EXPORT |\n",
        "\n",
        "# Define Kernel\n",
        "kernel = np.ones((3,3), np.uint8)\n",
        "\n",
        "for i in IDs:\n",
        "  print('id of interest is ' + i)\n",
        "  LOoi = [] # lake occurrence of interest\n",
        "  for f in files:\n",
        "    if f.split('_')[0] == i: # list of files for ID\n",
        "      LOoi.append(f)\n",
        "  if len(LOoi) == 1: # binarize, dilate, vectorize if there's only one\n",
        "    print('processing solo ' + i)\n",
        "    src = rasterio.open(LOdirectory + LOoi[0])\n",
        "    data = src.read(1)\n",
        "    meta = src.meta\n",
        "    trans = src.transform\n",
        "    # binary at 0.05\n",
        "    data = np.where(data > 0.05, 1, 0)\n",
        "    # dilate\n",
        "    dilated = scipy.ndimage.binary_dilation(data, kernel, iterations = 5).astype(np.int32)\n",
        "    # generate shapes\n",
        "    shape_gen = ((shape(s), v) for s, v in shapes(dilated, transform=src.transform))\n",
        "    # build a pd.DataFrame of shapes and convert to geodataframe\n",
        "    df = DataFrame(shape_gen, columns=['geometry', 'class'])\n",
        "    gdf = GeoDataFrame(df[\"class\"], geometry=df.geometry, crs=src.crs)\n",
        "    # filter out non-lake polygons\n",
        "    gdf = gdf[gdf['class'] == 1]\n",
        "    # write to file\n",
        "    gdf.to_file(outDirectory + i + '_seed.shp')\n",
        "  # For split up rasters\n",
        "  elif len(LOoi) == 2:\n",
        "    print('processing as mosaic ' + i)\n",
        "    # Open both rasters\n",
        "    LO1 = rasterio.open(LOdirectory + LOoi[0])\n",
        "    LO2 = rasterio.open(LOdirectory + LOoi[1])\n",
        "    # Merge and retrieve shape\n",
        "    LO, out_trans = merge([LO1, LO2])\n",
        "    LO = LO[0]\n",
        "    width = LO.shape[1]\n",
        "    height = LO.shape[0]\n",
        "    # Metadata and update shape for mosaic\n",
        "    meta = LO1.meta\n",
        "    meta.update({'width':width})\n",
        "    meta.update({'height':height})\n",
        "    LO = np.where(LO > 0.05, 1, 0)\n",
        "    dilated = scipy.ndimage.binary_dilation(LO, kernel, iterations = 5).astype(np.int32)\n",
        "    # generate shapes\n",
        "    shape_gen = ((shape(s), v) for s, v in shapes(dilated, transform=out_trans))\n",
        "    # build a pd.DataFrame of shapes and convert to geodataframe\n",
        "    df = DataFrame(shape_gen, columns=['geometry', 'class'])\n",
        "    gdf = GeoDataFrame(df[\"class\"], geometry=df.geometry, crs=LO1.crs)\n",
        "    # filter out non-lake polygons\n",
        "    gdf = gdf[gdf['class'] == 1]\n",
        "    # write to file\n",
        "    gdf.to_file(outDirectory + i + '_seed.shp')"
      ],
      "metadata": {
        "id": "blClG4CwzMgi"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# BINARIZE | VECTORIZE | EXPORT\n",
        "\n",
        "for i in IDs:\n",
        "  print('id of interest is ' + i)\n",
        "  LOoi = [] # lake occurrence of interest\n",
        "  for f in files:\n",
        "    if f.split('_')[0] == i: # list of files for ID\n",
        "      LOoi.append(f)\n",
        "  if len(LOoi) == 1: # binarize, dilate, vectorize if there's only one\n",
        "    print('processing solo ' + i)\n",
        "    src = rasterio.open(LOdirectory + LOoi[0])\n",
        "    data = src.read(1)\n",
        "    meta = src.meta\n",
        "    trans = src.transform\n",
        "    # binary at 0.05\n",
        "    data = np.where(data > 0.25, 1, 0).astype('int32')\n",
        "    # generate shapes\n",
        "    shape_gen = ((shape(s), v) for s, v in shapes(data, transform=src.transform))\n",
        "    # build a pd.DataFrame of shapes and convert to geodataframe\n",
        "    df = DataFrame(shape_gen, columns=['geometry', 'class'])\n",
        "    gdf = GeoDataFrame(df[\"class\"], geometry=df.geometry, crs=src.crs)\n",
        "    # filter out non-lake polygons\n",
        "    gdf = gdf[gdf['class'] == 1]\n",
        "    # write to file\n",
        "    gdf.to_file(outDirectory + i + '_lake75.shp')\n",
        "  # For split up rasters\n",
        "  elif len(LOoi) == 2:\n",
        "    print('processing as mosaic ' + i)\n",
        "    # Open both rasters\n",
        "    LO1 = rasterio.open(LOdirectory + LOoi[0])\n",
        "    LO2 = rasterio.open(LOdirectory + LOoi[1])\n",
        "    # Merge and retrieve shape\n",
        "    LO, out_trans = merge([LO1, LO2])\n",
        "    LO = LO[0]\n",
        "    width = LO.shape[1]\n",
        "    height = LO.shape[0]\n",
        "    # Metadata and update shape for mosaic\n",
        "    meta = LO1.meta\n",
        "    meta.update({'width':width})\n",
        "    meta.update({'height':height})\n",
        "    LO = np.where(LO > 0.25, 1, 0).astype('int32')\n",
        "    # generate shapes\n",
        "    shape_gen = ((shape(s), v) for s, v in shapes(LO, transform=out_trans))\n",
        "    # build a pd.DataFrame of shapes and convert to geodataframe\n",
        "    df = DataFrame(shape_gen, columns=['geometry', 'class'])\n",
        "    gdf = GeoDataFrame(df[\"class\"], geometry=df.geometry, crs=LO1.crs)\n",
        "    # filter out non-lake polygons\n",
        "    gdf = gdf[gdf['class'] == 1]\n",
        "    # write to file\n",
        "    gdf.to_file(outDirectory + i + '_lake75.shp')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "fgFWx-eDa_X6",
        "outputId": "fb7be856-cdb7-4f77-d9c9-98136069acf3"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "id of interest is 04WEA\n",
            "processing as mosaic 04WEA\n",
            "id of interest is 04WEB\n",
            "processing as mosaic 04WEB\n",
            "id of interest is 04WES\n",
            "processing as mosaic 04WES\n",
            "id of interest is 04WEU\n",
            "processing as mosaic 04WEU\n",
            "id of interest is 04WED\n",
            "processing as mosaic 04WED\n",
            "id of interest is 04WEV\n",
            "processing as mosaic 04WEV\n",
            "id of interest is 04VEN\n",
            "processing solo 04VEN\n",
            "id of interest is 04WEE\n",
            "processing as mosaic 04WEE\n",
            "id of interest is 04VEP\n",
            "processing as mosaic 04VEP\n",
            "id of interest is 04VEM\n",
            "processing solo 04VEM\n",
            "id of interest is 04WEC\n",
            "processing as mosaic 04WEC\n",
            "id of interest is 04VEQ\n",
            "processing as mosaic 04VEQ\n",
            "id of interest is 04VER\n",
            "processing as mosaic 04VER\n",
            "id of interest is 04WET\n",
            "processing as mosaic 04WET\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "#Grid\n",
        "grid = gpd.read_file('/content/drive/MyDrive/AGU/Lakes75/densityGrid_10km.shp') # read in grid\n",
        "grid = grid.set_crs('EPSG:32604') # project to utm\n",
        "grid['full_area'] = grid['geometry'].area # add area field"
      ],
      "metadata": {
        "id": "PeoXlBrhu77B"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# lakes\n",
        "lakes = gpd.read_file('/content/drive/MyDrive/AGU/Lakes75/transect_lake75_utm.shp')"
      ],
      "metadata": {
        "id": "VvfxbzTGvTzS"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "grid_clipped = gpd.clip(grid, lakes)"
      ],
      "metadata": {
        "id": "aJtIXxvhvxJ0"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}