HLS_Shortgrass_Explorer_v11.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import dask, concurrent.futures, time, warnings, os, re, pickle\n",
    "from osgeo import gdal\n",
    "import requests as r\n",
    "import panel as pn\n",
    "pn.extension('echarts')\n",
    "import param as pm\n",
    "import pandas as pd\n",
    "from collections import OrderedDict as odict\n",
    "import numpy as np\n",
    "from dask.distributed import LocalCluster, Client\n",
    "import xarray as xr\n",
    "import hvplot.pandas\n",
    "import hvplot.xarray\n",
    "#import warnings\n",
    "#warnings.filterwarnings('ignore')\n",
    "from urllib.request import urlopen\n",
    "from xml.etree.ElementTree import parse,fromstring\n",
    "from affine import Affine\n",
    "from pandas import to_datetime\n",
    "import jinja2 as jj2\n",
    "from rasterio.crs import CRS\n",
    "from tempfile import NamedTemporaryFile\n",
    "from datetime import datetime, timedelta\n",
    "from netrc import netrc\n",
    "from subprocess import Popen\n",
    "from pyproj import Proj\n",
    "from src.hls_funcs.masks import mask_hls\n",
    "from src.hls_funcs.predict import pred_cov, pred_bm, pred_bm_se, xr_cdf\n",
    "from src.hls_funcs.indices import ndvi_func\n",
    "import cartopy.crs as ccrs\n",
    "from bokeh.models.formatters import PrintfTickFormatter\n",
    "from bokeh.models import NumeralTickFormatter\n",
    "import stackstac\n",
    "from subprocess import Popen, DEVNULL, STDOUT\n",
    "from getpass import getpass\n",
    "from sys import platform\n",
    "import itertools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a LUT dict including the HLS product bands mapped to names\n",
    "lut = {'HLSS30':\n",
    "       {'B01': 'COASTAL-AEROSOL',\n",
    "        'B02': 'BLUE', \n",
    "        'B03': 'GREEN', \n",
    "        'B04': 'RED', \n",
    "        'B05': 'RED-EDGE1',\n",
    "        'B06': 'RED-EDGE2', \n",
    "        'B07': 'RED-EDGE3',\n",
    "        'B08': 'NIR-Broad',\n",
    "        'B8A': 'NIR1', \n",
    "        'B09': 'WATER-VAPOR',\n",
    "        'B10': 'CIRRUS',\n",
    "        'B11': 'SWIR1', \n",
    "        'B12': 'SWIR2', \n",
    "        'Fmask': 'FMASK'},\n",
    "       'HLSL30': \n",
    "       {'B01': 'COASTAL-AEROSOL',\n",
    "        'B02': 'BLUE', \n",
    "        'B03': 'GREEN', \n",
    "        'B04': 'RED', \n",
    "        'B05': 'NIR1',\n",
    "        'B06': 'SWIR1',\n",
    "        'B07': 'SWIR2', \n",
    "        'B09': 'CIRRUS', \n",
    "        'B10': 'TIR1', \n",
    "        'B11': 'TIR2', \n",
    "        'Fmask': 'FMASK'}}\n",
    "\n",
    "# List of all available/acceptable band names\n",
    "all_bands = ['ALL', 'COASTAL-AEROSOL', 'BLUE', 'GREEN', 'RED', 'RED-EDGE1', 'RED-EDGE2', 'RED-EDGE3', \n",
    "             'NIR1', 'SWIR1', 'SWIR2', 'CIRRUS', 'TIR1', 'TIR2', 'WATER-VAPOR', 'FMASK']\n",
    "\n",
    "needed_bands = ['BLUE', 'GREEN', 'RED', 'NIR1', 'SWIR1', 'SWIR2', 'FMASK']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def NASA_CMR_STAC(hls_data, aws):\n",
    "    stac = 'https://cmr.earthdata.nasa.gov/stac/' # CMR-STAC API Endpoint\n",
    "    stac_response = r.get(stac).json()            # Call the STAC API endpoint\n",
    "    stac_lp = [s for s in stac_response['links'] if 'LP' in s['title']]  # Search for only LP-specific catalogs\n",
    "\n",
    "    # LPCLOUD is the STAC catalog we will be using and exploring today\n",
    "    lp_cloud = r.get([s for s in stac_lp if s['title'] == 'LPCLOUD'][0]['href']).json()\n",
    "    lp_links = lp_cloud['links']\n",
    "    lp_search = [l['href'] for l in lp_links if l['rel'] == 'search'][0]  # Define the search endpoint\n",
    "    lim = 100\n",
    "    search_query = f\"{lp_search}?&limit={lim}\"    # Add in a limit parameter to retrieve 100 items at a time.\n",
    "    bbox_num=[-104.79107047,   40.78311181, -104.67687336,   40.87008987]\n",
    "    bbox = f'{bbox_num[0]},{bbox_num[1]},{bbox_num[2]},{bbox_num[3]}'  # Defined from ROI bounds\n",
    "    search_query2 = f\"{search_query}&bbox={bbox}\"                                                  # Add bbox to query\n",
    "    date_time = hls_data['date_range'][0]+'/'+hls_data['date_range'][1]  # Define start time period / end time period\n",
    "    search_query3 = f\"{search_query2}&datetime={date_time}\"  # Add to query that already includes bbox\n",
    "    collections = r.get(search_query3).json()['features']    \n",
    "    hls_collections = [c for c in collections if 'HLS' in c['collection']]\n",
    "    s30_items = [h for h in hls_collections if h['collection'] == 'HLSS30.v1.5']  # Grab HLSS30 collection\n",
    "    l30_items = [h for h in hls_collections if h['collection'] == 'HLSL30.v1.5']  # Grab HLSL30 collection\n",
    "    \n",
    "    if aws:\n",
    "        for stac in s30_items:\n",
    "            for band in stac['assets']:\n",
    "                stac['assets'][band]['href'] = stac['assets'][band]['href'].replace('https://lpdaac.earthdata.nasa.gov/lp-prod-protected', \n",
    "                                                                                    '/vsis3/lp-prod-protected')\n",
    "                stac['assets'][band]['href'] = stac['assets'][band]['href'].replace('https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected', \n",
    "                                                                                    '/vsis3/lp-prod-protected')\n",
    "                \n",
    "        for stac in l30_items:\n",
    "            for band in stac['assets']:\n",
    "                stac['assets'][band]['href'] = stac['assets'][band]['href'].replace('https://lpdaac.earthdata.nasa.gov/lp-prod-protected', \n",
    "                                                                                    '/vsis3/lp-prod-protected')\n",
    "                stac['assets'][band]['href'] = stac['assets'][band]['href'].replace('https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected', \n",
    "                                                                                    '/vsis3/lp-prod-protected')\n",
    "    return {'S30': s30_items,\n",
    "            'L30': l30_items}\n",
    "\n",
    "def setup_netrc(creds,aws):\n",
    "    urs = 'urs.earthdata.nasa.gov' \n",
    "    try:\n",
    "        netrcDir = os.path.expanduser(\"~/.netrc\")\n",
    "        netrc(netrcDir).authenticators(urs)[0]\n",
    "        del netrcDir\n",
    "\n",
    "    # Below, create a netrc file and prompt user for NASA Earthdata Login Username and Password\n",
    "    except FileNotFoundError:\n",
    "        homeDir = os.path.expanduser(\"~\")\n",
    "        Popen('touch {0}.netrc | chmod og-rw {0}.netrc | echo machine {1} >> {0}.netrc'.format(homeDir + os.sep, urs), shell=True)\n",
    "        Popen('echo login {} >> {}.netrc'.format(creds[0], homeDir + os.sep), shell=True)\n",
    "        Popen('echo password {} >> {}.netrc'.format(creds[1], homeDir + os.sep), shell=True)\n",
    "        del homeDir\n",
    "\n",
    "    # Determine OS and edit netrc file if it exists but is not set up for NASA Earthdata Login\n",
    "    except TypeError:\n",
    "        homeDir = os.path.expanduser(\"~\")\n",
    "        Popen('echo machine {1} >> {0}.netrc'.format(homeDir + os.sep, urs), shell=True)\n",
    "        Popen('echo login {} >> {}.netrc'.format(creds[0], homeDir + os.sep), shell=True)\n",
    "        Popen('echo password {} >> {}.netrc'.format(creds[1], homeDir + os.sep), shell=True)\n",
    "        del homeDir\n",
    "    del urs\n",
    "    if aws:\n",
    "        return(r.get('https://lpdaac.earthdata.nasa.gov/s3credentials').json())\n",
    "    else:\n",
    "        return('')\n",
    "\n",
    "def build_xr(stac_dict):\n",
    "    try:\n",
    "        s30_stack = stackstac.stack(stac_dict['S30'], epsg=32613, resolution=30, assets=[i for i in lut['HLSS30'] if lut['HLSS30'][i] in needed_bands],\n",
    "                                   chunksize=(4000, 4000)) #, reader=stackstac.reader_protocol.FakeReader\n",
    "        s30_stack['band'] = [lut['HLSS30'][b] for b in s30_stack['band'].values]\n",
    "        s30_stack['time'] = [datetime.fromtimestamp(t) for t in s30_stack.time.astype('int').values//1000000000]\n",
    "        s30_stack = s30_stack.to_dataset(dim='band').reset_coords(['end_datetime', 'start_datetime'], drop=True)\n",
    "    except ValueError:\n",
    "        s30_stack = None\n",
    "    try:\n",
    "        l30_stack = stackstac.stack(stac_dict['L30'], epsg=32613, resolution=30, assets=[i for i in lut['HLSL30'] if lut['HLSL30'][i] in needed_bands],\n",
    "                                   chunksize=(4000, 4000))\n",
    "        l30_stack['band'] = [lut['HLSL30'][b] for b in l30_stack['band'].values]\n",
    "        l30_stack['time'] = [datetime.fromtimestamp(t) for t in l30_stack.time.astype('int').values//1000000000]\n",
    "        l30_stack = l30_stack.to_dataset(dim='band').reset_coords(['name', 'end_datetime', 'start_datetime'], drop=True)\n",
    "    except ValueError:\n",
    "        l30_stack = None\n",
    "    if s30_stack is not None and l30_stack is not None:\n",
    "        hls_stack = xr.concat([s30_stack, l30_stack], dim='time')\n",
    "    elif s30_stack is not None:\n",
    "        hls_stack = s30_stack\n",
    "    elif l30_stack is not None:\n",
    "        hls_stack = l30_stack\n",
    "    else:\n",
    "        print('No data found for date range')\n",
    "    return hls_stack.chunk({'time': 1, 'y': -1, 'x': -1})\n",
    "    \n",
    "def get_hls(creds, hls_data={}, aws=False):\n",
    "    #Seteup creds\n",
    "    \n",
    "    s3_cred = setup_netrc(creds,aws=aws)\n",
    "    #define gdalenv\n",
    "    if aws:\n",
    "        \n",
    "        env = dict(GDAL_DISABLE_READDIR_ON_OPEN='FALSE', \n",
    "                   #AWS_NO_SIGN_REQUEST='YES',\n",
    "                   GDAL_MAX_RAW_BLOCK_CACHE_SIZE='200000000',\n",
    "                   GDAL_SWATH_SIZE='200000000',\n",
    "                   VSI_CURL_CACHE_SIZE='200000000',\n",
    "                   CPL_VSIL_CURL_ALLOWED_EXTENSIONS='TIF',\n",
    "                   GDAL_HTTP_UNSAFESSL='YES',\n",
    "                   GDAL_HTTP_COOKIEFILE=os.path.expanduser('~/cookies.txt'),\n",
    "                   GDAL_HTTP_COOKIEJAR=os.path.expanduser('~/cookies.txt'),\n",
    "                   AWS_REGION='us-west-2',\n",
    "                   AWS_SECRET_ACCESS_KEY=s3_cred['secretAccessKey'],\n",
    "                   AWS_ACCESS_KEY_ID=s3_cred['accessKeyId'],\n",
    "                   AWS_SESSION_TOKEN=s3_cred['sessionToken'])\n",
    "    else:\n",
    "        env = dict(GDAL_DISABLE_READDIR_ON_OPEN='EMPTY_DIR', \n",
    "                   AWS_NO_SIGN_REQUEST='YES',\n",
    "                   GDAL_MAX_RAW_BLOCK_CACHE_SIZE='200000000',\n",
    "                   GDAL_SWATH_SIZE='200000000',\n",
    "                   VSI_CURL_CACHE_SIZE='200000000',\n",
    "                   GDAL_HTTP_COOKIEFILE=os.path.expanduser('~/cookies.txt'),\n",
    "                   GDAL_HTTP_COOKIEJAR=os.path.expanduser('~/cookies.txt'))\n",
    "\n",
    "\n",
    "    os.environ.update(env)\n",
    "    \n",
    "    catalog = NASA_CMR_STAC(hls_data, aws)\n",
    "    da  = build_xr(catalog)\n",
    "    return da"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_dict = {'date_range': [str(datetime(2021, 1, 1).date()), str(datetime.now().date())]}\n",
    "tmp_data = get_hls(['spkearney', '1mrChamu'], hls_data=data_dict, aws=True)\n",
    "da = tmp_data.loc[dict(x=slice(517587.0, 527283.0), y=slice(4524402.0, 4514699.0))]\n",
    "da_mask = mask_hls(da['FMASK'])\n",
    "da = da.where(da_mask == 0)\n",
    "t0 = time.time()\n",
    "print(da.isel(time=1).mean().values)\n",
    "t1 = time.time()\n",
    "print(t1-t0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster = LocalCluster(threads_per_worker=1)\n",
    "cl = Client(cluster)\n",
    "cl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_dict = {'date_range': [str(datetime(2021, 1, 1).date()), str(datetime.now().date())]}\n",
    "tmp_data = get_hls(['spkearney', '1mrChamu'], hls_data=data_dict, aws=True)\n",
    "da = tmp_data.loc[dict(x=slice(517587.0, 527283.0), y=slice(4524402.0, 4514699.0))]\n",
    "da_mask = mask_hls(da['FMASK'])\n",
    "da = da.where(da_mask == 0)\n",
    "t0 = time.time()\n",
    "print(da.isel(time=1).mean().values)\n",
    "t1 = time.time()\n",
    "print(t1-t0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cartopy import crs\n",
    "import geoviews as gv\n",
    "import holoviews as hv\n",
    "from copy import deepcopy\n",
    "from holoviews import streams\n",
    "import affine\n",
    "from skimage.draw import polygon\n",
    "\n",
    "cov_dict = {'R': 'Dry veg',\n",
    "            'G': 'Green veg',\n",
    "            'B': 'Bare ground'}\n",
    "\n",
    "dat = da\n",
    "dat['time'] = dat.time.dt.floor(\"D\")\n",
    "dat = dat.rename(dict(time='date'))\n",
    "bm_mod = pickle.load(open('src/models/CPER_HLS_to_VOR_biomass_model_lr_simp.pk', 'rb'))\n",
    "cov_mod = pickle.load(open('src/models/CPER_HLS_to_LPI_cover_pls_binned_model.pk', 'rb'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class getData(pm.Parameterized):\n",
    "    bm_mod = pickle.load(open('src/models/CPER_HLS_to_VOR_biomass_model_lr_simp.pk', 'rb'))\n",
    "    cov_mod = pickle.load(open('src/models/CPER_HLS_to_LPI_cover_pls_binned_model.pk', 'rb'))\n",
    "    action = pm.Action(lambda x: x.param.trigger('action'), label='Load Data and Run Analysis')\n",
    "    get_stats = pm.Action(lambda y: y.param.trigger('get_stats'), label='Calculate static stats')\n",
    "    get_stats_ts = pm.Action(lambda y: y.param.trigger('get_stats_ts'), label='Calculate time-series stats')\n",
    "    username_input = pn.widgets.PasswordInput(name='NASA Earthdata Login', placeholder='Enter Username...')\n",
    "    password_input = pn.widgets.PasswordInput(name='', placeholder='Enter Password...')\n",
    "    year_picker = pn.widgets.Select(name='Year', options=[2021], value=2021)\n",
    "    thresh_picker = pn.widgets.IntSlider(name='Threshold', start=200, end=2000, step=25, value=500,\n",
    "                                 format=PrintfTickFormatter(format='%d kg/ha'))\n",
    "    \n",
    "    data_dict = {'date_range': [str(datetime(year_picker.value, 1, 1).date()),\n",
    "                                str(datetime(year_picker.value, 12, 31).date())]}\n",
    " \n",
    "    chunks = {'date': 1, 'y': 250, 'x': 250}\n",
    "    datCRS = crs.UTM(13)\n",
    "    mapCRS = crs.GOOGLE_MERCATOR\n",
    "    datProj = Proj(datCRS.proj4_init)\n",
    "    mapProj = Proj(mapCRS.proj4_init)\n",
    "    map_args = dict(crs=datCRS, rasterize=False, project=True, dynamic=True)\n",
    "    base_opts = dict(projection=mapCRS, backend='bokeh', xticks=None, yticks=None, width=900, height=700,\n",
    "                         padding=0, active_tools=['pan', 'wheel_zoom'], toolbar='left')\n",
    "    map_opts = dict(projection=mapCRS, responsive=False, xticks=None, yticks=None, width=900, height=700,\n",
    "                     padding=0, tools=['pan', 'wheel_zoom', 'box_zoom'],\n",
    "                     active_tools=['pan', 'wheel_zoom'], toolbar='left')\n",
    "    \n",
    "    poly_opts = dict(fill_color=['', ''], fill_alpha=[0.0, 0.0], line_color=['#1b9e77', '#d95f02'],\n",
    "                 line_width=[3, 3])    \n",
    "    bg_col='#ffffff'\n",
    "\n",
    "    css = '''\n",
    "    .bk.box1 {\n",
    "      background: #ffffff;\n",
    "      border-radius: 5px;\n",
    "      border: 1px black solid;\n",
    "    }\n",
    "    '''\n",
    "    \n",
    "    date = pn.widgets.DatePicker(name='Calendar', width=200)\n",
    "      \n",
    "    box_ctr_dat = [522435.0, 4519550.5]\n",
    "    box_ctr_coord = datProj(box_ctr_dat[0], box_ctr_dat[1], inverse=True)\n",
    "    box_ctr_map = mapProj(box_ctr_coord[0], box_ctr_coord[1], inverse=False)\n",
    "    \n",
    "    tiles = gv.tile_sources.EsriImagery.opts(**base_opts, level='glyph',\n",
    "                                                                   xlim=(box_ctr_map[0] - 5000,\n",
    "                                                                         box_ctr_map[0] + 5000),\n",
    "                                                                   ylim=(box_ctr_map[1] - 5000,\n",
    "                                                                         box_ctr_map[1] + 5000))\n",
    "    labels = gv.tile_sources.EsriReference.opts(**base_opts, level='overlay')\n",
    "       \n",
    "    basemap = gv.tile_sources.EsriImagery.opts(projection=mapCRS, backend='bokeh', level='glyph')\n",
    "    \n",
    "    base_rng = hv.streams.RangeXY(source=tiles)\n",
    "    \n",
    "    polys = hv.Polygons([])\n",
    "\n",
    "    poly_stream = streams.PolyDraw(source=polys, drag=True, num_objects=2,\n",
    "                                    show_vertices=True, styles=poly_opts)\n",
    "    edit_stream = streams.PolyEdit(source=polys, shared=True)\n",
    "    \n",
    "    gauge_obj = pn.indicators.Gauge(\n",
    "        name='Biomass', bounds=(0, 2500), format='{value} kg/ha',\n",
    "        colors=[(0.20, '#FF6E76'), (0.40, '#FDDD60'), (0.60, '#7CFFB2'), (1, '#58D9F9')],\n",
    "        num_splits=5, value=0, align='center', title_size=12,\n",
    "        start_angle=180, end_angle=0, height=200, width=300,\n",
    "    )\n",
    "    \n",
    "    def __init__(self, **params):\n",
    "        super(getData, self).__init__(**params)\n",
    "        self.da = None\n",
    "        self.da_sel = None\n",
    "        self.da_cov = None\n",
    "        self.da_bm = None\n",
    "        self.da_bm_se = None\n",
    "        self.da_thresh = None\n",
    "        self.all_maps = None\n",
    "        self.cov_map = None\n",
    "        self.bm_map = None\n",
    "        self.thresh_map = None\n",
    "        \n",
    "        self.cov_stats = None\n",
    "        self.bm_stats = None\n",
    "        self.thresh_stats = None\n",
    "        self.stats_titles = None\n",
    "                \n",
    "        self.view = self._create_view()\n",
    "        \n",
    "    \n",
    "    @pm.depends('action')\n",
    "    def proc_data(self):\n",
    "        message = 'Not yet launched'\n",
    "        if self.username_input.value != '':\n",
    "            message = 'username found'\n",
    "            try:\n",
    "                tmp_data = get_hls([self.username_input.value,self.password_input.value], \n",
    "                                   hls_data=self.data_dict, aws=True)\n",
    "                message = self.base_rng.x_range\n",
    "                coord_rng = self.mapProj(self.base_rng.x_range, self.base_rng.y_range, inverse=True)\n",
    "                message = coord_rng\n",
    "                dat_rng = np.round(self.datProj(coord_rng[0], coord_rng[1], inverse=False), 0)\n",
    "                message = dat_rng\n",
    "                self.da = tmp_data.loc[dict(x=slice(*dat_rng[0]), y=slice(*dat_rng[1][::-1]))]\n",
    "                self.da['time'] = self.da.time.dt.floor(\"D\")\n",
    "                self.da = self.da.rename(dict(time='date'))\n",
    "                self.da = self.da.chunk(self.chunks)\n",
    "                da_mask = mask_hls(self.da['FMASK'])\n",
    "                self.da = self.da.where(da_mask == 0)\n",
    "                self.da = self.da.sel(date=self.da['eo:cloud_cover'] < 80)\n",
    "                self.date.enabled_dates = [pd.Timestamp(x).to_pydatetime().date() for x in self.da['date'].values]\n",
    "                if self.date.value is None:\n",
    "                    self.da_sel = self.da.isel(date=-1).compute()\n",
    "                    self.date.value = pd.to_datetime(self.da_sel.date.values).date()\n",
    "                else:\n",
    "                    date_idx = self.date.enabled_dates.index(self.date.value)\n",
    "                    self.da_sel = self.da.isel(date=date_idx).compute()\n",
    "                message = 'Success!' + str(datetime.now())\n",
    "                return message\n",
    "            except:\n",
    "                return message + ': App Failure'\n",
    "        else:\n",
    "            return message\n",
    "      \n",
    "    @pm.depends('date.param')\n",
    "    def create_maps(self):\n",
    "        if self.da is not None and self.da_sel is not None:\n",
    "            if pd.to_datetime(self.da_sel.date.values).date() != self.date.value:\n",
    "                if self.date.value is not None:\n",
    "                    date_idx = self.date.enabled_dates.index(self.date.value)\n",
    "                    self.da_sel = self.da.isel(date=date_idx).compute()\n",
    "            if self.edit_stream.data is not None:\n",
    "                self.polys = self.edit_stream.element\n",
    "                self.poly_stream = streams.PolyDraw(source=self.polys, drag=True, num_objects=2,\n",
    "                                                    show_vertices=True, styles=self.poly_opts)\n",
    "                self.edit_stream = streams.PolyEdit(source=self.polys, shared=True)\n",
    "            elif self.poly_stream.data is not None:\n",
    "                self.polys = self.poly_stream.element\n",
    "                self.poly_stream = streams.PolyDraw(source=self.polys, drag=True, num_objects=2,\n",
    "                               show_vertices=True, styles=self.poly_opts)\n",
    "                self.edit_stream = streams.PolyEdit(source=self.polys, shared=True)\n",
    "            else:\n",
    "                self.polys = self.polys\n",
    "            da_bm = self.da_sel.map_blocks(pred_bm, template=self.da_sel['BLUE'],\n",
    "                                     kwargs=dict(model=self.bm_mod))\n",
    "            da_bm = da_bm.where(da_bm > 0)\n",
    "            da_bm.name = 'Biomass'\n",
    "            \n",
    "            da_cov_temp = self.da_sel[['BLUE', 'GREEN', 'RED', 'NIR1']].rename(dict(BLUE='SD', RED='BARE', NIR1='LITT'))\n",
    "            da_cov = self.da_sel.map_blocks(pred_cov, template=da_cov_temp,\n",
    "                         kwargs=dict(model=self.cov_mod))\n",
    "            da_cov = da_cov[['SD', 'GREEN', 'BARE']].to_array(dim='type')\n",
    "            da_cov = da_cov.where((da_cov < 1.0) | (da_cov.isnull()), 1.0)\n",
    "            da_cov.name = 'Cover'\n",
    "            \n",
    "            da_bm_se = self.da_sel.map_blocks(pred_bm_se, template=self.da_sel['BLUE'],\n",
    "                                         kwargs=dict(model=self.bm_mod))\n",
    "           \n",
    "            da_thresh_pre = (np.log(self.thresh_picker.value) - xr.ufuncs.log(da_bm)) / da_bm_se\n",
    "            da_thresh = da_thresh_pre.map_blocks(xr_cdf, template=self.da_sel['BLUE'])\n",
    "            #da_thresh = pred_bm_thresh(da_bm, da_bm_se, self.thresh_picker.value)\n",
    "            da_thresh.name = 'Threshold'\n",
    "            \n",
    "            da_ndvi = ndvi_func(self.da_sel)\n",
    "            \n",
    "            self.cov_map = da_cov.hvplot.rgb(x='x', y='y', bands='type',\n",
    "                                        **self.map_args).opts(**self.map_opts)\n",
    "            self.bm_map = da_bm.hvplot(x='x', y='y',\n",
    "                                  cmap='Inferno', clim=(100, 1000), colorbar=False,\n",
    "                                  **self.map_args).opts(**self.map_opts)\n",
    "            self.thresh_map = da_thresh.hvplot(x='x', y='y',\n",
    "                                          cmap='YlOrRd', clim=(0.05, 0.95), colorbar=False,\n",
    "                                          **self.map_args).opts(**self.map_opts)\n",
    "            \n",
    "            self.ndvi_map = da_ndvi.hvplot(x='x', y='y',\n",
    "                              cmap='Viridis', clim=(0.05, 0.50), colorbar=False,\n",
    "                              **self.map_args).opts(**self.map_opts)\n",
    "\n",
    "            self.da_cov = da_cov\n",
    "            self.da_bm = da_bm\n",
    "            self.da_bm_se = da_bm_se\n",
    "            self.da_thresh = da_thresh\n",
    "            self.da_ndvi = da_ndvi\n",
    "            \n",
    "            self.all_maps = pn.Tabs(\n",
    "                ('Cover', self.basemap * self.cov_map * self.polys),\n",
    "                ('Biomass', self.basemap * self.bm_map * self.polys),\n",
    "                ('Biomass threshold', self.basemap * self.thresh_map * self.polys), \n",
    "                ('Greenness (NDVI)', self.basemap * self.ndvi_map * self.polys))\n",
    "            \n",
    "            #self.poly_stream = streams.PolyDraw(source=self.polys, drag=True, num_objects=2,\n",
    "            #                        show_vertices=True, styles=self.poly_opts)\n",
    "        \n",
    "            return self.all_maps\n",
    "        else:\n",
    "            return pn.Column(self.tiles * self.labels)\n",
    "        \n",
    "    @pm.depends('get_stats')\n",
    "    def show_stats(self):\n",
    "        if self.poly_stream.data is not None:\n",
    "            markdown_list = []\n",
    "            thresh_list = []\n",
    "            bm_list = []\n",
    "            cov_list = []\n",
    "            stats_rows = []\n",
    "            for idx, ps_c in enumerate(self.poly_stream.data['line_color']):\n",
    "                r, c = polygon(self.poly_stream.data['xs'][idx]*-1, self.poly_stream.data['ys'][idx])\n",
    "                r = r * -1.0\n",
    "                \n",
    "                da_cov_tmp = self.cov_map[r][:,c].data\n",
    "                cov_factors = [k for k in app.cov_map.data.keys() if k not in ['y', 'x']]\n",
    "                cov_labels = [cov_dict[i] for i in cov_factors]\n",
    "                cov_labels.append('Litter')\n",
    "                cov_vals = [round(float(da_cov_tmp[f].mean()), 2) for f in cov_factors]\n",
    "                cov_vals.append(round(1.0 - np.sum(cov_vals), 2))\n",
    "                pct_fmt = NumeralTickFormatter(format=\"0%\")\n",
    "                cov_colors = hv.Cycle(['red', 'green', 'blue', 'orange'])\n",
    "                cov_scatter_tmp = hv.Overlay([hv.Scatter(f) for f in list(zip(cov_labels, cov_vals))]) \\\n",
    "                    .options({'Scatter': dict(xformatter=pct_fmt,\n",
    "                                              size=15,\n",
    "                                              fill_color=cov_colors,\n",
    "                                              line_color=cov_colors,\n",
    "                                              ylim=(0, 1))})\n",
    "                cov_spike_tmp = hv.Overlay([hv.Spikes(f) for f in cov_scatter_tmp])\\\n",
    "                    .options({'Spikes': dict(color=cov_colors, line_width=4,\n",
    "                                             labelled=[], invert_axes=True, color_index=None,\n",
    "                                             ylim=(0, 1))})\n",
    "                cov_list.append(pn.Column((cov_spike_tmp * cov_scatter_tmp).options(height=200,\n",
    "                                                                                    width=300,\n",
    "                                                                                    bgcolor=self.bg_col,\n",
    "                                                                                    toolbar=None),\n",
    "                                          css_classes=['box1'], margin=5, align='center'))\n",
    "                bm_dat_tmp = self.bm_map[r][:,c].data\n",
    "                bm_hist_tmp = bm_dat_tmp.hvplot.hist('Biomass', xlim=(0, 2000),\n",
    "                                                    bins=np.arange(0, 10000, 20))\\\n",
    "                    .opts(height=200, width=300, fill_color=ps_c, fill_alpha=0.6,\n",
    "                          line_color='black', line_width=0.5, line_alpha=0.6,\n",
    "                          bgcolor=self.bg_col, align='center').options(toolbar=None)\n",
    "                bm_gauge_obj = deepcopy(self.gauge_obj)\n",
    "                bm_gauge_obj.value = int(bm_dat_tmp.mean()['Biomass'])\n",
    "                bm_list.append(pn.Column(bm_gauge_obj,\n",
    "                                         css_classes=['box1'], margin=5, align='center'))\n",
    "                markdown = pn.pane.Markdown('# Region ' + str(idx+1) + ' :', align='center',\n",
    "                                            style={'font-family': \"serif\",\n",
    "                                                   'color': ps_c})\n",
    "                markdown_list.append(pn.Column(markdown, css_classes=['box1'], margin=5, align='center'))\n",
    "                thresh_pct = round(float(bm_dat_tmp.where(bm_dat_tmp < 500).count()['Biomass'])/\n",
    "                                   float(bm_dat_tmp.count()['Biomass']) * 100, 0)\n",
    "                thresh_text = pn.pane.Markdown(f'**{thresh_pct}%** of the region is estimated to have biomass ' +\n",
    "                                               f'less than {500} kg/ha.',\n",
    "                                               style={'font-family': \"Helvetica\"})\n",
    "                thresh_list.append(pn.Column(bm_hist_tmp * hv.VLine(x=500).opts(line_color='black', height=200),\n",
    "                                             thresh_text,\n",
    "                                             css_classes=['box1'], margin=5, align='center'))\n",
    "                stats_rows.append(pn.Row(markdown_list[idx], bm_list[idx], thresh_list[idx], cov_list[idx]))\n",
    "            self.stats_titles = pn.Column(*markdown_list)\n",
    "            self.cov_stats = pn.Column(*cov_list)\n",
    "            self.bm_stats = pn.Column(*bm_list)\n",
    "            self.thresh_stats = pn.Column(*thresh_list)\n",
    "            return pn.Column(*stats_rows)\n",
    "        else:\n",
    "            return pn.Row(None)\n",
    "    \n",
    "    def _create_view(self):\n",
    "        layout = pn.Column(pn.Row(pn.Column(self.username_input,\n",
    "                                  self.password_input,\n",
    "                                  pn.Param(self.param,\n",
    "                                           widgets={'action': pn.widgets.Button(name='Get data', width=200),\n",
    "                                                    'get_stats': pn.widgets.Button(name='Calculate static stats', width=200),\n",
    "                                                   'get_stats_ts': pn.widgets.Button(name='Calculate time-series stats', width=200)},\n",
    "                                  show_name=False),\n",
    "                         self.proc_data,\n",
    "                         self.date), pn.Column(self.create_maps)),\n",
    "                                               self.show_stats)\n",
    "        return layout\n",
    "                                        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "app = getData()\n",
    "app.view.servable()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from src.hls_funcs.masks import shp2mask\n",
    "from src.hls_funcs.smooth import despike_ts_xr, smooth_xr\n",
    "polymask_list = []\n",
    "for idx, ps_c in enumerate(app.poly_stream.data['line_color']):\n",
    "    row, col = polygon(app.poly_stream.data['xs'][idx]*-1, app.poly_stream.data['ys'][idx])\n",
    "    row = row * -1.0\n",
    "    r_map, c_map = app.mapProj(row, col, inverse=True)\n",
    "    r_dat, c_dat = app.datProj(r_map, c_map, inverse=False)\n",
    "    geometries = {\n",
    "        \"type\": \"Polygon\",\n",
    "        \"coordinates\": [list(map(list, zip(r_dat, c_dat)))]\n",
    "    }\n",
    "    ta = affine.Affine(30.0, 0.0, float(app.da['x'].min()),\n",
    "                       0.0, -30.0, float(app.da['y'].max()))\n",
    "    polymask_tmp = shp2mask([geometries], app.da,\n",
    "                         transform=ta, outshape=app.da['BLUE'].isel(date=0).shape, default_value=1)\n",
    "    polymask_tmp['region'] = idx + 1\n",
    "    polymask_tmp['line_color'] = ps_c\n",
    "    polymask_list.append(polymask_tmp)\n",
    "polymask = xr.concat(polymask_list, dim='region')\n",
    "polymask\n",
    "# don't want to fetch data twice, so need to get all masks at the same time with labels, then mean over date+label"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "da_sel_ts = app.da.sel(date=app.da['eo:cloud_cover'] < 80).where(polymask == 1).persist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "date_range_full = np.arange(da_sel_ts['date'].min().values, (da_sel_ts['date'].max() + pd.Timedelta(days=1)).values, timedelta(days=1))\n",
    "ndvi_ts_tmp = ndvi_func(da_sel_ts).chunk({'date': -1, 'y': 250, 'x': 250, 'region': -1})\n",
    "ndvi_ts = xr.DataArray(data=np.ones((len(date_range_full), len(da_sel_ts['y']), len(da_sel_ts['x']), len(da_sel_ts['region']))) * np.nan,\n",
    "                       dims=['date', 'y', 'x', 'region'],\n",
    "                       coords={'date': date_range_full,\n",
    "            'y': da_sel_ts['y'],\n",
    "            'x': da_sel_ts['x'],\n",
    "           'region': da_sel_ts['region']}, name='NDVI')\n",
    "ndvi_ts = ndvi_ts.combine_first(ndvi_ts_tmp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi_ts_mean = ndvi_ts.mean(dim=['y', 'x'])\n",
    "ndvi_ts_dspk = despike_ts_xr(ndvi_ts_mean, dat_thresh=0.025, mask_outliers=False, iters=2, dims=['date'])\n",
    "ndvi_ts_line = smooth_xr(ndvi_ts_dspk, dims=['date']).hvplot.line(x='date', y='NDVI', by='region', line_color=['#1b9e77', '#d95f02'])\n",
    "ndvi_ts_point = ndvi_ts_dspk.hvplot.scatter(x='date', y='NDVI', by='region', color=['#1b9e77', '#d95f02'])\n",
    "ndvi_ts_line * ndvi_ts_point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bm_ts = xr.DataArray(data=np.ones((len(date_range_full), len(da_sel_ts['y']), len(da_sel_ts['x']), len(da_sel_ts['region']))) * np.nan,\n",
    "                       dims=['date', 'y', 'x', 'region'],\n",
    "                       coords={'date': np.arange(da_sel_ts['date'].min().values, (da_sel_ts['date'].max() + pd.Timedelta(days=1)).values, timedelta(days=1)),\n",
    "            'y': da_sel_ts['y'],\n",
    "            'x': da_sel_ts['x'],\n",
    "           'region': da_sel_ts['region']}, name='Biomass')\n",
    "bm_ts_tmp = da_sel_ts.stack(tmp=('date', 'region')).map_blocks(pred_bm, template=da_sel_ts['BLUE'].stack(tmp=('date', 'region')),\n",
    "                                     kwargs=dict(model=app.bm_mod))\n",
    "bm_ts_tmp = bm_ts_tmp.unstack('tmp').transpose('date', 'y', 'x', 'region').chunk({'date': -1, 'y': 250, 'x': 250, 'region': -1})\n",
    "bm_ts_tmp = bm_ts_tmp.where(bm_ts_tmp > 0)\n",
    "bm_ts = bm_ts.combine_first(bm_ts_tmp)\n",
    "bm_ts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bm_ts_dspk = despike_ts_xr(bm_ts.mean(dim=['y', 'x']), dat_thresh=150, mask_outliers=False, iters=2, dims=['date'])\n",
    "bm_ts_line = smooth_xr(bm_ts_dspk, dims=['date']).hvplot.line(x='date', y='Biomass', by='region', line_color=['#1b9e77', '#d95f02'])\n",
    "bm_ts_point = bm_ts_dspk.hvplot.scatter(x='date', y='Biomass', by='region', color=['#1b9e77', '#d95f02'])\n",
    "bm_ts_line * bm_ts_point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lut['LS8T1_SR'] =  {'SR_B1.TIF': 'COASTAL-AEROSOL',\n",
    "               'SR_B2.TIF': 'BLUE', \n",
    "        'SR_B3.TIF': 'GREEN', \n",
    "        'SR_B4.TIF': 'RED', \n",
    "        'SR_B5.TIF': 'NIR1',\n",
    "        'SR_B6.TIF': 'SWIR1',\n",
    "        'SR_B7.TIF': 'SWIR2', \n",
    "        'QA_PIXEL.TIF': 'FMASK'}\n",
    "\n",
    "lut['LS7T1_SR'] =  {'SR_B1.TIF': 'BLUE', \n",
    "        'SR_B2.TIF': 'GREEN', \n",
    "        'SR_B3.TIF': 'RED', \n",
    "        'SR_B4.TIF': 'NIR1',\n",
    "        'SR_B5.TIF': 'SWIR1',\n",
    "        'SR_B7.TIF': 'SWIR2', \n",
    "        'QA_PIXEL.TIF': 'FMASK'}\n",
    "\n",
    "def fetch_sat_api(query):\n",
    "    '''\n",
    "    Queries the sat-api (STAC) backend.\n",
    "    This function handles pagination.\n",
    "    query is a python dictionary to pass as json to the request.\n",
    "    '''\n",
    "    headers = {\n",
    "        \"Content-Type\": \"application/json\",\n",
    "        \"Accept-Encoding\": \"gzip\",\n",
    "        \"Accept\": \"application/geo+json\",\n",
    "    }\n",
    "\n",
    "    url = f\"https://landsatlook.usgs.gov/sat-api/stac/search\"\n",
    "    data = r.post(url, headers=headers, json=query).json()\n",
    "    error = data.get(\"message\", \"\")\n",
    "    if error:\n",
    "        raise Exception(f\"SAT-API failed and returned: {error}\")\n",
    "\n",
    "    meta = data.get(\"meta\", {})\n",
    "    if not meta.get(\"found\"):\n",
    "        return []\n",
    "    print(meta)\n",
    "\n",
    "    features = data[\"features\"]\n",
    "    if data[\"links\"]:\n",
    "        curr_page = int(meta[\"page\"])\n",
    "        query[\"page\"] = curr_page + 1\n",
    "        query[\"limit\"] = meta[\"limit\"]\n",
    "\n",
    "        features = list(itertools.chain(features, fetch_sat_api(query)))\n",
    "\n",
    "    return features\n",
    "\n",
    "\n",
    "query = {\"bbox\":[-104.79107047,   40.78311181, -104.67687336,   40.87008987],\n",
    "         \"query\":{\"collections\":[\"landsat-c2l2-sr\"],\n",
    "                  \"eo:cloud_cover\":{\"lte\":75},\n",
    "                  \"platform\":{\"in\":[\"LANDSAT_7\", \"LANDSAT_8\"]},\n",
    "                  \"landsat:collection_category\":{\"in\":[\"T1\"]}},\n",
    "         \"time\":\"2018-08-21T00:00:00.000/2019-08-21T23:59:59.999\",\"limit\":700}\n",
    "\n",
    "features = fetch_sat_api(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for stac in features:\n",
    "            for band in stac['assets']:\n",
    "                stac['assets'][band]['href'] = stac['assets'][band]['href'].replace('https://landsatlook.usgs.gov/data/', \n",
    "                                                                                    's3://usgs-landsat/')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add AWS_REQUEST_PAYER='requester' to env?\n",
    "import boto3\n",
    "import rasterio as rio\n",
    "from rasterio.session import AWSSession\n",
    "gdal_env=stackstac.DEFAULT_GDAL_ENV.updated(always=dict(session=AWSSession(boto3.Session(), requester_pays=True, region_name='us_west-2')))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls8_stack = stackstac.stack([f for f in features if f['properties']['platform'] == 'LANDSAT_8'], \n",
    "                           epsg=32613, resolution=30, \n",
    "                           assets=[i for i in lut['LS8T1_SR'] if lut['LS8T1_SR'][i] in needed_bands],\n",
    "                           chunksize=(4000, 4000)) #, reader=stackstac.reader_protocol.FakeReader\n",
    "\n",
    "ls8_stack['band'] = [lut['LS8T1_SR'][b] for b in ls8_stack['band'].values]\n",
    "ls8_stack['time'] = [datetime.fromtimestamp(t) for t in ls8_stack.time.astype('int').values//1000000000]\n",
    "ls8_stack = ls8_stack.to_dataset(dim='band')#.reset_coords(['end_datetime', 'start_datetime'], drop=True)\n",
    "\n",
    "\n",
    "ls7_stack = stackstac.stack([f for f in features if f['properties']['platform'] == 'LANDSAT_7'], \n",
    "                           epsg=32613, resolution=30, \n",
    "                           assets=[i for i in lut['LS7T1_SR'] if lut['LS7T1_SR'][i] in needed_bands],\n",
    "                           chunksize=(4000, 4000)) #, reader=stackstac.reader_protocol.FakeReader\n",
    "\n",
    "ls7_stack['band'] = [lut['LS7T1_SR'][b] for b in ls7_stack['band'].values]\n",
    "ls7_stack['time'] = [datetime.fromtimestamp(t) for t in ls7_stack.time.astype('int').values//1000000000]\n",
    "ls7_stack = ls7_stack.to_dataset(dim='band')#.reset_coords(['end_datetime', 'start_datetime'], drop=True)\n",
    "\n",
    "\n",
    "if ls8_stack is not None and ls7_stack is not None:\n",
    "    ls_stack = xr.concat([ls8_stack, ls7_stack], dim='time')\n",
    "elif ls8_stack is not None:\n",
    "    ls_stack = ls8_stack\n",
    "elif l30_stack is not None:\n",
    "    ls_stack = ls7_stack\n",
    "\n",
    "ls_stack['time'] = ls_stack.time.dt.floor(\"D\")\n",
    "ls_stack = ls_stack.rename(dict(time='date'))\n",
    "ls_stack"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_stack = ls_stack.where(polymask == 1)\n",
    "#ls_stack = ls_stack.coarsen(date=7, boundary='pad').mean()\n",
    "ndvi_lta = ndvi_func(ls_stack).chunk({'date': -1, 'y': 250, 'x': 250, 'region': -1})\n",
    "ndvi_lta.name = 'NDVI'\n",
    "ndvi_wkly = ndvi_lta.groupby(ndvi_lta.date.dt.weekofyear).mean()\n",
    "#ndvi_wkly\n",
    "ndvi_wkly.hvplot.scatter(x='week', y='NDVI', by='region', color=['#1b9e77', '#d95f02'])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}