testGEEpull.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ee\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import shapefile\n",
    "import pandas as pd\n",
    "#ee.Authenticate(force=True)\n",
    "ee.Authenticate()\n",
    "#ee.Initialize(project = 'ee-fortschthomas52')\n",
    "ee.Initialize(project = 'ee-davidmkahler-limpopo')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Determine Area of Interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mamba, KNP\n",
    "site = \"Mamba\"\n",
    "sf = shapefile.Reader(\"/Volumes/dmk/gis/limpopo/kruger/logger_sites/reference_polygons/mamba/mamba.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Balule, KNP\n",
    "site = \"Balule\"\n",
    "sf = shapefile.Reader(\"/Volumes/dmk/gis/limpopo/kruger/logger_sites/reference_polygons/balule/balule.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Paul Kruger Gate, KNP\n",
    "site = \"Kruger Gate\"\n",
    "sf = shapefile.Reader(\"/Volumes/dmk/gis/limpopo/kruger/logger_sites/reference_polygons/kruger_gate/kruger_gate.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lower Sabie, KNP\n",
    "site = \"Lower Sabie\"\n",
    "sf = shapefile.Reader(\"/Volumes/dmk/gis/limpopo/kruger/logger_sites/reference_polygons/lower_sabie/lower_sabie.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hanover, VA, USA\n",
    "sf = shapefile.Reader(\"/Volumes/dmk/gis/limpopo/kruger/logger_sites/reference_polygons/virginia/hanover.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "shapes = sf.shapes()\n",
    "points = shapes[0].points\n",
    "aoi = ee.Geometry.Polygon(points)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Debug Code\n",
    "This code works one step at a time and produces an RGB image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "startDate = \"2023-10-10\"\n",
    "endDate = \"2023-10-11\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pull image\n",
    "image = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')\\\n",
    "        .filterDate(startDate, endDate)\\\n",
    "        .select('B2', 'B3', 'B4', 'B8', 'B11')\n",
    "         #.filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 10)\\\n",
    "mosaic = image.median().reproject(crs='EPSG:32736', scale=10) # This allows us to set the resolution.\n",
    "band_arrs = mosaic.sampleRectangle(region=aoi, defaultValue = 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Determine Dates of Acquisition\n",
    "Check out the [GEE Code Editor](https://code.earthengine.google.com/2bc6d0b39051c230c556370381c32c6f?noload=true) to determine the date of acquisition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get individual band arrays.\n",
    "band_arr_b2 = band_arrs.get('B2')   # Blue\n",
    "band_arr_b3 = band_arrs.get('B3')   # Green\n",
    "band_arr_b4 = band_arrs.get('B4')   # Red\n",
    "band_arr_b8 = band_arrs.get('B8')   # NIR\n",
    "band_arr_b11 = band_arrs.get('B11') # SWIR\n",
    "# Transfer the arrays from server to client and cast as np array.\n",
    "b2 = np.array(band_arr_b2.getInfo())   # b2  Blue\n",
    "b3 = np.array(band_arr_b3.getInfo())   # b3  Green\n",
    "b4 = np.array(band_arr_b4.getInfo())   # b4  Red\n",
    "b8 = np.array(band_arr_b8.getInfo())   # b8  NIR\n",
    "b11 = np.array(band_arr_b11.getInfo()) # b11 SWIR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWgAAAGeCAYAAAC5LpvRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAupUlEQVR4nO3db4wV933v8c9vZs45hnhZxU29a2SSUoWkTZEt/0kRaRpoUlailWXXT6o6ynVVXSkOthWuH7jFPMimUlmbB4hIOG7d5jq+qhCRbuM0DxoLpNrQCFkC18jIlixVoslW8RalcoFg2N2Z+d0HmHOzZs/vs8sf8zO8X9I+4Mycmd/5zZzvmWU/8z0hxhgFAMhOcbUHAACYHwUaADJFgQaATFGgASBTFGgAyBQFGgAyRYEGgExRoAEgUxRoAMhUdbUH8H5t2+pnP/uZhoaGFEK42sMBgMsuxqhTp05p+fLlKorEdXK8Qp5++un4a7/2a7HX68U777wzHjhwYEHPm5ycjJL44Ycffq75n8nJyWQ9vCJX0N/73ve0efNmffvb39bv/M7v6G/+5m+0ceNGvfnmm/r4xz+efO7Q0JAk6f8887SWLlky7zpxNqYHUNfJxdOz6eVn2ja5vK3S0xbq9JV/UHr7sUkuVlWlV2hjevuSdLZxv52k//erMM+PVfoYBTPHwfxyFzvJxSpmzfIifQ7UTXq5gvnfwdIM0Px2WLlj6A6fO4fN04vSrBAWUDpCeiNFNO8TO0iz+Ib0MYhmFiozB+4cazX4PXDmzBl99X892q93g4QYL3+zpDVr1ujOO+/UM88803/sN3/zN3XfffdpYmIi+dyTJ09qeHhY//e7/1tLly6ddx1boGfTM3epBbq5xAJdZFCgz1Cgk8tzL9Duf//ih6BAlzE9h6E17/PCvM9cgTafALZAz6SXpwr0u2fe1f946H/qxIkTWrZs2eB9pHexeDMzM3r11Vc1NjY25/GxsTEdPHjwgvWnp6d18uTJOT8AgCtQoH/+85+raRqNjIzMeXxkZERTU1MXrD8xMaHh4eH+z4oVKy73kADgQ+mKxezen8CIMc6bytiyZYtOnDjR/5mcnLxSQwKAD5XL/kfCj33sYyrL8oKr5ePHj19wVS1JvV5PvV7vcg8DAD70LnuB7na7uuuuu7Rv3z790R/9Uf/xffv26d57713k0OYfXmX+QNKU6f/8L2L6D0DVTPoXi3Ym/Ue6Wum/LlSl+SugeX0zM+av4+4PWJLUpP8A09g/NJpBtukxRPMHInXSyztufPbMNn+km03P8aw5h6I9B9L7b9r0ORLM+IsiPQFVYf6IGFx2wN+jULgUh9uHO43N82OdXl64PzS7P6Sn8suSqjD4GFbmD5z99Ra01iI99thj+spXvqK7775ba9eu1bPPPquf/vSneuihh67E7gDgmnRFCvQf//Ef67/+67/0l3/5l3r77be1evVq/dM//ZM+8YlPXIndAcA16Yrd6r1p0yZt2rTpSm0eAK55NEsCgExRoAEgUxRoAMgUBRoAMpVdP+jzgtqBTYVixzTqMU36im563z29m15h2uQ3Z0z+sjEZXtPlqzUZ3Na2wpFMvyW5PjWu0UzRmCy56fZTzqa3X5oMajAZ2dlEIxtJKsv08o7JGdsUuZkfmWMsk7OOpet2lN5/6Rp6mZy4JAVzP0JjNlGahkzufgKb5HY5Z5NVtw2z0htf0FpcQQNApijQAJApCjQAZIoCDQCZokADQKYo0ACQKQo0AGQq2xx0O1uomZ3/86MwGdWmdb1WTQaxk24U21E6ozrfN8fMYcbXuj6zTfobUesF9Oo136mqrkuR2qbcpp+xacZbmH65LmfcmH7XrTtHzDeGum+kNt/ZKplzyG2gVnp+K9OPuzK9kKOZn7Ljs/au3XNtcs42ZmxzzCYrbzYfzHsg2vsRUsv8/ElcQQNAtijQAJApCjQAZIoCDQCZokADQKYo0ACQKQo0AGQq2xx0Eys1A/o6N7MuxXtpfVzb2Es/v0hnZENlMo7B9bF1/a7T42sLn7Hs2Ha06W20Zg4KF4J1IVTTM9vlSG0vYHtpYvotm5x10aQzsrU5RiG6DK6ZwFlzDtnXf4nHb0HrpOcoxnTj9tbuIL39ski/D9uQDovPmnM0lSWfqd0Zeg5X0ACQKQo0AGSKAg0AmaJAA0CmKNAAkCkKNABkigINAJnKNgfdqlUb588ZtiZH7FOwJgftYsRmeWtCto0JobpevEWYSS63EWdJjenF25qcZmVyzsH1i3bHKJpXYXLELgru+kEXpen3bHqSV9Pp/dusu+ll3DPz15p+0XXr5tecowu4trPnockRh9Jk4StzDpmsem16jrv3ece8wtTTqwVeG3MFDQCZokADQKYo0ACQKQo0AGSKAg0AmaJAA0CmKNAAkKl8c9BtrbYd0M/V5BMLk+F0vYSLeGm9fKPrxWsysG10OW/z/Fn/uVu7OTAZ0eA+24OZI5dBNVnxwuZ4Tc9w10/ZnEOdxuSEO2b+2vT4YpN+fmUywo2Z/2DOsdJkfEOR7pUs+au/sk2Xn+jO0dZk8U3P8linX0OI6fsNmjI9/pA6h9z5f361Ba0FAPjAUaABIFMUaADIFAUaADJFgQaATFGgASBTFGgAyFS2OWgV7cCmvoULQpfpjGc0GVsTUZXrdBua9AaCGX9rMrau33RhMrKSVJmcrGI6I9o06TkOLmfcmhyxyekWruW3eX2lWV6b8bkctc2Rp6P2/sqpSr91u6Zfd2N6JRem53qQz0G7HLLbh1veFt308839Ahp0n8V7GjOHMll1JXuGuxP4HK6gASBTFGgAyBQFGgAyRYEGgExRoAEgUxRoAMgUBRoAMpVtDrqroG6YP2dYmwxha15WOWC7/edH13A6vViVyfiads+VyTmHxvQCTsdDJUlFa7LY7dn08mAyqKZXr6r08miOUTDnQGFy3DI9v7syWXpzjEsT1C5MDrttTNY+mF7Grl2262lu3gKty7lLKtz71PTMLtv0OVaant3RnIOlO4fMG70152hqjisb5D8/BgBAlijQAJApCjQAZIoCDQCZokADQKYo0ACQKQo0AGQq2xy0FN77uVDp+ry6jKnp9etaJZcmI9uajLFMzjm6fKfJUJYugywphvShD2U6JxxM0+za9axuXVDXHGOTo1ZtGi6X5tqkMDlq06vYPd+98WZb8/pNT3OXtS9cz3TXz9oFpSVF0y+56posd+1yzu59nD4H2p7JMZs5rtz7LLHc3idwfr0FrQUA+MBRoAEgUxRoAMgUBRoAMkWBBoBMUaABIFMUaADIVLY56DaEgf1Wg8mwluZzp2OCzq6fcm16DRcmIxxNjnvG9OqV62O7gI9dExFVYzKgtmGwGYTr5+y6Dde1OXVdRlbpCShNTDXahskuTJ9+hXb6NW2Wp3PmpTmHQ2nuJahNjl1SEU1P6xkzR1V6eajNWWLmuC3T59CgfvR9LkpfD57jatYE1Re2iwsdOHBA99xzj5YvX64Qgn7wgx/MWR5j1Pj4uJYvX64lS5Zo/fr1euONNxa7GwC47i26QJ8+fVq33367du3aNe/y7du3a8eOHdq1a5cOHTqk0dFRbdiwQadOnbrkwQLA9WTR/8WxceNGbdy4cd5lMUbt3LlTW7du1f333y9Jev755zUyMqLdu3frq1/96qWNFgCuI5f1j4THjh3T1NSUxsbG+o/1ej2tW7dOBw8enPc509PTOnny5JwfAMBlLtBTU1OSpJGRkTmPj4yM9Je938TEhIaHh/s/K1asuJxDAoAPrSsSswvv++tnjPGCx87bsmWLTpw40f+ZnJy8EkMCgA+dyxqzGx0dlXTuSvqWW27pP378+PELrqrP6/V66vV6l3MYAHBNuKwFeuXKlRodHdW+fft0xx13SJJmZma0f/9+PfXUU4vaVl2Uqov5c5SFixCa/GUT0jnlsjK/WLTpaWtDOkQb65n0/k08tDG/+DQL6NVre1a7X65MRLRtTI5W6eWuF2+nTc+haXet2vXjjSYnXaXPsdacg6U7RuYcKswBKEwO2h2/0pzj0YxP8llq15Lb9WV3OWmZOeiYnuSlawnu9l4NnsMqsWzOegta65f84he/0L/927/1/33s2DEdOXJEN910kz7+8Y9r8+bN2rZtm1atWqVVq1Zp27ZtWrp0qR544IHF7goArmuLLtCHDx/W7/3e7/X//dhjj0mSHnzwQX33u9/V448/rjNnzmjTpk165513tGbNGu3du1dDQ0OXb9QAcB1YdIFev369Yhz8600IQePj4xofH7+UcQHAdY9mSQCQKQo0AGSKAg0AmaJAA0CmKNAAkKl8G/Y3M2oHNGVvTcC8kLmTpWsi5nX6c6s1zeajTDN0E7AvzE0UzQKbfae4e3Ea09G/bc0YGnejhJljs/+mkz6GXXeji7nPwh7j2hwj86ULoTGv392lUZgvjUg/235jQ2NKg3kLSpK67ksbCne3k/tiCjcIs/9Oeo6jeZOUTXqOisQ5ULi70c6vt6C1AAAfOAo0AGSKAg0AmaJAA0CmKNAAkCkKNABkigINAJnKNgcd2lJhQNgymJxzUZiQq8mg1umIqP3GgDaml3cS3QAlyfXSN73iFUw+U5Jq03C9bdNfaqDGzHEnPchuk87AtubULEyOOpovZWjMFwb4byRIbz+Yc8geomjmpzE55jb9/Krspp9fpo9v5b4QQJLMed4M+Bq884Jpal90zZcmNOZ9anLOS0yW3sWwQ+IbCYK7EeH8Pha0FgDgA0eBBoBMUaABIFMUaADIFAUaADJFgQaATFGgASBT2eagZ4tGswPyxl3TJ7Z1Hzuml3HpMrZm+63r9Wt6BbsMbmn66LbRBbklFxXvddNjnDU528LkpAuTw3YZWjXp8bXBZGBNP+fW5XzNORDN+NxJGirz+s052pgceBHTr690GeDKn2PdRA5YklrzGmTOc/c2D6YneTTv09L0bA5dk9NuB/fsLlxQ/vx6C1oLAPCBo0ADQKYo0ACQKQo0AGSKAg0AmaJAA0CmKNAAkKlsc9DTM62KATnExvUKNv2iy8JkPIPJ2CaXSsHlqE2OO5qMcDD50GoBvWbb2fQ2pk2/4VYm52yW1yGdUQ3m2sHlfG1WPrlUKszrr+v0/NVxOr0D8/pnp9Pb73UHZ2wlSR2To1b69XXDkuTyojD7l2xWuzD3A0Qzxmiaapc9k3M2+w+l2b7Lyqdy4CYjfh5X0ACQKQo0AGSKAg0AmaJAA0CmKNAAkCkKNABkigINAJnKNgcdz9aKA7KuZ6t0/rBrUq4xpDOcLsPr+jXL5LBb1wt4QB/s8zp1+rC10Y1PmonpHHFr+gVXJosdzRw0Zoihkz6G3TY9vsb0e5bp9duac6iu06/PJa27mkkuPzOTHn/X5KAL1wu56iaXN+7azQXJJbWFKy/pnHNRun7QJqds3melGZ/tN216ZocwePyFrSELGwMA4CqhQANApijQAJApCjQAZIoCDQCZokADQKYo0ACQqWxz0LUa1YOytLPvpp87m/7cqbsuo2qmZdaEQE2GtzD5yxDTy2fM7guTsZUkE+G0Kc1g5qgN6Zxw4Za7Y+AGaPohhyY9R7XZQVWmtz/rctpm/93GHOSYzjGXZXr/RWH6QRemn/VCru1MFL0y22hNz2zXU9kl1UvTN92dYjGaF5i4VyAmMtJzNrGgtQAAHzgKNABkigINAJmiQANApijQAJApCjQAZIoCDQCZyjYHrVCf+5lHYfq8tiaAWbSm16/N8JpeuyaHHU0f3MZkJENIH7Ze6w9riOkcbqjNHJt+yoXS/YoHHdvzapMxrdwxdv2q6/QxLE0KtrwhfY6UremHPZ2e3455azZNen7rMr18SZneflWY49fxOd4QTRLZ5P0Ld/1oNl8uSR/jEE3WPJh+0q5tfOIc9j3nz68HAMgSBRoAMkWBBoBMUaABIFMUaADIFAUaADJFgQaATGWbg+4WQd0BWdbWZGRNu2VVMjlo16fW9OotzfhCx+VD0xnTwuU3TQZXkhozB+qanGtjXmMwOU+TMZVchjX9Gs2rU2Hm2IVclxY9s4fZ9NJu+vXPuHM0prcfoskxV6aft8mRV277khpzP4GCy9Knl1emfIU63TPbCZ30OdgW6fHFZvAcNollv4wraADIFAUaADJFgQaATFGgASBTFGgAyBQFGgAyRYEGgExlm4NuY6V2QL/YskrnD0Npev226V7ELuNbm37S0fRKLmM6n1mZ8UWTL3UZWkmqTUy5Mj2rZfoFt66fcmv6PbucsstJF2b7HZPjNhHabpVeoZ42PbsrkxEu0v26e63rt+16GZt7Ccz0t9GcH5Jad/3nblgw6sL11DYnuTlH0klzqWrT52CR6gdt5r+/3oLWes/ExIQ++9nPamhoSDfffLPuu+8+vfXWW3PWiTFqfHxcy5cv15IlS7R+/Xq98cYbi9kNAECLLND79+/Xww8/rFdeeUX79u1TXdcaGxvT6dOn++ts375dO3bs0K5du3To0CGNjo5qw4YNOnXq1GUfPABcyxb1O8aLL74459/PPfecbr75Zr366qv6whe+oBijdu7cqa1bt+r++++XJD3//PMaGRnR7t279dWvfvWCbU5PT2t6err/75MnT17M6wCAa84l/ZHwxIkTkqSbbrpJknTs2DFNTU1pbGysv06v19O6det08ODBebcxMTGh4eHh/s+KFSsuZUgAcM246AIdY9Rjjz2mz3/+81q9erUkaWpqSpI0MjIyZ92RkZH+svfbsmWLTpw40f+ZnJy82CEBwDXlov+M+sgjj+j111/Xj3/84wuWhfd9K3WM8YLHzuv1eur1XGcwALj+XNQV9KOPPqof/vCHeumll3Trrbf2Hx8dHZWkC66Wjx8/fsFVNQAgbVFX0DFGPfroo3rhhRf08ssva+XKlXOWr1y5UqOjo9q3b5/uuOMOSdLMzIz279+vp556alEDmymiygE5xW5rcsQmA9qYkOeAi/0+l+F1EcfW5JQL0ydXJkI7006nV5DUhvRvLbONydGaj3abQa3SOdrCthJOb780Oe2yMj23zTlSVOkMrOuW3NSm33LH/FZpMsSF6VVcmAPYujeBu5dAUnBh8pjOejfNDcnlZTBZeHOOBPN80xJcZZ0+SZP3QwSXsj5nUQX64Ycf1u7du/WP//iPGhoa6l8pDw8Pa8mSJQohaPPmzdq2bZtWrVqlVatWadu2bVq6dKkeeOCBxewKAK57iyrQzzzzjCRp/fr1cx5/7rnn9Kd/+qeSpMcff1xnzpzRpk2b9M4772jNmjXau3evhoaGLsuAAeB6sej/4nBCCBofH9f4+PjFjgkAIJolAUC2KNAAkCkKNABkigINAJnKth+0ZhqpHJAz7KUzmoP6SPeXN+mMaGVyyqYdtKLpU+s/FdNrxJjOXxZFOt8pSe2s2YbbQNeN0W3A5IzNMQgm5+v6IbtTvzD9qGdNL+BojoFrhewyuK3LAFcmx257hqeXT9f+LC676Tx+adLi5iXYnt+tOUZd8zZxOWl10u+hTuIU6vi3qCSuoAEgWxRoAMgUBRoAMkWBBoBMUaABIFMUaADIFAUaADKVbQ660OBPj3YmndGMJiNbula3VTpf2bjPtdbswOSkXb/oGZODPtv4plYuS+36PVemX3NrJrkwOejGZVDLdK9ht/3CbL82QeRojnE0/bBLcw4E0/TbxcAr8/oKcw63s+lezaF2OWpJhekHXaXLT5nqpyxptk2/ho6ZpGjC6O4c7yxJL1ccPL7CFaHz6y1oLQDAB44CDQCZokADQKYo0ACQKQo0AGSKAg0AmaJAA0Cmss1BNzGqaQdkSUM6Q1h1TMY0nE0uj6bXbTD5yNKFVE3OuQq1eX56APGsz6ialyATMVXdTR8DNwcup9yahtLloHPjPcHlsF2vYTc+szyYHLTLScfWHSGTdQ/pkG4wPdFdu+iwgMpRmX7NnZ4ZQkwPomuWlybn7M6xujLnsHmTpKZ4tllYQ2iuoAEgUxRoAMgUBRoAMkWBBoBMUaABIFMUaADIFAUaADKVbQ66VFA5IGtaFSajWs8ml0cT4rzBtGqNnXROOSb6wEpSa0K4tTksVdf1e0738pWkaPr5Nk16G3WRnoOyTec8Q0j3CrbdcoPJQZuMa9uarLnJ2rsBBvP8wvSTbt054FqSmyBzMFn7tpPefmEyxpJU2J7d6W0U5hzrlCZLXLksvLmfIrics3kfppab8/M8rqABIFMUaADIFAUaADJFgQaATFGgASBTFGgAyBQFGgAylW0OWtLArGk0GdZg8o1lme6V27h+zi7CaDK6Lr4ZmvTnZhHSG+h1b0jvQFI0OdfZ2XTPbJmMaFuZF2l6Bbugse2WbHr1FuYYuxyzYnoEIaZfv3v1ZSf9/K7J8reuabk5R13Et3JBaUllZcqLuzw0T4+2J7eZg5jeQeOy7Omtq0hk0Qt7BixsHwCAq4QCDQCZokADQKYo0ACQKQo0AGSKAg0AmaJAA0Cmss1B14qqB2QFWxNvdK1qY0z3Oi7M9ksTky476RWKupdc7jK+MhnYWC7gczeYMZpRlHV6HzaHbHK6rcmouhy2y1m35hiGxuSYq3QW3+asTU7a5ZQb0++5jK5htcnSmwxwSJ/C58bgqktrcsjmjVya+wWi2X7H9FUPZg5L17Q8UUiCef/1N7GgtQAAHzgKNABkigINAJmiQANApijQAJApCjQAZIoCDQCZyjYH3caodkBT2uAysC5iaHLQKrvJxa7fdGzTGdPGZFxDZTKwJuds85mS2jo9SU1t+jFHk6M2GdXQTfcTDo3rt2x6gsd0UDeWZvvtdPr5Nqft+iWb+YuXdg64LH23SJ/jMb3YxajfY7LktmW4ybKbvVemqXU0Pb3Lxjy/Y85xJSYxzCafex5X0ACQKQo0AGSKAg0AmaJAA0CmKNAAkCkKNABkigINAJnKNgcdVSgO+PwwKWaZVr3qVibHbHLEJl4pE69U6KRXqFo3vvRhCyajLMkGZSvz0W1i1GpNhlSN6TecypBKUmH6MZuD4PpBuwxuSMekFZemc65VkZ5gO3+z5iQ3PclNzF2aNb2WG186mq7rp5zeR9ccpGjuR3DXn6Vbbk5Bdw41iXPIHd7zuIIGgExRoAEgUxRoAMgUBRoAMkWBBoBMUaABIFMUaADIVL456HZmYF/lyoQ4W9Nr2LQalktX1kU6xFi1JqPbmoCl6fUb6/T2TUL2vX2kFxeF6Wd81uylNDuozUEo0yHT1mRoq9LkmF2G1pxjZTcdhG6a9PxV5q3Xmpx3LNLzF83+XQ66NSHfbrrd9rltmJbH7n6DsmPm0LXklusnbXqe277y6cV14oaJ2uTsz1vUFfQzzzyj2267TcuWLdOyZcu0du1a/ehHP+ovjzFqfHxcy5cv15IlS7R+/Xq98cYbi9kFAOA9iyrQt956q5588kkdPnxYhw8f1he/+EXde++9/SK8fft27dixQ7t27dKhQ4c0OjqqDRs26NSpU1dk8ABwLVtUgb7nnnv0B3/wB/rUpz6lT33qU/qrv/or3XjjjXrllVcUY9TOnTu1detW3X///Vq9erWef/55vfvuu9q9e/eVGj8AXLMu+o+ETdNoz549On36tNauXatjx45pampKY2Nj/XV6vZ7WrVungwcPDtzO9PS0Tp48OecHAHARBfro0aO68cYb1ev19NBDD+mFF17QZz7zGU1NTUmSRkZG5qw/MjLSXzafiYkJDQ8P939WrFix2CEBwDVp0QX605/+tI4cOaJXXnlFX/va1/Tggw/qzTff7C8PYe5fRmOMFzz2y7Zs2aITJ070fyYnJxc7JAC4Ji06ZtftdvXJT35SknT33Xfr0KFD+ta3vqU///M/lyRNTU3plltu6a9//PjxC66qf1mv11Ovt4DMDgBcZy45Bx1j1PT0tFauXKnR0VHt27dPd9xxhyRpZmZG+/fv11NPPbX47dZBcUBYszEX/h2lA4wuH9m0pl+zCXBG1+vXjL+YNX1wzVFrzf7PSedsm9n0TsrShVgv7RgopDOwwWTFXYQ1mBh2dDlw0zC7iukceJ34rVKSChNUj65ft8uZm57hhSkNtQ0JS1XHDNKcQk1temqbrL3tF22yyI05BxvTmb5InOPFwu5WWFyBfuKJJ7Rx40atWLFCp06d0p49e/Tyyy/rxRdfVAhBmzdv1rZt27Rq1SqtWrVK27Zt09KlS/XAAw8sZjcAAC2yQP/nf/6nvvKVr+jtt9/W8PCwbrvtNr344ovasGGDJOnxxx/XmTNntGnTJr3zzjtas2aN9u7dq6GhoSsyeAC4li2qQH/nO99JLg8haHx8XOPj45cyJgCAaJYEANmiQANApijQAJApCjQAZCrbftBq63M/83AZ0KY1K5TpG2OCzVGbnLLp99yWptevydh2TD5T0WRgJSmkxxh66W0Us6bXrsmYtoWZI5fTNce4KNLPr00/6dJdu9isu+lZ7nLKRTrjW5qsfjAx5XY2fQ6595gan+MNJgtfmjksi3R5imYONKCffH9xJ728MD2xg7khoUn0hW8a13X+vTEsaC0AwAeOAg0AmaJAA0CmKNAAkCkKNABkigINAJmiQANAprLNQRfVuZ95l5mMp80vmg2UwfTCNSHTaHLIZWMyyMmlUqt0n9y28hnLngm6uqx37foFuyG4fsJmfEXH5YzTz29MhDZ00iu4fsjmFFTZmH7Xmk4ujyYjrJDO+Jp21CqLdM65cPuXFFzfc9MSvDT9nEv3TnF92xuT1XdhciMVk3Y93c/jChoAMkWBBoBMUaABIFMUaADIFAUaADJFgQaATFGgASBT2eagQ1sptPMPrzI55mBywG2bzj9GzZjBpbdfKL19N+nB9EK2QekFhCwbtw9jxmXRXRbbBJHbKj2HM6fT1xZVx/RrNjnf1vSbrk075E7henKn56ep0jtoCpejdvt3/aTTrz9W/hwrXM9vE4SOZoy1u750WXCTZTcxbNszPNVyPCg9tv4YFrQWAOADR4EGgExRoAEgUxRoAMgUBRoAMkWBBoBMUaABIFMUaADIVL43qhSzCsX8Ye7YmgC6e1kmIN+6ALq5ySOYmxDamG64X8T0+FyzddcIXZLK0jQznzUbadKvQeYmBbkbZVwz9dbc6BLdOWKawbtm7Wb80/WZ5PJWveTyqjTLzet396kUnfT2g2nIH80XOkhSa274isUN6ee7Y+C+FcF8sYXa9M0+0dwQF8v0OVombvYJrf9SDYkraADIFgUaADJFgQaATFGgASBTFGgAyBQFGgAyRYEGgEzlm4NWUBjw+dGWJuNqMoaFaQZvc8QmwhhMw36Z/buMrOlDrs4CmvFHM0e1+VKD1rzGwoyh001fGzQuZxvNQYgmCBzTOe0Y0g3zg8noluYbDVpNJ5e7bvHRnQRlOuMrc/xqkzEubQZZquv0HJdtOqfcmi99KMyXIkTTkL82c1i561dzL0CbOIaz5rWfxxU0AGSKAg0AmaJAA0CmKNAAkCkKNABkigINAJmiQANAprLNQbehUhvmH15pgsimla1Ck/5ccgnPyuQzZSKqxYDX1d9+YzK4MT3+uIDPXTeHrt9xEUzONiwxy9M5YJtFL9OvsTZTYGLGCqbptutHHStzjF1Pb5dzNie5OcXVtulzrDI5afN0SVLHZMkb90ZxWXozBtdyuTBNs4M5xoXcAAYfo2IB83duHwCALFGgASBTFGgAyBQFGgAyRYEGgExRoAEgUxRoAMhUtjnoKrSqBvTcLYp0fvIS2zmr42bF5ZCjyXcG00vZRZRNP+kw63v1FuajueyYfsF1OiPamn630yaoW5ggdGH6MRetC0K7oLDbv8nSl+4guiy+yWHH9Pi6jblXwPTTjh2T8y7TPcslSSbvL5clr9PHeKZjziFzR0NHM+b5l3a/RFEmjlFq2ZwxAACyRIEGgExRoAEgUxRoAMgUBRoAMkWBBoBMUaABIFPZ5qA7oVBnQNa0NX1mY5vuRZzMJ0oKZrnrhdyajKyL4JpWxK5Nroquz0HHdATVZrFdv+bG9LR2PasLMwAzfHVMQ+ima3LGptevzXGbZsTR9Zs2B6gs0suni/Q52jG9mIto+k27k1BSXafn0BwC289ZJutemKy5O4d8P2nTkzwxRTPp2wR+aR8AgCxRoAEgUxRoAMgUBRoAMkWBBoBMUaABIFMUaADI1CXloCcmJvTEE0/o61//unbu3ClJijHqm9/8pp599lm98847WrNmjZ5++mn91m/91qK2PRu6qkJ3wKDTGc7WdGptzasOpo/trAkyF6bfs9p0RjW6Xsi2YXR6sSQFs1JhUqJFafpBmxyxy2G3zfzHvr/cBrXT54Dp2K2mSQdVY+vOQcOsYGLWCraftHn9Zn4Kk9MtXJhfUjBjmDb9nCtzjpYm616YrHZt7ldQla4D9n6GRL/raFPY7+1jQWvN49ChQ3r22Wd12223zXl8+/bt2rFjh3bt2qVDhw5pdHRUGzZs0KlTpy52VwBwXbqoAv2LX/xCX/7yl/W3f/u3+uhHP9p/PMaonTt3auvWrbr//vu1evVqPf/883r33Xe1e/fuyzZoALgeXFSBfvjhh/WHf/iH+v3f//05jx87dkxTU1MaGxvrP9br9bRu3TodPHhw3m1NT0/r5MmTc34AABfxf9B79uzRv/7rv+rQoUMXLJuampIkjYyMzHl8ZGREP/nJT+bd3sTEhL75zW8udhgAcM1b1BX05OSkvv71r+vv//7vdcMNNwxcL7yvEUyM8YLHztuyZYtOnDjR/5mcnFzMkADgmrWoK+hXX31Vx48f11133dV/rGkaHThwQLt27dJbb70l6dyV9C233NJf5/jx4xdcVZ/X6/XU6y3gG4IB4DqzqCvoL33pSzp69KiOHDnS/7n77rv15S9/WUeOHNGv//qva3R0VPv27es/Z2ZmRvv379fnPve5yz54ALiWLeoKemhoSKtXr57z2Ec+8hH9yq/8Sv/xzZs3a9u2bVq1apVWrVqlbdu2aenSpXrggQcWNbA6NqoHhGVdm9iiNBlRk0EdtN/zqpBeHoNJ2ZoMcTT5zWA/VxfwuWv6EUels9qF6ffsmum6Y9CYIG4o0jnpQf+ldl5pBlDEdAa3a3oRN6VLQqfnr1L69ak1GWFzDhWVOUfM+Gvz+iUpmKx1ZV6DYvocrM0Y22COQUwfg3bGjN9O4eDy2s6kn9vfx8JWW7jHH39cZ86c0aZNm/o3quzdu1dDQ0OXe1cAcE275AL98ssvz/l3CEHj4+MaHx+/1E0DwHWNXhwAkCkKNABkigINAJmiQANApijQAJCpyx6zu1xiExWb+XOSjYkZx8b0WjUfS64fs834mueHIj2+0vSpja5PrsmfSlI0We5yxvRzNh2Vo8kxx9rksGdMVtzltDumX7WZ46Jj3hqu37TJ0jeFOYkbE5Qt0q9fpelXbXLiLkfecQ29Jbl2y9E2BTfniLkhYtbcMVGat3kl0xT7EsYXTQb7PK6gASBTFGgAyBQFGgAyRYEGgExRoAEgUxRoAMgUBRoAMpVtDlqa1sDOz026V27jegGbz6XSBK1bE6A07ajVmGl37aRNAlbFAj5225geQ+Viur10hvXG6fQcTZt+yNOm129Tn04uL9qPJJfLfIlPYSaxlm1onVw8a955sTU9yV0O3PRCNhFguZRzcDlxSaXJKbvbFewgzAZclrsO6VlozQvomjoQEueI65V9HlfQAJApCjQAZIoCDQCZokADQKYo0ACQKQo0AGSKAg0Amco3B10HqZ7/86OtTK/ikM6IVm06g+ha9bp8ZWV6Lcs83zW6rU0v5rYx25fUdf2Cu+nP7m6bPgZNazKi5vkhpo9Rz+SMbU9v1wzYSo8vmAxtWaf3f9b0yw5luldxUabfA64Vc1maew1qn+Mt3Tlgrg9jcSa5fNbNkXkfVVW6/Lme4bMz6TloE2+xWd9OWxJX0ACQLQo0AGSKAg0AmaJAA0CmKNAAkCkKNABkigINAJnKNgfdKqoZkDV1KV+XY27MFkx8UoXJR8ZUAFJ+0t3z0wliqbQz5Hv1dlxT6Y7pSm36GYfZ9KuoTE5YhTmGHZPTNXPsjrHrl11GkxM/O5Pevpu/Nj2+1jQVD+b4hpgen21aLml2Np3VLsv0O6E0b+SyTs/Ru43ZfzRZf9fv2YyvjYOfP2s7cp/DFTQAZIoCDQCZokADQKYo0ACQKQo0AGSKAg0AmcouZhffi6acPXt28Epm1EWTjli5GFrhlpt4TWEiSJVpNViYF+gaPbrxS1JTpT+bO+7MMO0865np5PLpM+nlTeP6MZo5atPL3TGqTcyuNVPsYnZnTcxuxsTsCnNt5d4DLmZXDWj1+/934K/tQpOeg8LF7Mw52JxNv8YzZv+laZlau3alLs+bOEZnzpxrpRoTUTxJCtGt8QH7j//4D61YseJqDwMArrjJyUndeuutA5dnV6DbttXPfvYzDQ0NKYSgkydPasWKFZqcnNSyZcuu9vA+lJjDS8P8XTrmcK4Yo06dOqXly5erSPw2kt1/cRRFMe8nyrJlyziwl4g5vDTM36VjDv+/4eFhuw5/JASATFGgASBT2RfoXq+nb3zjG+r1eld7KB9azOGlYf4uHXN4cbL7IyEA4Jzsr6AB4HpFgQaATFGgASBTFGgAyBQFGgAylX2B/va3v62VK1fqhhtu0F133aV/+Zd/udpDytaBAwd0zz33aPny5Qoh6Ac/+MGc5TFGjY+Pa/ny5VqyZInWr1+vN9544+oMNjMTExP67Gc/q6GhId18882677779NZbb81Zh/lLe+aZZ3Tbbbf17xZcu3atfvSjH/WXM3+Ll3WB/t73vqfNmzdr69ateu211/S7v/u72rhxo376059e7aFl6fTp07r99tu1a9eueZdv375dO3bs0K5du3To0CGNjo5qw4YNOnXq1Ac80vzs379fDz/8sF555RXt27dPdV1rbGxMp0+f7q/D/KXdeuutevLJJ3X48GEdPnxYX/ziF3Xvvff2izDzdxFixn77t387PvTQQ3Me+43f+I34F3/xF1dpRB8ekuILL7zQ/3fbtnF0dDQ++eST/cfOnj0bh4eH41//9V9fhRHm7fjx41FS3L9/f4yR+btYH/3oR+Pf/d3fMX8XKdsr6JmZGb366qsaGxub8/jY2JgOHjx4lUb14XXs2DFNTU3Nmc9er6d169Yxn/M4ceKEJOmmm26SxPwtVtM02rNnj06fPq21a9cyfxcp2wL985//XE3TaGRkZM7jIyMjmpqaukqj+vA6P2fMpxdj1GOPPabPf/7zWr16tSTmb6GOHj2qG2+8Ub1eTw899JBeeOEFfeYzn2H+LlJ27UbfL4S532oQY7zgMSwc8+k98sgjev311/XjH//4gmXMX9qnP/1pHTlyRP/93/+tf/iHf9CDDz6o/fv395czf4uT7RX0xz72MZVlecGn6/Hjxy/4FIY3OjoqScyn8eijj+qHP/yhXnrppTl9yZm/hel2u/rkJz+pu+++WxMTE7r99tv1rW99i/m7SNkW6G63q7vuukv79u2b8/i+ffv0uc997iqN6sNr5cqVGh0dnTOfMzMz2r9/P/Opc1dyjzzyiL7//e/rn//5n7Vy5co5y5m/ixNj1PT0NPN3sa7iHyitPXv2xE6nE7/zne/EN998M27evDl+5CMfif/+7/9+tYeWpVOnTsXXXnstvvbaa1FS3LFjR3zttdfiT37ykxhjjE8++WQcHh6O3//+9+PRo0fjn/zJn8Rbbrklnjx58iqP/Or72te+FoeHh+PLL78c33777f7Pu+++21+H+UvbsmVLPHDgQDx27Fh8/fXX4xNPPBGLooh79+6NMTJ/FyPrAh1jjE8//XT8xCc+Ebvdbrzzzjv7sSdc6KWXXoqSLvh58MEHY4znomLf+MY34ujoaOz1evELX/hCPHr06NUddCbmmzdJ8bnnnuuvw/yl/dmf/Vn/vfqrv/qr8Utf+lK/OMfI/F0M+kEDQKay/T9oALjeUaABIFMUaADIFAUaADJFgQaATFGgASBTFGgAyBQFGgAyRYEGgExRoAEgUxRoAMjU/wNKP7XFv6mxUAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Export RGB\n",
    "np_arr_b4 = np.expand_dims(b4, 2)\n",
    "np_arr_b3 = np.expand_dims(b3, 2)\n",
    "np_arr_b2 = np.expand_dims(b2, 2)\n",
    "rgb_img = np.concatenate((np_arr_b4, np_arr_b3, np_arr_b2), 2)\n",
    "rgb_img = (255*((rgb_img)/3000)).astype('uint8')\n",
    "plt.imshow(rgb_img)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run\n",
    "This code runs the fuction over the date range.  Sentinel-2 flew over 2023-10-10, which appears to be the last day of the range; this is critical to determine the correlating *in situ* measurement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "start = pd.date_range(start= '2017-08-08' , end='2024-06-22' , \n",
    "              freq='5d')\n",
    "end = pd.date_range(start='2017-08-13' , end='2024-06-27' , \n",
    "              freq='5d')\n",
    "dates = pd.DataFrame ({'start': start ,  'end': end})\n",
    "#print(dates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define function\n",
    "def pulldata(startDate, endDate):\n",
    "    # Define source data\n",
    "    image = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')\\\n",
    "        .filterDate(startDate, endDate)\\\n",
    "        .select('B2', 'B3', 'B4', 'B8', 'B11')\n",
    "         #.filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 10)\\\n",
    "    \n",
    "    # CRS is not the same.\n",
    "    # proj = image.first().select('B2').projection() # EPSG:32656, UTM zone 56N (Siberia?)\n",
    "    # proj = balule.projection() # EPSG:4326\n",
    "\n",
    "    # Export arrays\n",
    "    # https://gist.github.com/jdbcode/f4d56d72f7fc5beeaa3859999b1f5c3d\n",
    "    # https://gist.github.com/jdbcode/f4d56d72f7fc5beeaa3859999b1f5c3d?permalink_comment_id=3355627#gistcomment-3355627\n",
    "    mosaic = image.median().reproject(crs='EPSG:32736', scale=10) # This allows us to set the resolution.\n",
    "    band_arrs = mosaic.sampleRectangle(region=aoi, defaultValue = 0)\n",
    "    \n",
    "    # Get individual band arrays.\n",
    "    band_arr_b2 = band_arrs.get('B2')   # Blue\n",
    "    band_arr_b3 = band_arrs.get('B3')   # Green\n",
    "    band_arr_b4 = band_arrs.get('B4')   # Red\n",
    "    band_arr_b5 = band_arrs.get('B5')   # Red edge\n",
    "    band_arr_b8 = band_arrs.get('B8')   # NIR\n",
    "    band_arr_b11 = band_arrs.get('B11') # SWIR\n",
    "\n",
    "    # Transfer the arrays from server to client and cast as np array.\n",
    "    b2 = np.array(band_arr_b2.getInfo())   # b2  Blue\n",
    "    b3 = np.array(band_arr_b3.getInfo())   # b3  Green\n",
    "    b4 = np.array(band_arr_b4.getInfo())   # b4  Red\n",
    "    b5 = np.array(band_arr_b5.getInfo())   # b4  Red edge\n",
    "    b8 = np.array(band_arr_b8.getInfo())   # b8  NIR\n",
    "    b11 = np.array(band_arr_b11.getInfo()) # b11 SWIR\n",
    "    if(np.max(b2)!=0):\n",
    "        # Normalized Difference Water Index (NDWI)\n",
    "        # NDWI = ( G - NIR ) / ( G + NIR )\n",
    "        ndwiG = (b3-b8)/(b3+b8) # Gao\n",
    "        # NDWI = ( NIR - SWIR ) / ( NIR + SWIR )\n",
    "        ndwiM = (b8-b11)/(b8+b11) # McFeeters\n",
    "        # NDWI = ( G - SWIR ) / ( G + SWIR )\n",
    "        mndwi = (b3-b11)/(b3+b11) # Modified NDWI\n",
    "        water = ndwiG > -0.02\n",
    "\n",
    "        TSS1 = np.NAN\n",
    "        Secchi1 = np.NAN\n",
    "        TSS2 = np.NAN\n",
    "        Secchi2 = np.NAN\n",
    "        TSS3 = np.NAN\n",
    "        Secchi3 = np.NAN\n",
    "        TSS4 = np.NAN\n",
    "        Ratio = np.NAN\n",
    "\n",
    "        TSS1 = (b3+b4)/2\n",
    "        TSS1 = TSS1 * water\n",
    "        TSS1 = np.sum(TSS1) / np.sum(TSS1>0)\n",
    "\n",
    "        Secchi1 = (b2/b4)\n",
    "        Secchi1 = Secchi1 * water\n",
    "        Secchi1 = np.sum(Secchi1)/ np.sum(Secchi1>0)\n",
    "\n",
    "        TSS2 = (b3/b4)\n",
    "        TSS2 = TSS2 * water\n",
    "        TSS2 = np.sum(TSS2) / np.sum(TSS2>0)\n",
    "\n",
    "        Secchi2 = (b4/b3)\n",
    "        Secchi2 = Secchi2 * water\n",
    "        Secchi2 = np.sum(Secchi2)/ np.sum(Secchi2>0)\n",
    "\n",
    "        TSS3 = (b8/b3 , b8/b4)\n",
    "        TSS3 = TSS3 * water\n",
    "        TSS3 = np.sum(TSS3) / np.sum(TSS3>0)\n",
    "\n",
    "        Secchi3 = (b4/b2)+b2\n",
    "        Secchi3 = Secchi3 * water\n",
    "        Secchi3 = np.sum(Secchi3)/ np.sum(Secchi3>0)\n",
    "\n",
    "        TSS4 = (b4/b3)+b8\n",
    "        TSS4 = TSS4 * water\n",
    "        TSS4 = np.sum(TSS4) / np.sum(TSS4>0)\n",
    "\n",
    "        Ratio = (ndwiG/ndwiM)\n",
    "        Ratio = Ratio * water\n",
    "        Ratio = np.sum(Ratio)/ np.sum(Ratio>0) \n",
    "\n",
    "\n",
    "        f = open(\"sentinel.csv\", \"a\")\n",
    "        f.write(str(site) + \", \" + str(startDate) + \", \" + str(endDate) + \", \" + \n",
    "                str(TSS1) + \", \" + str(Secchi1) + \", \" + \n",
    "                str(TSS2) + \", \" + str(Secchi2) + \", \" + \n",
    "                str(TSS3) + \", \" + str(Secchi3) + \", \" + \n",
    "                str(TSS4) + \",\" + str(Ratio) + '\\n') \n",
    "            \n",
    "        f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/3m/3_cnzbmx5r72vkz_1lrm0pyh0000gn/T/ipykernel_11409/1592690793.py:55: RuntimeWarning: invalid value encountered in scalar divide\n",
      "  TSS1 = np.sum(TSS1) / np.sum(TSS1>0)\n",
      "/var/folders/3m/3_cnzbmx5r72vkz_1lrm0pyh0000gn/T/ipykernel_11409/1592690793.py:59: RuntimeWarning: invalid value encountered in scalar divide\n",
      "  Secchi1 = np.sum(Secchi1)/ np.sum(Secchi1>0)\n",
      "/var/folders/3m/3_cnzbmx5r72vkz_1lrm0pyh0000gn/T/ipykernel_11409/1592690793.py:63: RuntimeWarning: invalid value encountered in scalar divide\n",
      "  TSS2 = np.sum(TSS2) / np.sum(TSS2>0)\n",
      "/var/folders/3m/3_cnzbmx5r72vkz_1lrm0pyh0000gn/T/ipykernel_11409/1592690793.py:67: RuntimeWarning: invalid value encountered in scalar divide\n",
      "  Secchi2 = np.sum(Secchi2)/ np.sum(Secchi2>0)\n",
      "/var/folders/3m/3_cnzbmx5r72vkz_1lrm0pyh0000gn/T/ipykernel_11409/1592690793.py:71: RuntimeWarning: invalid value encountered in scalar divide\n",
      "  TSS3 = np.sum(TSS3) / np.sum(TSS3>0)\n",
      "/var/folders/3m/3_cnzbmx5r72vkz_1lrm0pyh0000gn/T/ipykernel_11409/1592690793.py:75: RuntimeWarning: invalid value encountered in scalar divide\n",
      "  Secchi3 = np.sum(Secchi3)/ np.sum(Secchi3>0)\n",
      "/var/folders/3m/3_cnzbmx5r72vkz_1lrm0pyh0000gn/T/ipykernel_11409/1592690793.py:79: RuntimeWarning: invalid value encountered in scalar divide\n",
      "  TSS4 = np.sum(TSS4) / np.sum(TSS4>0)\n",
      "/var/folders/3m/3_cnzbmx5r72vkz_1lrm0pyh0000gn/T/ipykernel_11409/1592690793.py:81: RuntimeWarning: divide by zero encountered in divide\n",
      "  Ratio = (ndwiG/ndwiM)\n",
      "/var/folders/3m/3_cnzbmx5r72vkz_1lrm0pyh0000gn/T/ipykernel_11409/1592690793.py:82: RuntimeWarning: invalid value encountered in multiply\n",
      "  Ratio = Ratio * water\n",
      "/var/folders/3m/3_cnzbmx5r72vkz_1lrm0pyh0000gn/T/ipykernel_11409/1592690793.py:83: RuntimeWarning: invalid value encountered in scalar divide\n",
      "  Ratio = np.sum(Ratio)/ np.sum(Ratio>0)\n",
      "/var/folders/3m/3_cnzbmx5r72vkz_1lrm0pyh0000gn/T/ipykernel_11409/1592690793.py:83: RuntimeWarning: divide by zero encountered in scalar divide\n",
      "  Ratio = np.sum(Ratio)/ np.sum(Ratio>0)\n",
      "/Users/davidkahler/opt/miniconda3/lib/python3.9/site-packages/numpy/core/fromnumeric.py:88: RuntimeWarning: invalid value encountered in reduce\n",
      "  return ufunc.reduce(obj, axis, dtype, out, **passkwargs)\n"
     ]
    }
   ],
   "source": [
    "# Run over all dates\n",
    "for i in range(len(dates)):\n",
    "    pulldata(dates[\"start\"][i], dates[\"end\"][i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/3m/3_cnzbmx5r72vkz_1lrm0pyh0000gn/T/ipykernel_11409/1592690793.py:83: RuntimeWarning: divide by zero encountered in scalar divide\n",
      "  Ratio = np.sum(Ratio)/ np.sum(Ratio>0)\n"
     ]
    }
   ],
   "source": [
    "pulldata(startDate, endDate)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}