HLS_Shortgrass_Explorer_v8.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()\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\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, pred_bm_thresh\n",
    "import cartopy.crs as ccrs\n",
    "from bokeh.models.formatters import PrintfTickFormatter\n",
    "import stackstac\n",
    "from subprocess import Popen, DEVNULL, STDOUT\n",
    "from getpass import getpass\n",
    "from sys import platform"
   ]
  },
  {
   "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",
    "bm_mod = pickle.load(open('src/models/CPER_HLS_to_VOR_biomass_model_lr_simp.pk', 'rb'))\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",
    "da_bm = pred_bm(da, bm_mod)\n",
    "da_bm = da_bm.where(da_bm > 0)\n",
    "da_se = pred_bm_se(da, bm_mod)\n",
    "da_se = da_se.where(da_bm > 0)\n",
    "t0 = time.time()\n",
    "print(da_bm.isel(time=1).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, 15).date()), str(datetime.now().date())]}\n",
    "tmp_data = get_hls(['spkearney', '1mrChamu'], hls_data=data_dict, aws=True)\n",
    "bm_mod = pickle.load(open('src/models/CPER_HLS_to_VOR_biomass_model_lr_simp.pk', 'rb'))\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",
    "da_bm = pred_bm(da, bm_mod)\n",
    "da_bm = da_bm.where(da_bm > 0)\n",
    "da_se = pred_bm_se(da, bm_mod)\n",
    "da_se = da_se.where(da_bm > 0)\n",
    "t0 = time.time()\n",
    "print(da_bm.isel(time=1).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 src.objects.charts import gauge_obj\n",
    "from holoviews import streams\n",
    "import affine\n",
    "from src.hls_funcs.masks import shp2mask\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",
    "    \n",
    "    action = pm.Action(lambda x: x.param.trigger('action'), label='Load Data and Run Analysis')\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",
    "    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=True, project=False, dynamic=True)\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",
    "    poly_opts = dict(fill_color=['', ''], fill_alpha=[0.0, 0.0], line_color=['#1b9e77', '#d95f02'],\n",
    "                 line_width=[3, 3])\n",
    "    \n",
    "    date = pn.widgets.DatePicker(name='Calendar', width=200)\n",
    "    \n",
    "    tiles = gv.tile_sources.EsriImagery.opts(projection=mapCRS, backend='bokeh')\n",
    "    \n",
    "    polys = hv.Polygons([])\n",
    "    \n",
    "    def __init__(self, **params):\n",
    "        super(getData, self).__init__(**params)\n",
    "        self.da = None\n",
    "        self.da_cov = None\n",
    "        self.da_bm = None\n",
    "        self.da_bm_se = None\n",
    "        self.da_thresh = None\n",
    "        \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",
    "    \n",
    "    @pm.depends('action')\n",
    "    def proc_data(self):\n",
    "        message = 'Not yet launched'\n",
    "        if self.username_input.value != '':\n",
    "            try:\n",
    "                tmp_data = get_hls([self.username_input.value,self.password_input.value], \n",
    "                                   hls_data=self.data_dict, aws=True)\n",
    "                self.da = tmp_data.loc[dict(x=slice(517587.0, 527283.0), y=slice(4524402.0, 4514699.0))]\n",
    "                self.da['time'] = self.da.time.dt.floor(\"D\")\n",
    "                self.da = self.da.rename(dict(time='date'))\n",
    "                da_mask = mask_hls(self.da['FMASK'])\n",
    "                self.da = self.da.where(da_mask == 0)\n",
    "                #self.da_bm = pred_bm(self.da, self.bm_mod)\n",
    "                #self.da_bm = self.da_bm.where(self.da_bm > 0)\n",
    "                #self.da_bm.name = 'Biomass'\n",
    "                #self.da_cov = pred_cov(self.da, self.cov_mod)\n",
    "                #self.da_cov = self.da_cov[['SD', 'GREEN', 'BARE']].to_array(dim='type')\n",
    "                #self.da_cov = self.da_cov.where((self.da_cov < 1.0) | (self.da_cov.isnull()), 1.0)\n",
    "                #self.da_cov.name = 'Cover'\n",
    "                #self.da_bm_se = pred_bm_se(self.da, self.bm_mod)\n",
    "                #self.da_bm_se = self.da_bm_se.where(self.da_bm.notnull())\n",
    "                #self.da_thresh = pred_bm_thresh(self.da_bm, self.da_se, self.thresh_picker.value)\n",
    "                #self.da_thresh.name = 'Threshold'\n",
    "                self.date.enabled_dates = [pd.Timestamp(x).to_pydatetime().date() for x in self.da['date'].values]\n",
    "                message = 'Success!'\n",
    "                return message\n",
    "            except:\n",
    "                return message + ': App Failure'\n",
    "        else:\n",
    "            return message\n",
    "    \n",
    "    #def __init__(self, **params):\n",
    "    #    super(getMaps, self).__init__(**params)\n",
    "        #self.date_init = self.da['date'][0].values\n",
    "        #self.date.default=pd.Timestamp(self.date_init).to_pydatetime()\n",
    "      \n",
    "    @pm.depends('date.param')\n",
    "    def create_maps(self):\n",
    "        if self.da is not None and self.date.value is not None:\n",
    "            date_idx = self.date.enabled_dates.index(self.date.value)\n",
    "            da_sel = self.da.isel(date=date_idx).persist()\n",
    "            da_bm = pred_bm(da_sel, self.bm_mod)\n",
    "            da_bm = da_bm.where(da_bm > 0)\n",
    "            da_bm.name = 'Biomass'\n",
    "            da_cov = pred_cov(da_sel, 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",
    "            da_bm_se = pred_bm_se(da_sel, bm_mod)\n",
    "            da_bm_se = da_bm_se.where(da_bm.notnull())\n",
    "            da_thresh = pred_bm_thresh(da_bm, da_bm_se, self.thresh_picker.value)\n",
    "            da_thresh.name = 'Threshold'\n",
    "            cov_map = da_cov.hvplot.rgb(x='x', y='y', bands='type',\n",
    "                                        **self.map_args).opts(**self.map_opts)\n",
    "            bm_map = da_bm.hvplot(x='x', y='y',\n",
    "                                  cmap='Viridis', clim=(100, 1000), colorbar=False,\n",
    "                                  **self.map_args).opts(**self.map_opts)\n",
    "            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.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",
    "            \n",
    "            all_maps = pn.Tabs(('Cover', self.tiles * cov_map * self.polys),\n",
    "                       ('Biomass', self.tiles * bm_map * self.polys),\n",
    "                       ('Threshold', self.tiles * thresh_map * self.polys))\n",
    "        else:\n",
    "            all_maps = None\n",
    "        return all_maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_step = getData()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "layout = pn.Row(pn.Column(data_step.username_input,\n",
    "                          data_step.password_input,\n",
    "                          pn.Param(data_step.param,\n",
    "                                   widgets={'action': pn.widgets.Button(name='Compute',\n",
    "                                                                           width=200)}),\n",
    "                         data_step.proc_data,\n",
    "                         data_step.date), data_step.create_maps)\n",
    "layout.servable()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t0 = time.time()\n",
    "#print(data_step.da_bm.sel(date=np.datetime64(data_step.date.value)).values)\n",
    "print(data_step.da_bm.sel(date=np.datetime64(datetime(2021, 1, 18).date())).values)\n",
    "#print(data_step.da_bm.isel(date=data_step.date.enabled_dates.index(data_step.date.value)).values)\n",
    "#print(data_step.da_bm.isel(date=6).values)\n",
    "t1 = time.time()\n",
    "print(t1-t0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "class App(pm.Parameterized):\n",
    "    date_init = dat['date'][0].values\n",
    "    thresh_init = 500\n",
    "    bm_mean1 = 0\n",
    "\n",
    "    da_bm = pred_bm(dat, bm_mod)\n",
    "    da_bm.name = 'Biomass'\n",
    "    da_bm = da_bm.where(da_bm > 0)\n",
    "    da_bm_sel = da_bm.sel(date=date_init).persist()\n",
    "    da_cov = pred_cov(dat, 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_sel = da_cov.sel(date=date_init).persist()\n",
    "    da_se = pred_bm_se(dat, bm_mod)\n",
    "    da_se = da_se.where(da_bm.notnull())\n",
    "    da_thresh = pred_bm_thresh(da_bm, da_se, thresh_init)\n",
    "    da_thresh.name = 'Threshold'\n",
    "    da_thresh_sel = da_thresh.sel(date=date_init).persist()\n",
    "\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=True, project=False, dynamic=True)\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",
    "    poly_opts = dict(fill_color=['', ''], fill_alpha=[0.0, 0.0], line_color=['#1b9e77', '#d95f02'],\n",
    "                     line_width=[3, 3])\n",
    "    gauge_opts = dict(height=200, width=300)\n",
    "\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",
    "    pn.extension(raw_css=[css])\n",
    "\n",
    "    basemap = pm.ObjectSelector(default=\"Satellite\", objects=[\"Satellite\", \"Map\"])\n",
    "    alpha = pm.Number(default=1.0)\n",
    "    date = pm.CalendarDate(default=pd.Timestamp(date_init).to_pydatetime())\n",
    "    thresh = pm.Integer(default=thresh_init)\n",
    "    action = pm.Action(lambda self: self.param.trigger('action'), 'Compute')\n",
    "\n",
    "    #action = pn.widgets.Button(name='Save regions and \\ncompute stats',\n",
    "    #                           width=200)\n",
    "    #action.on_click(lambda self: self.param.trigger('action'))\n",
    "\n",
    "    #action_val = action\n",
    "\n",
    "\n",
    "    cov_map = da_cov_sel.hvplot.rgb(x='x', y='y', bands='type',\n",
    "                                              **map_args).opts(**map_opts)\n",
    "    bm_map = da_bm_sel.hvplot(x='x', y='y',\n",
    "                                        cmap='Viridis', clim=(100, 1000), colorbar=False,\n",
    "                                        **map_args).opts(**map_opts)\n",
    "    thresh_map = da_thresh_sel.hvplot(x='x', y='y',\n",
    "                                                cmap='YlOrRd', clim=(0.05, 0.95), colorbar=False,\n",
    "                                                **map_args).opts(**map_opts)\n",
    "\n",
    "    tiles = gv.tile_sources.EsriImagery.opts(projection=mapCRS, backend='bokeh')\n",
    "\n",
    "    polys = hv.Polygons([])\n",
    "\n",
    "    def __init__(self, **params):\n",
    "        super(App, self).__init__(**params)\n",
    "\n",
    "        self.gauge_obj = deepcopy(gauge_obj)\n",
    "\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",
    "        self.select_stream = streams.Selection1D(source=self.polys)\n",
    "\n",
    "        self.startX, self.endX = (float(self.da_bm['x'].min().values), float(self.da_bm['x'].max().values))\n",
    "        self.startY, self.endY = (float(self.da_bm['y'].min().values), float(self.da_bm['y'].max().values))\n",
    "        self.cov_stats = ''\n",
    "        self.bm_stats = ''\n",
    "        self.thresh_stats = ''\n",
    "        #self.stats_dict = {0: self.cov_stats,\n",
    "        #         1: self.bm_stats,\n",
    "        #         2: self.thresh_stats}\n",
    "        #self.stats = self.stats_dict[self.active_tab]\n",
    "        #self.text2.jscallback(args={'gauge1': self.gauge_pane1}, code=\"\"\"\n",
    "        #gauge1.data.series[0].data[0].value = cb_obj.value\n",
    "        #gauge1.properties.data.change.emit()\n",
    "        #\"\"\")\n",
    "        self.all_maps = pn.Tabs(('Cover', pn.Row(self.tiles * self.cov_map * self.polys, self.cov_stats)),\n",
    "                           ('Biomass', pn.Row(self.tiles * self.bm_map * self.polys, self.bm_stats)),\n",
    "                           ('Threshold', pn.Row(self.tiles * self.thresh_map * self.polys, self.thresh_stats)))\n",
    "        self.active_tab = self.all_maps.active\n",
    "\n",
    "    def keep_zoom(self, x_range, y_range):\n",
    "        map_x_range, map_y_range = self.mapProj(x_range, y_range, inverse=True)\n",
    "        (self.startX, self.endX), (self.startY, self.endY) = self.datProj(map_x_range, map_y_range, inverse=False)\n",
    "\n",
    "    @pm.depends('basemap')\n",
    "    def map_base(self):\n",
    "        if self.basemap == \"Satellite\":\n",
    "            self.tiles = gv.tile_sources.EsriImagery(projection=self.mapCRS, backend='bokeh')\n",
    "        elif self.basemap == \"Map\":\n",
    "            self.tiles = gv.tile_sources.Wikipedia(projection=self.mapCRS, backend='bokeh')\n",
    "        return self.tiles\n",
    "\n",
    "    @pm.depends('date', watch=True)\n",
    "    def bm_date(self):\n",
    "        self.date_init = np.datetime64(self.date)\n",
    "        self.da_bm_sel = self.da_bm.sel(date=self.date_init).persist()\n",
    "        self.bm_map = self.da_bm_sel.hvplot(x='x', y='y',\n",
    "                                                               #xlim=(self.startX, self.endX),\n",
    "                                                               #ylim=(self.startY, self.endY),\n",
    "                                                               cmap='Viridis', clim=(100, 1000),\n",
    "                                                               colorbar=False,\n",
    "                                                               **self.map_args).opts(alpha=self.alpha,\n",
    "                                                                                     **self.map_opts)\n",
    "        self.bm_map.streams[-1].add_subscriber(self.keep_zoom)\n",
    "        return self.bm_map\n",
    "\n",
    "    @pm.depends('date', watch=True)\n",
    "    def cov_date(self):\n",
    "        self.date_init = np.datetime64(self.date)\n",
    "        self.da_cov_sel = self.da_cov.sel(date=self.date_init).persist()\n",
    "        self.cov_map = self.da_cov_sel.hvplot.rgb(x='x', y='y',\n",
    "                                                    #xlim=(self.startX, self.endX),\n",
    "                                                    #ylim=(self.startY, self.endY),\n",
    "                                                    bands='type',\n",
    "                                                    **self.map_args).opts(alpha=self.alpha,\n",
    "                                                                          **self.map_opts)\n",
    "        return self.cov_map\n",
    "\n",
    "    @pm.depends('date', 'thresh', watch=True)\n",
    "    def thresh_date(self):\n",
    "        self.date_init = np.datetime64(self.date)\n",
    "        self.thresh_init = self.thresh\n",
    "        self.da_thresh = pred_bm_thresh(self.da_bm, self.da_se, self.thresh_init)\n",
    "        self.da_thresh_sel = self.da_thresh.sel(date=self.date_init).persist()\n",
    "        self.thresh_map = self.da_thresh_sel.hvplot(x='x', y='y',\n",
    "                                           #xlim=(self.startX, self.endX),\n",
    "                                           #ylim=(self.startY, self.endY),\n",
    "                                           cmap='YlOrRd', clim=(0.05, 0.95), colorbar=False,\n",
    "                                           **self.map_args).opts(alpha=self.alpha,\n",
    "                                                                 **self.map_opts)\n",
    "        self.bm_map.streams[-1].add_subscriber(self.keep_zoom)\n",
    "        return self.thresh_map\n",
    "\n",
    "    @pm.depends('alpha', watch=True)\n",
    "    def cov_alpha(self):\n",
    "        return self.cov_map.opts(alpha=self.alpha,\n",
    "                                 xlim=self.bm_map.streams[-1].x_range,\n",
    "                                 ylim=self.bm_map.streams[-1].y_range,\n",
    "                               **self.map_opts)\n",
    "    @pm.depends('alpha', watch=True)\n",
    "    def bm_alpha(self):\n",
    "        self.bm_map = self.bm_map.opts(alpha=self.alpha,\n",
    "                                 xlim=self.bm_map.streams[-1].x_range,\n",
    "                                 ylim=self.bm_map.streams[-1].y_range,\n",
    "                               **self.map_opts)\n",
    "        self.bm_map.streams[-1].add_subscriber(self.keep_zoom)\n",
    "        return self.bm_map\n",
    "\n",
    "    @pm.depends('alpha', watch=True)\n",
    "    def thresh_alpha(self):\n",
    "        return self.thresh_map.opts(alpha=self.alpha,\n",
    "                                 xlim=self.bm_map.streams[-1].x_range,\n",
    "                                 ylim=self.bm_map.streams[-1].y_range,\n",
    "                               **self.map_opts)\n",
    "\n",
    "    @pm.depends('action', watch=True)\n",
    "    def show_hist(self):\n",
    "        if self.poly_stream.data is None:\n",
    "            self.cov_stats = ''\n",
    "            self.bm_stats = ''\n",
    "            self.thresh_stats = ''\n",
    "        else:\n",
    "            self.bm_stats = 'Yes'\n",
    "            thresh_list = []\n",
    "            bm_list = []\n",
    "            cov_list = []\n",
    "            ts_yr_list = []\n",
    "            ts_avg_list = []\n",
    "            for idx, ps_c in enumerate(self.poly_stream.data['line_color']):\n",
    "                xs_map, ys_map = self.mapProj(self.poly_stream.data['xs'][idx],\n",
    "                                              self.poly_stream.data['ys'][idx], inverse=True)\n",
    "                xs_dat, ys_dat = self.datProj(xs_map, ys_map, inverse=False)\n",
    "                geometries = {\n",
    "                    \"type\": \"Polygon\",\n",
    "                    \"coordinates\": [\n",
    "                        list(map(list, zip(xs_dat, ys_dat)))\n",
    "                    ]\n",
    "                }\n",
    "                ta = affine.Affine(30.0, 0.0, float(self.da_bm_sel['x'].min()),\n",
    "                                   0.0, -30.0, float(self.da_bm_sel['y'].max()))\n",
    "                poly_mask = shp2mask([geometries], self.da_bm_sel,\n",
    "                                     transform=ta, outshape=self.da_bm_sel.shape, default_value=1)\n",
    "                da_bm_tmp = self.da_bm_sel.where(poly_mask == 1)\n",
    "                bm_hist_tmp = da_bm_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).options(toolbar=None)\n",
    "                markdown = pn.pane.Markdown('## Region stats', height=50,\n",
    "                                            style={'font-family': \"serif\",\n",
    "                                                   'color': ps_c})\n",
    "                thresh_pct = round(float(da_bm_tmp.where(da_bm_tmp < self.thresh_init).count())/\n",
    "                                   float(da_bm_tmp.count()) * 100, 0)\n",
    "                thresh_text = pn.pane.Markdown(f'**{thresh_pct}%** of the region is estimated to have biomass ' +\n",
    "                                               f'less than {self.thresh_init} kg/ha.',\n",
    "                                               style={'font-family': \"Helvetica\"})\n",
    "                thresh_list.append(pn.Column(pn.Row(pn.layout.HSpacer(), markdown, pn.layout.HSpacer()),\n",
    "                                             bm_hist_tmp * hv.VLine(x=self.thresh_init).opts(line_color='black'),\n",
    "                                             thresh_text,\n",
    "                                             css_classes=['box1'], margin=5))\n",
    "                bm_gauge_obj = deepcopy(self.gauge_obj)\n",
    "                bm_gauge_obj['series'][0]['data'][0]['value'] = int(da_bm_tmp.mean().values)\n",
    "                bm_gauge_pane = pn.pane.ECharts(bm_gauge_obj, **self.gauge_opts)\n",
    "                bm_list.append(pn.Column(pn.Row(pn.layout.HSpacer(),markdown, pn.layout.HSpacer()),\n",
    "                                         bm_gauge_pane,\n",
    "                                         css_classes=['box1'], margin=5))\n",
    "                #yr = int(self.da_bm_sel.YEAR.values)\n",
    "                #ts_bm_yr_tmp = self.da_bm.where(poly_mask == 1).sel(date=slice(datetime(yr, 5, 1),\n",
    "                #                                                                datetime(yr, 10, 31)))\n",
    "                #ts_yr_list.append(ts_bm_yr_tmp.mean(dim=['x', 'y']).hvplot.line(x='date',\n",
    "                #                                                               y='Biomass').opts(line_color=ps_c))\n",
    "                #ts_bm_avg_tmp = self.da_bm.where(poly_mask == 1).groupby(da_bm.date.dt.dayofyear).mean(dim=['x', 'y'])\n",
    "                #ts_avg_list.append(ts_bm_avg_tmp.hvplot.line(x='date', y='Biomass').opts(line_color=ps_c))\n",
    "                da_cov_tmp = self.da_cov_sel.where(poly_mask == 1)\n",
    "                cov_factors = list(da_cov_tmp.type.values)\n",
    "                cov_vals = [round(float(da_cov_tmp.sel(type=f).mean().values), 2) for f in cov_factors]\n",
    "                from bokeh.models import NumeralTickFormatter\n",
    "                pct_fmt = NumeralTickFormatter(format=\"0%\")\n",
    "                cov_colors = hv.Cycle(['red', 'green', 'blue'])\n",
    "                cov_scatter_tmp = hv.Overlay([hv.Scatter(f) for f in list(zip(cov_factors, 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(pn.Row(pn.layout.HSpacer(), markdown, pn.layout.HSpacer()),\n",
    "                                          (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))\n",
    "\n",
    "\n",
    "            self.polys=self.poly_stream.element.opts(xlim=(self.startX, self.endX),\n",
    "                                                     ylim=(self.startY, self.endY))\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",
    "            #self.gauge_pane1.object['series'][0]['data'][0]['value'] = int(da_bm_tmp.mean().values)\n",
    "            #self.gauge_pane2.object['series'][0]['data'][0]['value'] = int(da_bm_tmp.mean().values) + 100\n",
    "            self.thresh_stats = pn.Column(*thresh_list)\n",
    "            self.bm_stats = pn.Column(*bm_list)\n",
    "            self.cov_stats = pn.Column(*cov_list)\n",
    "\n",
    "\n",
    "    def view_all(self):\n",
    "        #self.da_bm_sel.name = 'Biomass'\n",
    "        self.active_tab = self.all_maps.active\n",
    "        self.bm_map.streams[-1].add_subscriber(self.keep_zoom)\n",
    "        #self.stats = self.stats_dict[self.active_tab]\n",
    "        base = hv.DynamicMap(self.map_base)\n",
    "        cov = self.cov_map\n",
    "        bm = self.bm_map\n",
    "        thresh = self.thresh_map\n",
    "        self.all_maps = pn.Tabs(('Cover', pn.Row(base * cov * self.polys, self.cov_stats)),\n",
    "                                ('Biomass', pn.Row(base * bm * self.polys, self.bm_stats)),\n",
    "                                ('Threshold', pn.Row(base * thresh * self.polys, self.thresh_stats)),\n",
    "                                active=self.active_tab)\n",
    "        return pn.Column(self.all_maps)\n",
    "\n",
    "\n",
    "\n",
    "viewer = App()\n",
    "layout = pn.Row(pn.Column(pn.Param(viewer.param,\n",
    "                                      widgets={'date': pn.widgets.DatePicker(name='Calendar',\n",
    "                                                                             #value=datetime(2021, 1, 1).date(),\n",
    "                                                                             enabled_dates = [pd.Timestamp(x).to_pydatetime().date() for x in dat['date'].values],\n",
    "                                                                             width=200),\n",
    "                                               'alpha': pn.widgets.FloatSlider(name='Map transparency',\n",
    "                                                                               value=1.0,\n",
    "                                                                               start=0.0, end=1.0,\n",
    "                                                                               step=0.1,\n",
    "                                                                               width=200),\n",
    "                                               'thresh': pn.widgets.IntSlider(name='Threshold',\n",
    "                                                                              start=200, end=2000,\n",
    "                                                                              step=25, value=500,\n",
    "                                                                              format=PrintfTickFormatter(\n",
    "                                                                                  format='%d kg/ha'),\n",
    "                                                                              width=200),\n",
    "                                               'basemap': pn.widgets.Select(name=\"Change basemap\",\n",
    "                                                                            options=[\"Satellite\", \"Map\"],\n",
    "                                                                            width=200),\n",
    "                                               'action': pn.widgets.Button(name='Save regions and \\ncompute stats',\n",
    "                                                                           width=200)\n",
    "                                               })),\n",
    "                viewer.view_all)\n",
    "layout.servable()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4519550.5"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "np.mean([4524402.0, 4514699.0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}