find_lakes.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "accelerator": "GPU",
    "colab": {
      "name": "find_lakes.ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "amKmhTOZqnlK",
        "outputId": "c7a6482a-e12c-41f9-d4d9-3b179c7670fb"
      },
      "source": [
        "import numpy as np\n",
        "from glob import glob\n",
        "import matplotlib.pyplot as plt\n",
        "import gdal\n",
        "import imageio\n",
        "from PIL import Image\n",
        "import json\n",
        "import os\n",
        "from scipy.ndimage import gaussian_filter\n",
        "from skimage import data, img_as_float, measure, exposure\n",
        "from skimage.morphology import reconstruction\n",
        "from skimage.transform import rescale, resize\n",
        "import torch\n",
        "from torchvision import datasets, models, transforms\n",
        "from torch.autograd import Variable\n",
        "from torch import nn, optim, cuda\n",
        "import cv2\n",
        "from cv2 import contourArea\n",
        "from shapely.geometry import Polygon, MultiPolygon, MultiPoint, Point\n",
        "from shapely.ops import nearest_points\n",
        "from shapely.ops import unary_union\n",
        "import scipy as sp\n",
        "import scipy.ndimage\n",
        "import itertools\n",
        "from itertools import combinations\n",
        "\n",
        "!pip install geojson\n",
        "import geojson\n",
        "from geojson import FeatureCollection, Feature, dump\n",
        "\n",
        "!pip install geopandas\n",
        "import geopandas\n",
        "\n",
        "!pip install rasterio\n",
        "from rasterio import mask\n",
        "import rasterio\n",
        "\n",
        "!pip install alphashape\n",
        "import alphashape\n",
        "\n",
        "import gc\n",
        "from shapely import affinity"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Collecting geojson\n",
            "  Downloading https://files.pythonhosted.org/packages/e4/8d/9e28e9af95739e6d2d2f8d4bef0b3432da40b7c3588fbad4298c1be09e48/geojson-2.5.0-py2.py3-none-any.whl\n",
            "Installing collected packages: geojson\n",
            "Successfully installed geojson-2.5.0\n",
            "Collecting geopandas\n",
            "\u001b[?25l  Downloading https://files.pythonhosted.org/packages/f7/a4/e66aafbefcbb717813bf3a355c8c4fc3ed04ea1dd7feb2920f2f4f868921/geopandas-0.8.1-py2.py3-none-any.whl (962kB)\n",
            "\u001b[K     |████████████████████████████████| 972kB 8.7MB/s \n",
            "\u001b[?25hCollecting fiona\n",
            "\u001b[?25l  Downloading https://files.pythonhosted.org/packages/37/94/4910fd55246c1d963727b03885ead6ef1cd3748a465f7b0239ab25dfc9a3/Fiona-1.8.18-cp36-cp36m-manylinux1_x86_64.whl (14.8MB)\n",
            "\u001b[K     |████████████████████████████████| 14.8MB 204kB/s \n",
            "\u001b[?25hCollecting pyproj>=2.2.0\n",
            "\u001b[?25l  Downloading https://files.pythonhosted.org/packages/e4/ab/280e80a67cfc109d15428c0ec56391fc03a65857b7727cf4e6e6f99a4204/pyproj-3.0.0.post1-cp36-cp36m-manylinux2010_x86_64.whl (6.4MB)\n",
            "\u001b[K     |████████████████████████████████| 6.5MB 50.7MB/s \n",
            "\u001b[?25hRequirement already satisfied: shapely in /usr/local/lib/python3.6/dist-packages (from geopandas) (1.7.1)\n",
            "Requirement already satisfied: pandas>=0.23.0 in /usr/local/lib/python3.6/dist-packages (from geopandas) (1.1.4)\n",
            "Requirement already satisfied: attrs>=17 in /usr/local/lib/python3.6/dist-packages (from fiona->geopandas) (20.3.0)\n",
            "Collecting munch\n",
            "  Downloading https://files.pythonhosted.org/packages/cc/ab/85d8da5c9a45e072301beb37ad7f833cd344e04c817d97e0cc75681d248f/munch-2.5.0-py2.py3-none-any.whl\n",
            "Requirement already satisfied: six>=1.7 in /usr/local/lib/python3.6/dist-packages (from fiona->geopandas) (1.15.0)\n",
            "Collecting click-plugins>=1.0\n",
            "  Downloading https://files.pythonhosted.org/packages/e9/da/824b92d9942f4e472702488857914bdd50f73021efea15b4cad9aca8ecef/click_plugins-1.1.1-py2.py3-none-any.whl\n",
            "Requirement already satisfied: click<8,>=4.0 in /usr/local/lib/python3.6/dist-packages (from fiona->geopandas) (7.1.2)\n",
            "Collecting cligj>=0.5\n",
            "  Downloading https://files.pythonhosted.org/packages/42/1e/947eadf10d6804bf276eb8a038bd5307996dceaaa41cfd21b7a15ec62f5d/cligj-0.7.1-py3-none-any.whl\n",
            "Requirement already satisfied: certifi in /usr/local/lib/python3.6/dist-packages (from fiona->geopandas) (2020.11.8)\n",
            "Requirement already satisfied: python-dateutil>=2.7.3 in /usr/local/lib/python3.6/dist-packages (from pandas>=0.23.0->geopandas) (2.8.1)\n",
            "Requirement already satisfied: pytz>=2017.2 in /usr/local/lib/python3.6/dist-packages (from pandas>=0.23.0->geopandas) (2018.9)\n",
            "Requirement already satisfied: numpy>=1.15.4 in /usr/local/lib/python3.6/dist-packages (from pandas>=0.23.0->geopandas) (1.18.5)\n",
            "Installing collected packages: munch, click-plugins, cligj, fiona, pyproj, geopandas\n",
            "Successfully installed click-plugins-1.1.1 cligj-0.7.1 fiona-1.8.18 geopandas-0.8.1 munch-2.5.0 pyproj-3.0.0.post1\n",
            "Collecting rasterio\n",
            "\u001b[?25l  Downloading https://files.pythonhosted.org/packages/33/1a/51baddc8581ead98fcef591624b4b2521b581943a9178912a2ac576e0235/rasterio-1.1.8-1-cp36-cp36m-manylinux1_x86_64.whl (18.3MB)\n",
            "\u001b[K     |████████████████████████████████| 18.3MB 174kB/s \n",
            "\u001b[?25hRequirement already satisfied: click-plugins in /usr/local/lib/python3.6/dist-packages (from rasterio) (1.1.1)\n",
            "Requirement already satisfied: cligj>=0.5 in /usr/local/lib/python3.6/dist-packages (from rasterio) (0.7.1)\n",
            "Requirement already satisfied: click<8,>=4.0 in /usr/local/lib/python3.6/dist-packages (from rasterio) (7.1.2)\n",
            "Collecting affine\n",
            "  Downloading https://files.pythonhosted.org/packages/ac/a6/1a39a1ede71210e3ddaf623982b06ecfc5c5c03741ae659073159184cd3e/affine-2.3.0-py2.py3-none-any.whl\n",
            "Requirement already satisfied: attrs in /usr/local/lib/python3.6/dist-packages (from rasterio) (20.3.0)\n",
            "Collecting snuggs>=1.4.1\n",
            "  Downloading https://files.pythonhosted.org/packages/cc/0e/d27d6e806d6c0d1a2cfdc5d1f088e42339a0a54a09c3343f7f81ec8947ea/snuggs-1.4.7-py3-none-any.whl\n",
            "Requirement already satisfied: numpy in /usr/local/lib/python3.6/dist-packages (from rasterio) (1.18.5)\n",
            "Requirement already satisfied: pyparsing>=2.1.6 in /usr/local/lib/python3.6/dist-packages (from snuggs>=1.4.1->rasterio) (2.4.7)\n",
            "Installing collected packages: affine, snuggs, rasterio\n",
            "Successfully installed affine-2.3.0 rasterio-1.1.8 snuggs-1.4.7\n",
            "Collecting alphashape\n",
            "  Downloading https://files.pythonhosted.org/packages/b3/54/e96d3c4730eb7676b7617b2d8a2aa63b29942f6aee6c58d06e8c7661b4e9/alphashape-1.1.0-py2.py3-none-any.whl\n",
            "Requirement already satisfied: scipy>=1.0.0 in /usr/local/lib/python3.6/dist-packages (from alphashape) (1.4.1)\n",
            "Requirement already satisfied: shapely>=1.4.0 in /usr/local/lib/python3.6/dist-packages (from alphashape) (1.7.1)\n",
            "Collecting click-log>=0.3.2\n",
            "  Downloading https://files.pythonhosted.org/packages/38/52/a9dbb622f40ceeb09df141d855062cc9fbb38011f0ad8caee0cd840f399c/click_log-0.3.2-py2.py3-none-any.whl\n",
            "Requirement already satisfied: Click>=6.0 in /usr/local/lib/python3.6/dist-packages (from alphashape) (7.1.2)\n",
            "Requirement already satisfied: numpy>=1.13.3 in /usr/local/lib/python3.6/dist-packages (from scipy>=1.0.0->alphashape) (1.18.5)\n",
            "Installing collected packages: click-log, alphashape\n",
            "Successfully installed alphashape-1.1.0 click-log-0.3.2\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "fmRo_blmu_1p",
        "outputId": "61b45f2f-2c09-4422-868f-cd4b6d264672"
      },
      "source": [
        "device = torch.device(\"cuda\" if torch.cuda.is_available() \n",
        "                                  else \"cpu\")\n",
        "print(device)\n",
        "\n",
        "from google.colab import drive\n",
        "\n",
        "drive.mount('/content/drive')"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "cuda\n",
            "Mounted at /content/drive\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "n_WQhY4MPHXO"
      },
      "source": [
        "model_version = \"alex_8\"\n",
        "params_file = model_version + '_params.pth' #model version to use\n",
        "\n",
        "#get model\n",
        "model_dir = '/content/drive/My Drive/Lake detection/models/'\n",
        "\n",
        "params = torch.load(model_dir + params_file)\n",
        "model_name = params['name']\n",
        "\n",
        "preprocess = transforms.Compose([transforms.ToTensor(),\n",
        "        transforms.Normalize([0.485, 0.456, 0.406], \n",
        "                             [0.229, 0.224, 0.225])\n",
        "                                     ])"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "W89Rgn3tvPKI"
      },
      "source": [
        "def make_folder(folder): #creates a new directory if it does not exist\n",
        "    directory = os.path.dirname(folder)\n",
        "    if not os.path.exists(directory):\n",
        "        os.makedirs(directory)\n",
        "\n",
        "def re_orderbands(b1, b2, b3): #reprders image bands\n",
        "    img_new = np.zeros((b1.shape[0], b1.shape[1],3))\n",
        "    img_new[:,:,0] = b1\n",
        "    img_new[:,:,1] = b2\n",
        "    img_new[:,:,2] = b3\n",
        "    return img_new\n",
        "\n",
        "def normalize_band(data,dmin, dmax): #Normalize: min = 0, max = 1\n",
        "    data_n = np.zeros(data.shape)\n",
        "    data_n = (data - dmin)/(dmax - dmin)\n",
        "    return data_n\n",
        "\n",
        "def loadImages(file): #load S1 and S2 data, return geographic info and data\n",
        "        raster = gdal.Open(file)\n",
        "        data = raster.ReadAsArray()\n",
        "        gray = data[0,:,:] #grayscale\n",
        "        HV = data[1,:,:] #HV\n",
        "        \n",
        "        gray[gray > 1] = 1\n",
        "        HV[HV > 1] = 1\n",
        "        HV[HV < 0] = 0\n",
        "        gray[gray == 0] = np.nan\n",
        "        HV[np.isnan(gray)] = np.nan\n",
        "        HV[(np.isnan(HV)) & (~np.isnan(gray))] = np.nanmean(HV) #gets rid of nan strips in S1 data\n",
        "\n",
        "        datas = np.zeros((HV.shape[0],HV.shape[1],2))\n",
        "        datas[:,:,0] = gray  \n",
        "        datas[:,:,1] = HV    \n",
        "        \n",
        "        ulx, xres, xskew, uly, yskew, yres  = raster.GetGeoTransform()\n",
        "        lrx = ulx + (raster.RasterXSize * xres)\n",
        "        lry = uly + (raster.RasterYSize * yres)\n",
        "        x_minmax = (ulx, lrx)\n",
        "        y_minmax = (uly, lry)\n",
        "\n",
        "        #image shape\n",
        "        m = datas.shape[0]\n",
        "        n = datas.shape[1]\n",
        "\n",
        "        #x and y coordinates for each pixel\n",
        "        xx = np.linspace(x_minmax[0], x_minmax[1], n) \n",
        "        yy = np.linspace(y_minmax[0], y_minmax[1], m)\n",
        "        \n",
        "        return datas, m, n, xx, yy\n",
        "\n",
        "def get_mtns(all_data, mt):\n",
        "    mtns = np.zeros(all_data[:,:,0].shape)\n",
        "    for item in mt:\n",
        "        i = item[0]\n",
        "        j = item[1]\n",
        "        dx = item[2]\n",
        "\n",
        "        temp_img = all_data[i:i+dx, j:j+dx, 0]\n",
        "        temp_mtns = np.zeros(temp_img.shape)\n",
        "        temp_mtns[temp_img < 0.4] = 1\n",
        "\n",
        "        mtns[i:i+dx, j:j+dx] = temp_mtns\n",
        "\n",
        "    if len(mt) > 5:\n",
        "        mtns[all_data[:,:,0] < 0.4] = 1\n",
        "    #mtns = open_space(mtns, np.ones((5,5)), 1)\n",
        "    \n",
        "    mtns = cv2.dilate(mtns,np.ones((20,20)),iterations = 1)\n",
        "    return mtns\n",
        "\n",
        "\n",
        "def fill_holes(img, kernel,i): #morphological closing operation\n",
        "    closing = cv2.morphologyEx(img, cv2.MORPH_CLOSE, kernel, iterations = i)\n",
        "    return closing\n",
        "\n",
        "def open_space(img, kernel,i): #morphological opening operation\n",
        "    opening = cv2.morphologyEx(img, cv2.MORPH_OPEN, kernel, iterations = i)\n",
        "    return opening\n",
        "\n",
        "def filter_lake(img): #perform morphological operations\n",
        "    kernel1 = np.ones((3,3),np.uint8)\n",
        "    kernel2 = np.ones((6,6),np.uint8)\n",
        "    kernel3 = np.ones((12,12),np.uint8)\n",
        "\n",
        "    th2 = fill_holes(img, kernel1, 1)\n",
        "    th2 = open_space(th2, kernel1, 1)\n",
        "    th3 = fill_holes(th2, kernel2, 1)\n",
        "    th3 = open_space(th3, kernel2, 1)\n",
        "    th4 = open_space(th3, kernel3, 1) \n",
        "    th4 = fill_holes(th4, kernel3, 1) \n",
        "\n",
        "    return th4\n",
        "\n",
        "def thresholding(img):#create binary image based on threshold\n",
        "      v2 = np.percentile(img, 5)\n",
        "      new_img = np.zeros(img.shape)\n",
        "      new_img[img < v2] = 1\n",
        "      return new_img\n",
        "\n",
        "def convert_poly(co, xx, yy): #convert polygon of indexes to geo-located polygon\n",
        "    geo_contour = list()\n",
        "    for ii in range(len(co)):\n",
        "        y,x = co[ii]\n",
        "        if int(y) < 1:\n",
        "            new_y = yy[0]\n",
        "        elif int(y) >= len(yy)-1:\n",
        "            new_y = yy[-1]\n",
        "        elif int(y) == y:\n",
        "            new_y = yy[int(y)]\n",
        "        else:\n",
        "            temp1 = int(y - 0.5)\n",
        "            temp2 = int(y + 0.5)\n",
        "            new_y = (yy[temp1] + yy[temp2])/2\n",
        "\n",
        "        if int(x) < 1:\n",
        "            new_x = xx[0]\n",
        "        elif int(x) >= len(xx)-1:\n",
        "            new_x = xx[-1]\n",
        "        elif int(x) == x:\n",
        "            new_x = xx[int(x)]\n",
        "        else:\n",
        "            temp1 = int(x - 0.5)\n",
        "            temp2 = int(x + 0.5)\n",
        "            new_x = (xx[temp1] + xx[temp2])/2\n",
        "            \n",
        "        tup = new_x, new_y\n",
        "        geo_contour.append(tup)\n",
        "    return geo_contour\n",
        "    \n",
        "def buffer_poly(poly, file, b, band): #create a buffer around the polygon\n",
        "    multipoly = MultiPolygon([poly])\n",
        "    new_multipoly = MultiPolygon([Polygon(poly.buffer(b), [poly.exterior.coords])])\n",
        "    with rasterio.open(file) as src:\n",
        "        out_image1, out_transform1 = rasterio.mask.mask(src, multipoly, crop = True)\n",
        "        out_image1 = np.squeeze(out_image1)\n",
        "        out_image1[out_image1 == 0] = np.nan\n",
        "        out_image2, out_transform2 = rasterio.mask.mask(src, new_multipoly, crop = True)\n",
        "        out_image2 = np.squeeze(out_image2)\n",
        "        out_image2[out_image2 == 0] = np.nan\n",
        "    if len(out_image1.shape) > 2:\n",
        "        return out_image1[band,:,:], out_image2[band,:,:]\n",
        "    else:\n",
        "        return None, None\n",
        "\n",
        "def edge_check(ilake): #check to see if lake is too close to edge of image\n",
        "    good = False\n",
        "    maxy = int(np.max(ilake[:,0]))\n",
        "    maxx = int(np.max(ilake[:,1]))\n",
        "    miny = int(np.min(ilake[:,0]))\n",
        "    minx = int(np.min(ilake[:,1]))\n",
        "\n",
        "    r = np.min([maxy-miny, maxx-minx, 100])\n",
        "    temp_img = all_data[miny-r:maxy+r, minx-r:maxx+r,1]\n",
        "\n",
        "    if np.mean(temp_img) > -5:\n",
        "        good = True  \n",
        "\n",
        "    return good\n",
        "\n",
        "def mtn_check(ilake, properties):\n",
        "    good = True\n",
        "    l = convert_poly(ilake, xx, yy)\n",
        "    p = Polygon(l)\n",
        "    im, b = buffer_poly(p, file, 1000,1)\n",
        "    im2, b2 = buffer_poly(p, file, 1000,0)\n",
        "\n",
        "    if len(mt) > 10 and np.nanmin(b2) < 0.4: good = False\n",
        "   \n",
        "    maxy = int(np.max(ilake[:,0]))\n",
        "    maxx = int(np.max(ilake[:,1]))\n",
        "    miny = int(np.min(ilake[:,0]))\n",
        "    minx = int(np.min(ilake[:,1]))\n",
        "    r = 5\n",
        "\n",
        "    temp_img = mtns[miny-r:maxy+r, minx-r:maxx+r]\n",
        "\n",
        "    if np.mean(temp_img) > 0:\n",
        "        good = False\n",
        "\n",
        "    ran = np.nanmax(im2) - np.nanmin(im2)\n",
        "    if ran > 0.75:\n",
        "        good = False\n",
        "\n",
        "    # fig, ax = plt.subplots(1,5,figsize = (20,5))\n",
        "    # ax[0].imshow(im, cmap = 'gray', vmin = 0, vmax = 1)\n",
        "    # ax[1].imshow(b, cmap = 'gray', vmin = 0, vmax = 1)\n",
        "    # ax[2].imshow(im2, cmap = 'gray', vmin = 0, vmax = 1)\n",
        "    # ax[3].imshow(b2, cmap = 'gray', vmin = 0, vmax = 1)\n",
        "    # ax[4].imshow(temp_img, cmap = 'gray', vmin = 0, vmax = 1)\n",
        "\n",
        "    return good\n",
        "\n",
        "def get_means(im, b):\n",
        "\n",
        "    HV_mean = np.nanmean(im)\n",
        "    im_non_nans = (~np.isnan(im)).sum()\n",
        "    total_non_nans = (~np.isnan(b)).sum()\n",
        "    b_non_nans = total_non_nans - im_non_nans\n",
        "    total_mean = np.nanmean(b)\n",
        "    buffer_mean = (total_non_nans*total_mean - im_non_nans*HV_mean)/b_non_nans\n",
        "    return HV_mean, buffer_mean\n",
        "\n",
        "\n",
        "def list_to_json(lst, pth): #write list to json file\n",
        "    _json = {'type': 'FeatureCollection', 'features':lst}\n",
        "    json_string = json.dumps(_json)\n",
        "    with open(pth, 'w') as json_file:\n",
        "        json.dump(_json, json_file)\n",
        "\n",
        "def refine_lakes1(contoured_lakes):\n",
        "    full_im = np.zeros(all_data[:,:,1].shape)\n",
        "    lakes = list()\n",
        "    remaining_lakes = list()\n",
        "\n",
        "    for ilake in contoured_lakes:\n",
        "        #ilake = contoured_lakes[i]\n",
        "        l = convert_poly(ilake, xx, yy)\n",
        "        p = Polygon(l)\n",
        "        im, b = buffer_poly(p, file, 1000,1)\n",
        "\n",
        "        if im is not None:\n",
        "            HV_mean, buffer_mean = get_means(im, b)\n",
        "\n",
        "            maxy = int(np.max(ilake[:,0]))+15\n",
        "            maxx = int(np.max(ilake[:,1]))+15\n",
        "            miny = int(np.min(ilake[:,0]))-15\n",
        "            minx = int(np.min(ilake[:,1]))-15\n",
        "\n",
        "            xxtemp = xx[minx:maxx]\n",
        "            yytemp = yy[miny:maxy]\n",
        "            image = all_data[miny:maxy, minx:maxx,1]\n",
        "\n",
        "            th = np.zeros(image.shape)\n",
        "            th[(image - buffer_mean) < -0.2] = 1\n",
        "\n",
        "            kernel1 = cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(3,3))\n",
        "            kernel2 = cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(6,6))\n",
        "            th2 = fill_holes(th, kernel1,1)\n",
        "            th3 = open_space(th2, kernel1,1)\n",
        "            th4 = fill_holes(th3, kernel1,1)\n",
        "            th5 = open_space(th4, kernel1,1)\n",
        "            th5 = cv2.blur(th5,(3,3))\n",
        "            th5[th5 > 0] = 1\n",
        "            th6 = fill_holes(th5, kernel2,1)\n",
        "            th7 = open_space(th6, kernel2,1)\n",
        "            th7 = cv2.blur(th7,(6,6))\n",
        "            th7[th7 > 0] = 1\n",
        "\n",
        "            full_im[miny-2:maxy-2, minx-2:maxx-2] = th7\n",
        "\n",
        "    full_im[full_im < 1] = 0\n",
        "    full_im[full_im >= 1] = 1\n",
        "\n",
        "     #contour lakes\n",
        "    contoured_lakes = measure.find_contours(full_im, 0.5) #contour 'lakes'\n",
        "    del full_im\n",
        "\n",
        "    for co in contoured_lakes:\n",
        "        properties = get_stats(co, 0)\n",
        "        l = convert_poly(co, xx, yy)\n",
        "        p = Polygon(l)\n",
        "        p2 = affinity.scale(p, xfact = 0.9, yfact = 0.9)\n",
        "        coords = list(zip(*p2.exterior.coords.xy))\n",
        "        d = int(len(coords)/25)\n",
        "        if d > 1:\n",
        "            coords = coords[0::d]\n",
        "            coords.append(coords[0])\n",
        "\n",
        "        if check_lake(properties, 0) and mtn_check(co, properties) and edge_check(co):\n",
        "            coords = [coords]\n",
        "            if properties['area'] > 250000 and properties['cvh_ratio'] > 0.95 and properties['Q4'] < -0.15 and properties ['Q3'] < -0.17 and properties['background_diff'] < -0.15 and properties['w_l_ratio'] < 2.5: properties['code'] = 'round'\n",
        "            elif properties['area'] > 250000 and properties['cvh_ratio'] > 0.9 and properties['Q4'] < -0.16 and properties['background_diff'] < -0.16 and properties['w_l_ratio'] < 2.5: properties['code'] = 'round'\n",
        "            elif properties['area'] > 1500000 and properties['background_diff'] < -0.15 and properties['Q4'] < -0.14 and properties['Q1'] < -0.2:  properties['code'] = 'big'\n",
        "            elif properties['area'] > 1000000 and properties['background_diff'] < -0.19 and properties['Q4'] < -0.19 and properties['Q1'] < -0.2:  properties['code'] = 'big'\n",
        "            elif properties['area'] < 150000 and properties ['area'] > 50000 and properties['Q4'] < -0.2 and properties['cvh_ratio'] > 0.8 and properties['w_l_ratio'] < 2.5: properties['code'] = 'small'\n",
        "            #elif properties['w_l_ratio'] > 2.5 and properties['background_diff'] < -0.2 and properties['Q4'] < -0.2 and properties['area'] > 150000: properties['code'] = 'long'\n",
        "            elif properties['HV_mean'] < 0.3 and properties['Q4'] < -0.16 and properties['background_diff'] < -0.18 and properties['area'] > 50000 and properties['w_l_ratio'] < 2.5: properties['code'] = 'dark'\n",
        "            else: properties['code']  = 'good'\n",
        "            poly_dict = {'type':'Polygon', 'coordinates':coords}\n",
        "            geom = {'type':'Feature', 'geometry':poly_dict, 'properties':properties}\n",
        "            lakes.append(geom)\n",
        "        elif properties['area'] > 50000 and properties['background_diff'] < -0.1:\n",
        "            remaining_lakes.append(co)\n",
        "  \n",
        "    del contoured_lakes\n",
        "    return lakes, remaining_lakes\n",
        "\n",
        "def refine_lakes2(contoured_lakes):\n",
        "    full_im = np.zeros(all_data[:,:,1].shape)\n",
        "    lakes = list()\n",
        "    remaining_lakes = list()\n",
        "\n",
        "    for ilake in contoured_lakes:\n",
        "        l = convert_poly(ilake, xx, yy)\n",
        "        p = Polygon(l)\n",
        "        im, b = buffer_poly(p, file, 1000,1)\n",
        "\n",
        "        if im is not None:\n",
        "            HV_mean, buffer_mean = get_means(im, b)\n",
        "\n",
        "            maxy = int(np.max(ilake[:,0]))+15\n",
        "            maxx = int(np.max(ilake[:,1]))+15\n",
        "            miny = int(np.min(ilake[:,0]))-15\n",
        "            minx = int(np.min(ilake[:,1]))-15\n",
        "\n",
        "            xxtemp = xx[minx:maxx]\n",
        "            yytemp = yy[miny:maxy]\n",
        "            image = all_data[miny:maxy, minx:maxx,1]\n",
        "\n",
        "            th = np.zeros(image.shape)\n",
        "            th[(image - buffer_mean) < -0.2] = 1\n",
        "\n",
        "            kernel1 = cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(3,3))\n",
        "            kernel2 = cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(6,6))\n",
        "            kernel3 = cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(12,12))\n",
        "            th2 = fill_holes(th, kernel1,1)\n",
        "            th3 = open_space(th2, kernel1,1)\n",
        "            th4 = fill_holes(th3, kernel1,1)\n",
        "            th5 = open_space(th4, kernel1,1)\n",
        "            th6 = fill_holes(th5, kernel2,1)\n",
        "            th7 = open_space(th6, kernel2,1)\n",
        "            if len(xxtemp) * len(yytemp) > 30000:\n",
        "                th7 = fill_holes(th7, kernel3,1)\n",
        "                th7 = open_space(th7, kernel3,1)\n",
        "\n",
        "        full_im[miny-2:maxy-2, minx-2:maxx-2] = th7\n",
        "\n",
        "    full_im[full_im < 1] = 0\n",
        "    full_im[full_im >= 1] = 1\n",
        "\n",
        "     #contour lakes\n",
        "    contoured_lakes = measure.find_contours(full_im, 0.5) #contour 'lakes'\n",
        "    del full_im\n",
        "\n",
        "    for co in contoured_lakes:\n",
        "        properties = get_stats(co,1)\n",
        "        l = convert_poly(co, xx, yy)\n",
        "        p = Polygon(l)\n",
        "        im, b = buffer_poly(p, file, 1000, 1)\n",
        "\n",
        "        coords = list(zip(*p.exterior.coords.xy))\n",
        "        d = int(len(coords)/25)\n",
        "        if d > 1:\n",
        "            coords = coords[0::d]\n",
        "            coords.append(coords[0])\n",
        "        if check_lake(properties, 1) and mtn_check(co, properties) and edge_check(co):\n",
        "            coords = [coords]\n",
        "            if properties['area'] > 250000 and properties['cvh_ratio'] > 0.95 and properties['Q4'] < -0.15 and properties ['Q3'] < -0.17 and properties['background_diff'] < -0.15 and properties['w_l_ratio'] < 2.5: properties['code'] = '1 round'\n",
        "            elif properties['area'] > 250000 and properties['cvh_ratio'] > 0.9 and properties['Q4'] < -0.16 and properties['background_diff'] < -0.16 and properties['w_l_ratio'] < 2.5: properties['code'] = '1 round'\n",
        "            elif properties['area'] > 1500000 and properties['background_diff'] < -0.15 and properties['Q4'] < -0.14 and properties['Q1'] < -0.2:  properties['code'] = '1 big'\n",
        "            elif properties['area'] > 1000000 and properties['background_diff'] < -0.19 and properties['Q4'] < -0.19 and properties['Q1'] < -0.2:  properties['code'] = '1 big'\n",
        "            elif properties['area'] < 150000 and properties ['area'] > 50000 and properties['Q4'] < -0.2 and properties['cvh_ratio'] > 0.8 and properties['w_l_ratio'] < 2.5: properties['code'] = '1 small'\n",
        "            else: properties['code']  = '1 good'\n",
        "            poly_dict = {'type':'Polygon', 'coordinates':coords}\n",
        "            geom = {'type':'Feature', 'geometry':poly_dict, 'properties':properties}\n",
        "            lakes.append(geom)\n",
        "\n",
        "    return lakes\n",
        "\n",
        "def get_stats(ilake, code):\n",
        "    properties = {}\n",
        "\n",
        "    lake = convert_poly(ilake, xx, yy)\n",
        "    poly = Polygon(lake) \n",
        "    if code == 0:\n",
        "        poly = affinity.scale(poly, xfact = 0.9, yfact = 0.9)\n",
        "    properties['area'] = poly.area #7) area\n",
        "    convex_hull = poly.convex_hull\n",
        "    cvh_area = convex_hull.area\n",
        "    properties['cvh_ratio'] = poly.area/cvh_area #8) poly to convex hull poly area ratio\n",
        "    \n",
        "    im, b = buffer_poly(poly, file, 1000,1)\n",
        "    if im is not None:\n",
        "        HV_mean, buffer_mean = get_means(im, b)\n",
        "        properties['HV_mean'] = HV_mean.astype(float) #1) average HV value\n",
        "        properties['background_diff'] = (HV_mean - buffer_mean).astype(float) #difference of lake from background\n",
        "\n",
        "        maxy = int(np.max(ilake[:,0]))\n",
        "        maxx = int(np.max(ilake[:,1]))\n",
        "        miny = int(np.min(ilake[:,0]))\n",
        "        minx = int(np.min(ilake[:,1]))\n",
        "        rangex = maxx-minx\n",
        "        rangey = maxy-miny\n",
        "\n",
        "        Q1 = np.nanmean(all_data[miny-rangey:maxy-rangey, minx:maxx,1])\n",
        "        Q2 = np.nanmean(all_data[miny+rangey:maxy+rangey, minx:maxx,1])\n",
        "        Q3 = np.nanmean(all_data[miny:maxy, minx-rangex:maxx-rangex,1])\n",
        "        Q4 = np.nanmean(all_data[miny:maxy, minx+rangex:maxx+rangex,1])\n",
        "        Qs = [HV_mean - Q1, HV_mean - Q2, HV_mean - Q3, HV_mean - Q4]\n",
        "        Qs.sort()\n",
        "        properties['Q1'] = Qs[0].astype(float) #4-7 quadrant differences\n",
        "        properties['Q2'] = Qs[1].astype(float)\n",
        "        properties['Q3'] = Qs[2].astype(float)\n",
        "        properties['Q4'] = Qs[3].astype(float)\n",
        "\n",
        "        box = poly.minimum_rotated_rectangle\n",
        "        xbox, ybox = box.exterior.coords.xy\n",
        "        edge_length = (Point(xbox[0], ybox[0]).distance(Point(xbox[1], ybox[1])), Point(xbox[1], ybox[1]).distance(Point(xbox[2], ybox[2]))) # get length of bounding box edges\n",
        "        length = max(edge_length)\n",
        "        width = min(edge_length)\n",
        "        properties['w_l_ratio'] = length/width\n",
        "\n",
        "    else:\n",
        "        properties['HV_mean'] = 0\n",
        "        properties['background_diff'] = 0\n",
        "        properties['Q1'] = 0\n",
        "        properties['Q2'] = 0\n",
        "        properties['Q3'] = 0\n",
        "        properties['Q4'] = 0\n",
        "        properties['w_l_ratio'] = 0\n",
        "\n",
        "    return properties\n",
        "\n",
        "def check_lake(properties, t):\n",
        "    round_lake = False\n",
        "    if properties['area'] > 250000 and properties['cvh_ratio'] > 0.95 and properties['Q4'] < -0.1575 and properties ['Q3'] < -0.17 and properties['background_diff'] < -0.1575 and properties['w_l_ratio'] < 2.5: round_lake = True\n",
        "    elif properties['area'] > 250000 and properties['cvh_ratio'] > 0.9 and properties['Q4'] < -0.168 and properties['background_diff'] < -0.168 and properties['w_l_ratio'] < 2.5: round_lake = True\n",
        "    big_lake = False\n",
        "    if properties['area'] > 1500000 and properties['background_diff'] < -0.1575 and properties['Q4'] < -0.147 and properties['Q1'] < -0.2:  big_lake = True\n",
        "    elif properties['area'] > 1000000 and properties['background_diff'] < -0.1995 and properties['Q4'] < -0.1995 and properties['Q1'] < -0.2:  big_lake = True\n",
        "    small_lake = False\n",
        "    if properties['area'] < 150000 and properties ['area'] > 50000 and properties['Q4'] < -0.21 and properties['cvh_ratio'] > 0.8 and properties['w_l_ratio'] < 2.5: small_lake = True\n",
        "    dark_lake = False\n",
        "    if t == 0:\n",
        "        #if properties['w_l_ratio'] > 2.5 and properties['background_diff'] < -0.2 and properties['Q4'] < -0.2 and properties['area'] > 150000: long_lake = True\n",
        "        if properties['HV_mean'] < 0.3 and properties['Q4'] < -0.16 and properties['background_diff'] < -0.18 and properties['area'] > 50000 and properties['w_l_ratio'] < 2.5: dark_lake = True\n",
        "    #cond1 = properties['area'] > 750000 or (properties['area'] < 750000 and properties['area'] > 150000 and properties['cvh_ratio'] > 0.7)\n",
        "    cond2 = properties['Q4'] < -0.1785 and properties['Q3'] < -0.20 and properties['background_diff'] < -0.21 and properties ['area'] > 150000# and properties['w_l_ratio'] < 2.5\n",
        "    cond3 = True\n",
        "    if big_lake and properties['background_diff'] > -0.2 and (properties['Q4'] - properties['background_diff']) < -0.01: cond3 = False\n",
        "    if not big_lake and properties['background_diff'] > -0.22 and (properties['Q4'] - properties['background_diff']) < -0.01: cond3 = False\n",
        "    \n",
        "    \n",
        "    return (cond2 or round_lake or big_lake or small_lake or dark_lake) and cond3\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "90dz9nyxCl6I",
        "outputId": "51022062-c531-4190-ab0b-aae70571076f"
      },
      "source": [
        "regions = ['CW', 'SW', 'NW', 'NO', 'NE', 'SE']\n",
        "year = '2019'\n",
        "\n",
        "region = 'SE'\n",
        "if year == '2018' or year == '2019':\n",
        "    tifs_folder = '/content/drive/My Drive/Lake detection/greenland/' + year + '/' + region[0:2] + '/original_data/' #original data location\n",
        "    tiles_folder = '/content/drive/My Drive/Lake detection/greenland/' + year + '/' + region[0:2]  + '/lake_tiles/' #location to store tiles\n",
        "else:\n",
        "    tifs_folder = '/content/drive/My Drive/Lake detection/greenland/' + year + '/original_data/' #original data location\n",
        "    tiles_folder = '/content/drive/My Drive/Lake detection/greenland/' + year + '/lake_tiles/' #location to store tiles\n",
        "make_folder(tiles_folder) #make folder to store subsurface lake image tiles\n",
        "dxs = [256,512]\n",
        "\n",
        "files = sorted(glob(tifs_folder + region + '*np.tif'))\n",
        "#files = [files[-1]]\n",
        "lakes = list()\n",
        "num = 0\n",
        "\n",
        "for file in files:\n",
        "    mt = list()\n",
        "    print(file)\n",
        "    print('loading model')\n",
        "    model = torch.load(model_dir + model_name)\n",
        "    for param in model.parameters():\n",
        "        param.requires_grad = False\n",
        "    model.eval();\n",
        "\n",
        "    print(\"loading data\")\n",
        "    all_data, m, n, xx, yy = loadImages(file) #get original image, pixel dimensions (m,n) and x and y coordinates for each pixel (xx, yy)\n",
        "    full_im = np.zeros((m,n))\n",
        "\n",
        "    for dx in dxs:\n",
        "        print(f'{dx}: classifying images')\n",
        "        for k in range(0,m-dx,int(dx/2)):\n",
        "            y1 = yy[k+dx]\n",
        "            y2 = yy[k]\n",
        "            for j in range(0,n-dx,int(dx/2)):\n",
        "                f_name = region + '_' + str(num) + '_' + str(dx) + '.jpg'                    \n",
        "                x1 = xx[j]\n",
        "                x2 = xx[j+dx]\n",
        "                temp_img = all_data[k:k+dx,j:j+dx,:] #tile\n",
        "                temp_img = resize(temp_img, (256, 256),anti_aliasing=True) #resize\n",
        "\n",
        "                if np.mean(temp_img) > -5:                      \n",
        "                    HV_norm = normalize_band(temp_img[:,:,1], np.nanmin(temp_img[:,:,1]), np.nanmax(temp_img[:,:,1])) #normalize HV band\n",
        "                    img= re_orderbands(temp_img[:,:,0], HV_norm, temp_img[:,:,1]) *255 #rearrange bands\n",
        "                    img = img.astype(np.uint8) #convert values to int\n",
        "                    img_ob = Image.fromarray(img) #convert to PIL Image object\n",
        "                    img2 = preprocess(img_ob)\n",
        "                    img2 = np.expand_dims(img2, 0) #Convert 2D image to 1D vector\n",
        "                    img2 = torch.from_numpy(img2)\n",
        "\n",
        "                    #classify image\n",
        "                    inputs = Variable(img2).to(device)\n",
        "                    inputs = inputs.float()\n",
        "                    outputs = model(inputs) #labels\n",
        "                    probabilities = torch.softmax(outputs.data, 1) #get probability for each class\n",
        "\n",
        "                    if probabilities.cpu()[:,4] > 0.7: #classify as subsurface lake\n",
        "                        img_ob.save(tiles_folder + f_name) #save image tile\n",
        "                        num  = num + 1\n",
        "                        lake_tile = temp_img[:,:,1] #get HV band\n",
        "                        lake_tile = resize(lake_tile, (dx, dx),anti_aliasing=True) #resize   \n",
        "                        th_img = thresholding(lake_tile) #thresholding\n",
        "                        filtered_img = filter_lake(th_img) #filter lake\n",
        "                        full_im[k:k+dx, j:j+dx] += filtered_img\n",
        "                    elif probabilities.cpu()[:,6] > 0.75 or probabilities.cpu()[:,5] > 0.75: #there are mtns in image\n",
        "                        mt.append([k,j,dx])\n",
        "\n",
        "\n",
        "    del model\n",
        "    cuda.empty_cache()\n",
        "    gc.collect()\n",
        "        \n",
        "    full_im[full_im < 1] = 0\n",
        "    full_im[full_im >= 1] = 1\n",
        "    mtns = get_mtns(all_data, mt)\n",
        "\n",
        "    #contour lakes\n",
        "    contoured_lakes = measure.find_contours(full_im, 0.5) #contour 'lakes'\n",
        "    if len(contoured_lakes) > 0:  \n",
        "        print(f'There are potential lakes in this region')\n",
        "        good_lakes, remaining_lakes = refine_lakes1(contoured_lakes)\n",
        "        print(f'good lakes: {len(good_lakes)}, remaining: {len(remaining_lakes)}')\n",
        "        good_lakes2 = refine_lakes2(remaining_lakes)\n",
        "        print(f'new good lakes: {len(good_lakes2)}')\n",
        "        lakes = lakes + good_lakes + good_lakes2\n",
        "        print(f'final lakes: {len(lakes)}')\n",
        "        print('done')\n",
        "        del good_lakes\n",
        "        del good_lakes2\n",
        "        del remaining_lakes\n",
        "    else: print(f'There are no lakes in this region')\n",
        "\n",
        "    del full_im\n",
        "    gc.collect()\n",
        "    del mtns\n",
        "    del all_data\n",
        "  \n",
        "   \n",
        "intersection_lakes = list()\n",
        "new_lakes = list()\n",
        "for (i1, i2) in combinations(np.arange(len(lakes)), 2):  # 2 for pairs, 3 for triplets, etc\n",
        "    l1 = Polygon(lakes[i1]['geometry']['coordinates'][0])\n",
        "    l2 = Polygon(lakes[i2]['geometry']['coordinates'][0])\n",
        "    if l1.intersects(l2):\n",
        "        intersection_lakes.append(i1)\n",
        "        intersection_lakes.append(i2)\n",
        "        if l1.area < l2.area:\n",
        "            new_lakes.append(lakes[i1])\n",
        "        else:\n",
        "            new_lakes.append(lakes[i2])\n",
        "\n",
        "for i in range(len(lakes)):\n",
        "    if i not in intersection_lakes:\n",
        "        new_lakes.append(lakes[i])\n",
        "\n",
        "print(len(new_lakes))\n",
        "path = '/content/drive/My Drive/Lake detection/greenland/lakes/decrease/' + region + '_' + year + '.geojson'\n",
        "list_to_json(new_lakes, path)\n",
        "    "
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "/content/drive/My Drive/Lake detection/greenland/2019/SE/original_data/SE10_2019_np.tif\n",
            "loading model\n",
            "loading data\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:24: RuntimeWarning: invalid value encountered in greater\n",
            "/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:25: RuntimeWarning: invalid value encountered in greater\n",
            "/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:26: RuntimeWarning: invalid value encountered in less\n"
          ],
          "name": "stderr"
        },
        {
          "output_type": "stream",
          "text": [
            "256: classifying images\n",
            "512: classifying images\n",
            "There are potential lakes in this region\n",
            "good lakes: 0, remaining: 10\n",
            "new good lakes: 2\n",
            "final lakes: 2\n",
            "done\n",
            "/content/drive/My Drive/Lake detection/greenland/2019/SE/original_data/SE11_2019_np.tif\n",
            "loading model\n",
            "loading data\n",
            "256: classifying images\n",
            "512: classifying images\n",
            "There are potential lakes in this region\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:195: RuntimeWarning: Mean of empty slice\n"
          ],
          "name": "stderr"
        },
        {
          "output_type": "stream",
          "text": [
            "good lakes: 1, remaining: 13\n",
            "new good lakes: 0\n",
            "final lakes: 3\n",
            "done\n",
            "/content/drive/My Drive/Lake detection/greenland/2019/SE/original_data/SE12_2019_np.tif\n",
            "loading model\n",
            "loading data\n",
            "256: classifying images\n",
            "512: classifying images\n",
            "There are potential lakes in this region\n",
            "good lakes: 11, remaining: 18\n",
            "new good lakes: 0\n",
            "final lakes: 14\n",
            "done\n",
            "/content/drive/My Drive/Lake detection/greenland/2019/SE/original_data/SE13_2019_np.tif\n",
            "loading model\n",
            "loading data\n",
            "256: classifying images\n",
            "512: classifying images\n",
            "There are potential lakes in this region\n",
            "good lakes: 2, remaining: 4\n",
            "new good lakes: 0\n",
            "final lakes: 16\n",
            "done\n",
            "/content/drive/My Drive/Lake detection/greenland/2019/SE/original_data/SE14_2019_np.tif\n",
            "loading model\n",
            "loading data\n",
            "256: classifying images\n",
            "512: classifying images\n",
            "There are potential lakes in this region\n",
            "good lakes: 5, remaining: 6\n",
            "new good lakes: 1\n",
            "final lakes: 22\n",
            "done\n",
            "/content/drive/My Drive/Lake detection/greenland/2019/SE/original_data/SE15_2019_np.tif\n",
            "loading model\n",
            "loading data\n",
            "256: classifying images\n",
            "512: classifying images\n",
            "There are potential lakes in this region\n",
            "good lakes: 9, remaining: 4\n",
            "new good lakes: 0\n",
            "final lakes: 31\n",
            "done\n",
            "/content/drive/My Drive/Lake detection/greenland/2019/SE/original_data/SE16_2019_np.tif\n",
            "loading model\n",
            "loading data\n",
            "256: classifying images\n",
            "512: classifying images\n",
            "There are potential lakes in this region\n",
            "good lakes: 1, remaining: 1\n",
            "new good lakes: 1\n",
            "final lakes: 33\n",
            "done\n",
            "/content/drive/My Drive/Lake detection/greenland/2019/SE/original_data/SE17_2019_np.tif\n",
            "loading model\n",
            "loading data\n",
            "256: classifying images\n",
            "512: classifying images\n",
            "There are potential lakes in this region\n",
            "good lakes: 1, remaining: 0\n",
            "new good lakes: 0\n",
            "final lakes: 34\n",
            "done\n",
            "/content/drive/My Drive/Lake detection/greenland/2019/SE/original_data/SE18_2019_np.tif\n",
            "loading model\n",
            "loading data\n",
            "256: classifying images\n",
            "512: classifying images\n",
            "There are potential lakes in this region\n",
            "good lakes: 0, remaining: 1\n",
            "new good lakes: 0\n",
            "final lakes: 34\n",
            "done\n",
            "/content/drive/My Drive/Lake detection/greenland/2019/SE/original_data/SE19_2019_np.tif\n",
            "loading model\n",
            "loading data\n",
            "256: classifying images\n",
            "512: classifying images\n",
            "There are potential lakes in this region\n",
            "good lakes: 1, remaining: 0\n",
            "new good lakes: 0\n",
            "final lakes: 35\n",
            "done\n",
            "/content/drive/My Drive/Lake detection/greenland/2019/SE/original_data/SE1_2019_np.tif\n",
            "loading model\n",
            "loading data\n",
            "256: classifying images\n",
            "512: classifying images\n",
            "There are potential lakes in this region\n",
            "good lakes: 3, remaining: 21\n",
            "new good lakes: 1\n",
            "final lakes: 39\n",
            "done\n",
            "/content/drive/My Drive/Lake detection/greenland/2019/SE/original_data/SE20_2019_np.tif\n",
            "loading model\n",
            "loading data\n",
            "256: classifying images\n",
            "512: classifying images\n",
            "There are potential lakes in this region\n",
            "good lakes: 2, remaining: 8\n",
            "new good lakes: 1\n",
            "final lakes: 42\n",
            "done\n",
            "/content/drive/My Drive/Lake detection/greenland/2019/SE/original_data/SE21_2019_np.tif\n",
            "loading model\n",
            "loading data\n",
            "256: classifying images\n",
            "512: classifying images\n",
            "There are potential lakes in this region\n",
            "good lakes: 3, remaining: 39\n",
            "new good lakes: 1\n",
            "final lakes: 46\n",
            "done\n",
            "/content/drive/My Drive/Lake detection/greenland/2019/SE/original_data/SE2_2019_np.tif\n",
            "loading model\n",
            "loading data\n",
            "256: classifying images\n",
            "512: classifying images\n",
            "There are no lakes in this region\n",
            "/content/drive/My Drive/Lake detection/greenland/2019/SE/original_data/SE3_2019_np.tif\n",
            "loading model\n",
            "loading data\n",
            "256: classifying images\n",
            "512: classifying images\n",
            "There are potential lakes in this region\n",
            "good lakes: 5, remaining: 25\n",
            "new good lakes: 4\n",
            "final lakes: 55\n",
            "done\n",
            "/content/drive/My Drive/Lake detection/greenland/2019/SE/original_data/SE4_2019_np.tif\n",
            "loading model\n",
            "loading data\n",
            "256: classifying images\n",
            "512: classifying images\n",
            "There are potential lakes in this region\n",
            "good lakes: 4, remaining: 15\n",
            "new good lakes: 0\n",
            "final lakes: 59\n",
            "done\n",
            "/content/drive/My Drive/Lake detection/greenland/2019/SE/original_data/SE5_2019_np.tif\n",
            "loading model\n",
            "loading data\n",
            "256: classifying images\n",
            "512: classifying images\n",
            "There are potential lakes in this region\n",
            "good lakes: 5, remaining: 19\n",
            "new good lakes: 0\n",
            "final lakes: 64\n",
            "done\n",
            "/content/drive/My Drive/Lake detection/greenland/2019/SE/original_data/SE6_2019_np.tif\n",
            "loading model\n",
            "loading data\n",
            "256: classifying images\n",
            "512: classifying images\n",
            "There are potential lakes in this region\n",
            "good lakes: 0, remaining: 10\n",
            "new good lakes: 0\n",
            "final lakes: 64\n",
            "done\n",
            "/content/drive/My Drive/Lake detection/greenland/2019/SE/original_data/SE7_2019_np.tif\n",
            "loading model\n",
            "loading data\n",
            "256: classifying images\n",
            "512: classifying images\n",
            "There are potential lakes in this region\n",
            "good lakes: 1, remaining: 18\n",
            "new good lakes: 0\n",
            "final lakes: 65\n",
            "done\n",
            "/content/drive/My Drive/Lake detection/greenland/2019/SE/original_data/SE8_2019_np.tif\n",
            "loading model\n",
            "loading data\n",
            "256: classifying images\n",
            "512: classifying images\n",
            "There are potential lakes in this region\n",
            "good lakes: 3, remaining: 2\n",
            "new good lakes: 0\n",
            "final lakes: 68\n",
            "done\n",
            "/content/drive/My Drive/Lake detection/greenland/2019/SE/original_data/SE9_2019_np.tif\n",
            "loading model\n",
            "loading data\n",
            "256: classifying images\n",
            "512: classifying images\n",
            "There are potential lakes in this region\n",
            "good lakes: 1, remaining: 4\n",
            "new good lakes: 0\n",
            "final lakes: 69\n",
            "done\n",
            "69\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "BTAsGpuTq7bh",
        "outputId": "e00b2484-adc4-4f42-9983-8b09cc86e3a5"
      },
      "source": [
        "areas = 0\n",
        "for lake in new_lakes:\n",
        "    areas = areas + lake['properties']['area']\n",
        "print(areas/(1000**2))"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "46.04304289858684\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "z4N3x5JjE-el"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}