hls_biomass_aws_v2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'satsearch'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-2-1936424463ca>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mos\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[1;31m#import fsspec, os, netrc, aiohttp,dask\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 6\u001b[1;33m \u001b[1;32mfrom\u001b[0m \u001b[0msatsearch\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mSearch\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      7\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mhvplot\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpandas\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      8\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mhvplot\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mxarray\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mModuleNotFoundError\u001b[0m: No module named 'satsearch'"
     ]
    }
   ],
   "source": [
    "from dask.distributed import LocalCluster, Client\n",
    "import xarray as xr\n",
    "import dask\n",
    "import intake\n",
    "import os\n",
    "#import fsspec, os, netrc, aiohttp,dask\n",
    "from satsearch import Search\n",
    "import hvplot.pandas\n",
    "import hvplot.xarray\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "import gdal\n",
    "import requests\n",
    "import concurrent.futures\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 time\n",
    "import jinja2 as jj2\n",
    "from rasterio.crs import CRS\n",
    "from tempfile import NamedTemporaryFile\n",
    "import re\n",
    "from datetime import datetime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# AUTHENTICATION CONFIGURATION\n",
    "from netrc import netrc\n",
    "from subprocess import Popen\n",
    "from getpass import getpass\n",
    "\n",
    "urs = 'urs.earthdata.nasa.gov'    # Earthdata URL to call for authentication\n",
    "prompts = ['Enter NASA Earthdata Login Username \\n(or create an account at urs.earthdata.nasa.gov): ',\n",
    "           'Enter NASA Earthdata Login Password: ']\n",
    "\n",
    "# Determine if netrc file exists, and if so, if it includes NASA Earthdata Login Credentials\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(getpass(prompt=prompts[0]), homeDir + os.sep), shell=True)\n",
    "    Popen('echo password {} >> {}.netrc'.format(getpass(prompt=prompts[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(getpass(prompt=prompts[0]), homeDir + os.sep), shell=True)\n",
    "    Popen('echo password {} >> {}.netrc'.format(getpass(prompt=prompts[1]), homeDir + os.sep), shell=True)\n",
    "    del homeDir\n",
    "del urs, prompts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'accessKeyId': 'ASIAZLX6ZES47TNY5DOQ',\n",
       " 'secretAccessKey': 'Z91jrRUNmCFTcpXTe74OS98WxRmEw1SByvnf0pka',\n",
       " 'sessionToken': 'FwoGZXIvYXdzEBgaDPUdJGBFljfs8pOkWyLaAaQCfzu5fjmUtPOJrwrflKHfLp76MgArjLp6MiSpyMwoiG5xzqDMsssJGISuKLrvYAc1EbEXNdhQJEmyE91q/OSBVAs5SZK6LjvU1DCJkP/GEeyNEaZK3B6x4Ap36UKFBsyD5Q8RdH0IEMdZ2PVLQC8WVMpj0A9bHXcxbCOSaAr3KBJTMwk/sykyLpiggS6Uh0OtWWoT++QHl6EgsCBIp55kf5YaJorTSeAfJqtmQIVHi4CXSaYlF0i6kEjLlPSXbiC75qQgQ/xPqV+vFC1sZf107pzP6cO941s2KLin24EGMi1hABLAFBdL873Oazyq59acC1pEl8+oDhNlF9AQUTYpcbpoLmvKM2XZXuryLtA=',\n",
       " 'expiration': '2021-02-24 23:31:20+00:00'}"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s3_cred = requests.get('https://lpdaac.earthdata.nasa.gov/s3credentials').json()\n",
    "s3_cred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Setup GDAL Env for optimum performance\n",
    "env = dict(GDAL_DISABLE_READDIR_ON_OPEN='YES', \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",
    "\n",
    "\n",
    "os.environ.update(env)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a LUT dict including the HLS product bands mapped to names\n",
    "lut = {'HLSS30': \n",
    "       {'COASTAL-AEROSOL':'B01', 'BLUE':'B02', 'GREEN':'B03', 'RED':'B04', \n",
    "        'RED-EDGE1':'B05', 'RED-EDGE2':'B06', 'RED-EDGE3':'B07', 'NIR-Broad':'B08', 'NIR1':'B8A', \n",
    "        'WATER-VAPOR':'B09', 'CIRRUS':'B10', 'SWIR1':'B11', 'SWIR2':'B12', 'FMASK':'Fmask'},\n",
    "       'HLSL30': \n",
    "       {'COASTAL-AEROSOL':'B01', 'BLUE':'B02', 'GREEN':'B03', 'RED':'B04', \n",
    "        'NIR1':'B05', 'SWIR1':'B06','SWIR2':'B07', \n",
    "        'CIRRUS':'B09', 'TIR1':'B10', 'TIR2':'B11', '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']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 items found!\n"
     ]
    }
   ],
   "source": [
    "import requests as r\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_collections = [l['href'] for l in lp_links if l['rel'] == 'collections'][0]  # Set collections endpoint to variable\n",
    "collections_response = r.get(f\"{lp_collections}\").json()    \n",
    "collections = collections_response['collections']\n",
    "hls_collections = [c for c in collections if 'HLS' in c['title']]\n",
    "s30 = [h for h in hls_collections if h['short_name'] == 'HLSS30'][0]  # Grab HLSS30 collection\n",
    "s30_id = s30['id']\n",
    "l30 = [h for h in hls_collections if h['short_name'] == 'HLSL30'][0]  # Grab HLSL30 collection\n",
    "l30_id = l30['id']\n",
    "\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",
    "\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 = '2021-01-01/2021-01-31'  # Define start time period / end time period\n",
    "search_query3 = f\"{search_query2}&datetime={date_time}\"  # Add to query that already includes bbox\n",
    "\n",
    "# Search for the HLSS30 items of interest:\n",
    "search_query4 = f\"{search_query3}&collections={s30_id}\"\n",
    "s30_items = r.get(search_query4).json()['features']\n",
    "\n",
    "# Search for the HLSL30 items of interest:\n",
    "search_query5 = f\"{search_query3}&collections={l30_id}\"\n",
    "l30_items = r.get(search_query5).json()['features']\n",
    "\n",
    "# Combine the S30 ad L30 items:\n",
    "hls_items = s30_items + l30_items\n",
    "\n",
    "print(f\"{len(hls_items)} items found!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Function to build metadata for Stac Catalog\n",
    "def open_hls_meta(stac_id):\n",
    "    var_url = urlopen('https://cmr.earthdata.nasa.gov/search/concepts/'+stac_id)\n",
    "    xmldoc = parse(var_url)\n",
    "    res={stac_id:{}}\n",
    "    for child in xmldoc.findall('.//AdditionalAttributes/AdditionalAttribute'):\n",
    "        if child.find('Name').text in ['ULX',\n",
    "                                       'ULY',\n",
    "                                       'NROWS',\n",
    "                                       'NCOLS',\n",
    "                                       'SPATIAL_RESOLUTION',\n",
    "                                       'HORIZONTAL_CS_CODE',\n",
    "                                       'SENSING_TIME',\n",
    "                                       'MGRS_TILE_ID',\n",
    "                                       'CLOUD_COVERAGE',\n",
    "                                       'FILLVALUE',\n",
    "                                       'QA_FILLVALUE',\n",
    "                                       'MEAN_SUN_AZIMUTH_ANGLE',\n",
    "                                       'MEAN_SUN_ZENITH_ANGLE',\n",
    "                                       'MEAN_VIEW_AZIMUTH_ANGLE',\n",
    "                                       'MEAN_VIEW_ZENITH_ANGLE',\n",
    "                                       'NBAR_SOLAR_ZENITH',\n",
    "                                       'IDENTIFIER_PRODUCT_DOI',\n",
    "                                       'IDENTIFIER_PRODUCT_DOI_AUTHORITY',\n",
    "                                       'REF_SCALE_FACTOR',\n",
    "                                       'ADD_OFFSET']:\n",
    "            res[stac_id][child.find('Name').text]=child.find('.//Value').text\n",
    "    return(res)\n",
    "#Convert STAC catalog to VRT file(s) (1 vrt file per band)\n",
    "#Outputs as single xarray dataset object (dims = x,y,t ; variables = bands)\n",
    "def build_xr(catalog,bands,):\n",
    "    #Retreive Metadata using threads - not cpu bound, so works well\n",
    "    l_meta={}\n",
    "    t1 = time.time()\n",
    "    with concurrent.futures.ThreadPoolExecutor(5) as executor:\n",
    "        futures = []\n",
    "        for stac in catalog:\n",
    "            stac_id = stac['id']\n",
    "            futures.append(executor.submit(open_hls_meta, stac_id=stac_id))\n",
    "        for future in concurrent.futures.as_completed(futures):\n",
    "            l_meta.update(future.result())\n",
    "        for stac in catalog:\n",
    "            stac_id = stac['id']\n",
    "            for band in bands:\n",
    "                if bool(re.search('S30', l_meta[stac_id]['IDENTIFIER_PRODUCT_DOI'])):\n",
    "                        b = lut['HLSS30'][band]\n",
    "                elif bool(re.search('L30', l_meta[stac_id]['IDENTIFIER_PRODUCT_DOI'])):\n",
    "                        b = lut['HLSL30'][band]\n",
    "                l_meta[stac_id][b+'_url'] = '/vsicurl/'+stac['assets'][b]['href']#.replace('https://lpdaac.earthdata.nasa.gov/lp-prod-protected', '/vsis3/lp-prod-protected')\n",
    "    \n",
    "    t2 = time.time()\n",
    "    print('Get File Meta',t2-t1)\n",
    "    \n",
    "    #Setup template file for each raster. Use Jinja to quickly fill in metadata.\n",
    "    vrt_template = jj2.Template('''\n",
    "    <VRTDataset rasterXSize=\"{{rasterXSize}}\" rasterYSize=\"{{rasterYSize}}\">\n",
    "      <SRS>{{SRS}}</SRS>\n",
    "      <GeoTransform>{{GeoTransform}}</GeoTransform>\n",
    "      <VRTRasterBand dataType=\"{{dtype}}\" band=\"1\">\n",
    "        <NoDataValue>{{nodata}}</NoDataValue>\n",
    "        <Scale>{{scale}}</Scale>\n",
    "        <Metadata>\n",
    "          <MDI key=\"obs_date\">{{obs_date}}</MDI>\n",
    "        </Metadata>\n",
    "        <SimpleSource>\n",
    "          <SourceFilename relativeToVRT=\"1\">{{SourceFilename}}</SourceFilename>\n",
    "          <SourceBand>1</SourceBand>\n",
    "          <SourceProperties RasterXSize=\"{{rasterXSize}}\" RasterYSize=\"{{rasterYSize}}\" DataType=\"{{dtype}}\" BlockXSize=\"1024\" BlockYSize=\"1024\" />\n",
    "          <SrcRect xOff=\"0\" yOff=\"0\" xSize=\"{{rasterXSize}}\" ySize=\"{{rasterYSize}}\" />\n",
    "          <DstRect xOff=\"0\" yOff=\"0\" xSize=\"{{rasterXSize}}\" ySize=\"{{rasterYSize}}\" />\n",
    "          <NODATA>{{nodata}}</NODATA>\n",
    "        </SimpleSource>\n",
    "      </VRTRasterBand>\n",
    "    </VRTDataset>\n",
    "    ''')\n",
    "    \n",
    "    #Enumerate the stac catalog (by band) and create vrt objects in a dictionary (l_vrt)\n",
    "    l_vrt={}\n",
    "    for i, band in enumerate(bands):             \n",
    "        l_tmp = []\n",
    "        for k in l_meta.keys():\n",
    "            item = l_meta[k]\n",
    "            if bool(re.search('S30', item['IDENTIFIER_PRODUCT_DOI'])):\n",
    "                b = lut['HLSS30'][band]\n",
    "            elif bool(re.search('L30', item['IDENTIFIER_PRODUCT_DOI'])):\n",
    "                b = lut['HLSL30'][band]  \n",
    "            vrt = vrt_template.render(rasterXSize = int(item['NCOLS']),\n",
    "                                      rasterYSize = int(item['NROWS']),\n",
    "                                      SourceFilename = item[b+'_url'],\n",
    "                                      SRS=CRS.from_epsg('32613').wkt,#CRS.from_epsg(item['HORIZONTAL_CS_CODE'].split(':')[1]).wkt,\n",
    "                                      GeoTransform=item['ULX']+', '+item['SPATIAL_RESOLUTION']+', 0, '+item['ULY']+', 0, -'+item['SPATIAL_RESOLUTION'],\n",
    "                                      band=band,\n",
    "                                      obs_date=item['SENSING_TIME'],\n",
    "                                      dtype='int16',\n",
    "                                      nodata=item['FILLVALUE'],\n",
    "                                      scale=item['REF_SCALE_FACTOR'])\n",
    "            l_tmp.append(vrt)\n",
    "        l_vrt[band] = l_tmp\n",
    "    \n",
    "    t3 = time.time()\n",
    "    print('Build VRT',t3-t2)\n",
    "    \n",
    "    #Use GDAL to merge vrts from the same bands\n",
    "    vrt_files={}\n",
    "    for b in bands:\n",
    "        with NamedTemporaryFile() as tmpfile:\n",
    "            vrt_options = gdal.BuildVRTOptions(separate=True, bandList=[1])\n",
    "            my_vrt = gdal.BuildVRT(tmpfile.name, l_vrt[b], options=vrt_options)\n",
    "            my_vrt = None\n",
    "            f = tmpfile.read().decode(\"utf-8\")\n",
    "            vrt_files[b]=f\n",
    "    \n",
    "    t4 = time.time()\n",
    "    print('Build VRT',t4-t3)\n",
    "    \n",
    "    #Extract metadata and lazy-load a nd merge into single xarray object\n",
    "    if type(bands) == str:\n",
    "        bands = [bands]\n",
    "    b_l = []\n",
    "    for b in bands:\n",
    "        v = vrt_files[b]\n",
    "        xmldoc = fromstring(v)\n",
    "        dates = []\n",
    "        for child in xmldoc.findall('''.//VRTRasterBand/ComplexSource/SourceFilename'''):\n",
    "            dates.append(to_datetime(child.text[child.text.find('obs_date')+10:\n",
    "                                                child.text.find('obs_date')+36]).date())\n",
    "        ds_tmp = xr.open_rasterio(v,chunks={'band':1,\n",
    "                                            'x':'auto',\n",
    "                                            'y':'auto'}).to_dataset(name=b)\n",
    "        ds_tmp = ds_tmp.rename({'band':'time'})\n",
    "        ds_tmp['time'] = to_datetime(dates)\n",
    "        b_l.append(ds_tmp)\n",
    "    \n",
    "    t5 = time.time()\n",
    "    print('Create/Merge Xarray',t5-t4)\n",
    "    \n",
    "    ds = xr.merge(b_l)\n",
    "    return(ds,vrt_files)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dask.config.set({'distributed.dashboard.link':'http://localhost:8888/proxy/8787/status'})#'https://localhost:8787/status'})\n",
    "cluster = LocalCluster(threads_per_worker=2)\n",
    "client = Client(cluster)\n",
    "client"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "da, vrt = build_xr(hls_items,  ['BLUE', 'GREEN', 'RED', 'NIR1', 'SWIR1', 'SWIR2', 'FMASK'])\n",
    "da"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'\\nh = hls_items[0]\\n\\nfrom datetime import datetime\\nimport re\\ndef stac_to_xr(item, band):\\n    if band not in all_bands:\\n        sys.exit(f\"Band: {band} is not a valid input option.\" +\\n                 \"Valid inputs are ALL, COASTAL-AEROSOL, BLUE, GREEN, RED, RED-EDGE1, RED-EDGE2, RED-EDGE3,\" + \\n                 \" NIR1, SWIR1, SWIR2, CIRRUS, TIR1, TIR2, WATER-VAPOR, FMASK.\" + \\n                 \"To request multiple layers, provide them in comma separated format with no spaces.\" +\\n                 \"Unsure of the names for your bands?\" +\\n                 \"--check out the README which contains a table of all bands and band names.\")\\n    if item[\\'collection\\'] == s30_id:\\n        href = item[\\'assets\\'][lut[\\'HLSS30\\'][band]][\\'href\\']\\n    if item[\\'collection\\'] == l30_id:\\n        href = item[\\'assets\\'][lut[\\'HLSL30\\'][band]][\\'href\\']\\n    ds = xr.open_rasterio(re.sub(\\'https://lpdaac.earthdata.nasa.gov/\\', \\'/vsis3/\\', href), chunks={\\'band\\':1,\\n                                            \\'x\\':-1,\\n                                            \\'y\\':-1}).to_dataset(name=band)\\n    ds = ds.rename({\\'band\\': \\'time\\'})\\n    ds[\\'time\\'] = [item[\\'properties\\'][\\'datetime\\']]\\n    ds[\\'time\\'] = ds[\\'time\\'].astype(\\'datetime64\\').dt.floor(\\'D\\')\\n    #ds = ds.assign_coords(t = datetime.strptime(item[\\'properties\\'][\\'datetime\\'], \\'%Y-%m-%dT%H:%M:%S.%fZ\\'))\\n    return ds\\n\\ndef stac_to_ds(items, bands):\\n    ds_list = []\\n    for item in items:\\n        ds_item_list = []\\n        for band in bands:\\n            ds_item_list.append(stac_to_xr(item, band))\\n            ds_item = xr.merge(ds_item_list)\\n        ds_list.append(ds_item)\\n    ds_out = xr.concat(ds_list, \\'time\\').chunk(dict(time=1, y=-1, x=-1))\\n    return ds_out\\n\\nda = stac_to_ds(hls_items, [\\'BLUE\\', \\'GREEN\\', \\'RED\\', \\'NIR1\\', \\'SWIR1\\', \\'SWIR2\\', \\'FMASK\\'])\\nda\\n'"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "h = hls_items[0]\n",
    "\n",
    "from datetime import datetime\n",
    "import re\n",
    "def stac_to_xr(item, band):\n",
    "    if band not in all_bands:\n",
    "        sys.exit(f\"Band: {band} is not a valid input option.\" +\n",
    "                 \"Valid inputs are ALL, COASTAL-AEROSOL, BLUE, GREEN, RED, RED-EDGE1, RED-EDGE2, RED-EDGE3,\" + \n",
    "                 \" NIR1, SWIR1, SWIR2, CIRRUS, TIR1, TIR2, WATER-VAPOR, FMASK.\" + \n",
    "                 \"To request multiple layers, provide them in comma separated format with no spaces.\" +\n",
    "                 \"Unsure of the names for your bands?\" +\n",
    "                 \"--check out the README which contains a table of all bands and band names.\")\n",
    "    if item['collection'] == s30_id:\n",
    "        href = item['assets'][lut['HLSS30'][band]]['href']\n",
    "    if item['collection'] == l30_id:\n",
    "        href = item['assets'][lut['HLSL30'][band]]['href']\n",
    "    ds = xr.open_rasterio(re.sub('https://lpdaac.earthdata.nasa.gov/', '/vsis3/', href), chunks={'band':1,\n",
    "                                            'x':-1,\n",
    "                                            'y':-1}).to_dataset(name=band)\n",
    "    ds = ds.rename({'band': 'time'})\n",
    "    ds['time'] = [item['properties']['datetime']]\n",
    "    ds['time'] = ds['time'].astype('datetime64').dt.floor('D')\n",
    "    #ds = ds.assign_coords(t = datetime.strptime(item['properties']['datetime'], '%Y-%m-%dT%H:%M:%S.%fZ'))\n",
    "    return ds\n",
    "\n",
    "def stac_to_ds(items, bands):\n",
    "    ds_list = []\n",
    "    for item in items:\n",
    "        ds_item_list = []\n",
    "        for band in bands:\n",
    "            ds_item_list.append(stac_to_xr(item, band))\n",
    "            ds_item = xr.merge(ds_item_list)\n",
    "        ds_list.append(ds_item)\n",
    "    ds_out = xr.concat(ds_list, 'time').chunk(dict(time=1, y=-1, x=-1))\n",
    "    return ds_out\n",
    "\n",
    "da = stac_to_ds(hls_items, ['BLUE', 'GREEN', 'RED', 'NIR1', 'SWIR1', 'SWIR2', 'FMASK'])\n",
    "da\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from src.hls_funcs.masks import mask_hls\n",
    "da_mask = mask_hls(da['FMASK'])\n",
    "da_mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyproj import Proj\n",
    "utmProj = Proj(\"+proj=utm +zone=13U, +north +ellps=WGS84 +datum=WGS84 +units=m +no_defs\")\n",
    "bbox_utm = utmProj([bbox_num[i] for i in [0, 2]], [bbox_num[i] for i in [3, 1]]) \n",
    "tuple(bbox_utm[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "da_sub = da.loc[dict(x=slice(*tuple(bbox_utm[0])), y=slice(*tuple(bbox_utm[1])))].where(da_mask == 0)\n",
    "da_sub"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "da_stacked = da_sub.stack(z=('y', 'x')).chunk(dict(time=1, z=-1))\n",
    "da_stacked"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from src.hls_funcs.predict import pred_bm\n",
    "import pickle\n",
    "bm_mod = pickle.load(open('src/models/CPER_HLS_to_VOR_biomass_model_lr_simp.pk', 'rb'))\n",
    "da_bm = pred_bm(da_stacked, bm_mod, dim='z')\n",
    "da_bm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "da_bm = da_bm.unstack('z').persist()\n",
    "da_bm = da_bm.where(da_bm > 0)\n",
    "da_bm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from src.hls_funcs.predict import pred_cov\n",
    "ends_dict = {\n",
    "    'SD': {\n",
    "        'ndvi': 0.30,\n",
    "        'dfi': 16,\n",
    "        'bai_126': 155},\n",
    "    'GREEN': {\n",
    "        'ndvi': 0.55,\n",
    "        'dfi': 10,\n",
    "        'bai_126': 160},\n",
    "    'BARE': {\n",
    "        'ndvi': 0.10,\n",
    "        'dfi': 8,\n",
    "        'bai_126': 140}}\n",
    "da_cov = pred_cov(da_stacked, ends_dict, dim='z')\n",
    "da_cov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "da_cov = da_cov.to_array(dim='type')\n",
    "da_cov = da_cov.where((da_cov < 1.0) | (da_cov.isnull()), 1.0)\n",
    "da_cov = da_cov.where(~(da_cov.any(dim='time').isnull()))\n",
    "da_cov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import param\n",
    "import panel as pn\n",
    "import datetime as dt\n",
    "import cartopy.crs as ccrs\n",
    "import holoviews as hv\n",
    "import numpy as np\n",
    "hv.extension('bokeh', logo=False)\n",
    "pn.extension()\n",
    "def load_maps(date):\n",
    "    bm_map = da_bm.sel(time=date).hvplot(x='x',y='y', tiles='EsriImagery', crs=ccrs.UTM(13),\n",
    "                         cmap='inferno', clim=(100, 1000), colorbar=False).opts(height=500, width=750)\n",
    "    cv_map = da_cov.sel(time=date).hvplot.rgb(x='x',y='y', \n",
    "                                            bands='type', tiles='EsriImagery', crs=ccrs.UTM(13)).opts(height=500, width=750)\n",
    "    all_maps = pn.Tabs(('Biomass', bm_map), ('Cover', cv_map))\n",
    "    return all_maps\n",
    "    \n",
    "#date_slider = pn.widgets.IntSlider(name='Date Slider',\n",
    "#                                    start=0, end=len(da_bm.time), value=0)\n",
    "\n",
    "date_picker = pn.widgets.DatePicker(name='Calendar',\n",
    "                                    value=datetime.utcfromtimestamp(da_bm.time[-1].data.astype('int') * 1e-9).date(),\n",
    "                                   enabled_dates = [datetime.utcfromtimestamp(x).date() for\n",
    "                                                    x in da_bm.time.data.astype('int') * 1e-9])\n",
    "\n",
    "@pn.depends(date=date_picker.param.value)\n",
    "def load_maps_cb(date):\n",
    "    return load_maps(np.datetime64(date))\n",
    "\n",
    "pn.Row(pn.WidgetBox('Click on the date below to see avialable dates', date_picker, width=150), \n",
    "       load_maps_cb).servable()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"import param\n",
    "import panel as pn\n",
    "import datetime as dt\n",
    "import cartopy.crs as ccrs\n",
    "import holoviews as hv\n",
    "import numpy as np\n",
    "def load_map(date):\n",
    "    return da_bm.sel(time=date).hvplot(x='x',y='y',rasterize=True,tiles='EsriImagery', crs=ccrs.UTM(13),\n",
    "                         cmap='inferno', clim=(100, 1000))\n",
    "\n",
    "#date_slider = pn.widgets.IntSlider(name='Date Slider',\n",
    "#                                    start=0, end=len(da_bm.time), value=0)\n",
    "\n",
    "date_picker = pn.widgets.DatePicker(name='Date Picker',\n",
    "                                    value=datetime.utcfromtimestamp(da_bm.time[0].data.astype('int') * 1e-9).date(),\n",
    "                                   enabled_dates = [datetime.utcfromtimestamp(x).date() for x in da_bm.time.data.astype('int') * 1e-9])\n",
    "\n",
    "@pn.depends(date=date_picker.param.value)\n",
    "def load_map_cb(date):\n",
    "    return load_map(np.datetime64(date))\n",
    "\n",
    "pn.panel(pn.Row(pn.WidgetBox('Select date', date_picker), load_map_cb)).servable()\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"import param\n",
    "import panel as pn\n",
    "import datetime as dt\n",
    "import cartopy.crs as ccrs\n",
    "import holoviews as hv\n",
    "def load_map(date): \n",
    "    return da_cov.isel(time=date).hvplot.rgb(x='x',y='y', bands='type', tiles='EsriImagery', crs=ccrs.UTM(13))\n",
    "\n",
    "date_slider = pn.widgets.IntSlider(name='Date Slider',\n",
    "                                    start=0, end=len(da_cov.time), value=0)\n",
    "\n",
    "@pn.depends(date=date_slider.param.value)\n",
    "def load_map_cb(date):\n",
    "    return load_map(date)\n",
    "\n",
    "#pn.Row(pn.WidgetBox('Select date', date_slider), load_map_cb)\"\"\""
   ]
  }
 ],
 "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
}