{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3284bd7f",
   "metadata": {},
   "source": [
    "# Santa Barbara County Vineyards Analysis 2023  \n",
    "\n",
    "> This notebook is part of a numbered series (01-05) that focuses specifically on combining EMIT VSWIR and ECOSTRESS TIR data for integrated land surface analysis.\n",
    "\n",
    "**Authors:**  \n",
    "Gregory Halverson and Claire Villanueva-Weeks  \n",
    "Jet Propulsion Laboratory, California Institute of Technology\n",
    "\n",
    "This research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, and was sponsored by ECOSTRESS and the National Aeronautics and Space Administration (80NM0018D0004). \n",
    "\n",
    "© 2023. All rights reserved.\n",
    "\n",
    "### Summary and Background\n",
    " Santa Barbara County is a well-known for its viticulture. Positioned near the ocean, the unique topography of Santa Barbara allows for unique microclimates that yield diverse crops and wines. In this notebook, we aim to visualize and analyze Vineyard health using Earth observation data from two different satellite products: the EMIT L2A Reflectance product stored in NETCDF4 format and an ECOSTRESS L2 Land Surface Temperature product stored in GEOTIFF format. Our goal is to visualize Normalized Difference Vegetation Index (NDVI), Equivalent Water Thickness (EWT), Evapotranspiration (ET), and Surface Temperature (ST) interactively in `holoviews` maps and generate linear regression plots of point source data.\n",
    "\n",
    "<figure align=\"center\">\n",
    "  <img src=\"../../img/vineyard_stock_image.jpg\" width=\"200\" height=\"200\"/>\n",
    "    <figcaption>Figure 1. Vineyard.</figcaption> \n",
    "</figure>\n",
    "\n",
    "### Requirements  \n",
    " - [NASA Earthdata Account](https://urs.earthdata.nasa.gov/home)  \n",
    " - *No Python setup requirements if connected to the workshop cloud instance!*  \n",
    " - **Local Only** Set up Python Environment - See **setup_instructions.md** in the `/setup/` folder to set up a local compatible Python environment  \n",
    " - Downloaded necessary files. This is done at the end of the [01_Finding_Concurrent_Data](01_Finding_Concurrent_Data.ipynb) notebook.\n",
    "\n",
    "### Tutorial references\n",
    "01_Finding_Concurrent_Data.ipynb \n",
    "- Finding and downloading concurrent EMIT and ECOSTRESS data with earthaccess\n",
    "\n",
    "03_EMIT_CWC_from_Reflectance.ipynb\n",
    "- Calculating EWT of ROI\n",
    "\n",
    "03_EMIT_Inreractive_Points.ipynb\n",
    "- Generating a map with clickable points\n",
    "\n",
    "### Required datasets  \n",
    "- EMIT Reflectance\n",
    "    - EMIT_L2A_RFL_001_20230401T203803_2309114_003.nc\n",
    "\n",
    "- EMIT Canopy Water Content\n",
    "    - EMIT_L2A_RFL_001_20230401T203751_2309114_002_roi_bbox_cwc_merged.tif\n",
    "\n",
    "- ECOSTRESS Land Surface Temperature\n",
    "    - ECOv002_L2T_LSTE_26860_001_10SGD_20230401T203733_0710_01_LST.tif\n",
    "\n",
    "- ECOSTRESS Evapotranspiration\n",
    "    - ECOv002_L3T_JET_26860_001_10SGD_20230401T203732_0700_01_ETdaily.tif\n",
    "    \n",
    "- SB_ROI_ag.geojson: This file includes polygons of delineated agricultural field boundaries. It has been clipped to our ROI and fields were edited for ease of usability. \n",
    "\n",
    "#### Tutorial Outline  \n",
    "5.1 Set up  \n",
    "5.2 Locating Data Sources and Loading Data  \n",
    "5.3 Color Hex-Codes for Raster Plots  \n",
    "5.4 Subsetting the data over our ROI  \n",
    "5.5 Interactive Maps  \n",
    "5.6 Scatter Plots  \n",
    "\n",
    "#### References\n",
    "1. Jackson, M. (2018). Santa Barbara County (Images of America). Arcadia Publishing.\n",
    "2. [CalAgData](https://cosantabarbara.app.box.com/s/jdt95fy7gst3g8649l9t3ukrorr5xeh9)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a47dcebb",
   "metadata": {},
   "source": [
    "## 5.1 Set up \n",
    "\n",
    "These are some built-in Python functions we need for this notebook, including functions for handling filenames and dates."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f311dff9",
   "metadata": {},
   "source": [
    "We're using the `rioxarray` package for loading raster data from a GeoTIFF file, and we're importing it as `rxr`. We're using the `numpy` library to handle arrays, and we're importing it as `np`. We're using the `rasterio` package to subset the data.\n",
    "\n",
    "We're using the `geopandas` library to load vector data from GeoJSON files, and we're importing it as `gpd`. We're using the `shapely` library to handle vector data and the `pyproj` library to handle projections.\n",
    "\n",
    "Import the `emit_tools` module and call use from emit_tools import emit_xarray\n",
    "help(emit_xarray) the help function to see how it can be used.\n",
    "\n",
    "> Note: This function currently works with L1B Radiance and L2A Reflectance Data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e35d8d2f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import os\n",
    "from os.path import join, expanduser, splitext, basename\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "from osgeo import gdal\n",
    "import rasterio as rio\n",
    "import rioxarray as rxr\n",
    "import holoviews as hv\n",
    "import hvplot.xarray\n",
    "import hvplot.pandas\n",
    "import geoviews as gv\n",
    "import earthaccess\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import LinearSegmentedColormap\n",
    "\n",
    "from holoviews.streams import SingleTap\n",
    "from scipy.stats import linregress\n",
    "\n",
    "from vitals.emit_tools import emit_xarray, ortho_xr"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "049bde85",
   "metadata": {},
   "source": [
    "We can define these constants to prescribe the dimensions of our figures. Feel free to adjust these to fit your display. We're setting the alpha to make the raster semi-transparent on top of the basemap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8034a909",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "FIG_WIDTH_PX = 1080\n",
    "FIG_HEIGHT_PX = 720\n",
    "FIG_WIDTH_IN = 16\n",
    "FIG_HEIGHT_IN = 9\n",
    "FIG_ALPHA = 0.7\n",
    "BASEMAP = \"EsriImagery\"\n",
    "SEABORN_STYLE = \"whitegrid\"\n",
    "FILL_COLOR = \"none\"\n",
    "LINE_COLOR = \"red\"\n",
    "LINE_WIDTH = 3\n",
    "WEBMAP_PROJECTION = \"EPSG:4326\"\n",
    "sns.set_style(SEABORN_STYLE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62394514",
   "metadata": {},
   "source": [
    "## 5.2 Locating Data Sources and Loading Data\n",
    "\n",
    "Here we are defining location of EMIT and ECOSTRESS raster files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30187057",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "data_dir = \"../../data/\" \n",
    "data_dir"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d37057cb",
   "metadata": {},
   "source": [
    "Here we are loading in our vector data source. We can image the delineated agricultural field boundaries and take a look at the fields."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20266fcc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ag_file = join(data_dir,\"SB_ROI_ag.geojson\")\n",
    "# cut out the agroi section \n",
    "ag_latlon = gpd.read_file(ag_file)\n",
    "\n",
    "ag_fig = ag_latlon.to_crs(WEBMAP_PROJECTION).hvplot.polygons(\n",
    "    tiles=BASEMAP,\n",
    "    line_color = LINE_COLOR,\n",
    "    line_width = LINE_WIDTH,\n",
    "    fill_color = FILL_COLOR,\n",
    "    crs='EPSG:4326'\n",
    ") \n",
    "\n",
    "ag_fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c431928f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ag_latlon"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca6ba107",
   "metadata": {},
   "source": [
    "### Loading in EMIT reflectance data\n",
    "\n",
    "This notebook requires downloading an EMIT scene if this was not done in a previous notebook.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ec9e8f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get requests https Session using Earthdata Login Info\n",
    "url = 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/EMITL2ARFL.001/EMIT_L2A_RFL_001_20230401T203803_2309114_003/EMIT_L2A_RFL_001_20230401T203803_2309114_003.nc'\n",
    "earthaccess.login()\n",
    "fs = earthaccess.get_requests_https_session()\n",
    "# Retrieve granule asset ID from URL (to maintain existing naming convention)\n",
    "granule_asset_id = url.split('/')[-1]\n",
    "# Define Local Filepath\n",
    "fp = join(data_dir, granule_asset_id)\n",
    "# Download the Granule Asset if it doesn't exist\n",
    "if not os.path.isfile(fp):\n",
    "    with fs.get(url,stream=True) as src:\n",
    "        with open(fp,'wb') as dst:\n",
    "            for chunk in src.iter_content(chunk_size=64*1024*1024):\n",
    "                dst.write(chunk)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22f5026e",
   "metadata": {},
   "source": [
    "To open up the EMIT `.nc` file we will use the `netCDF4`, `xarray` and `emit_tools` libraries. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a73fbe21",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "EMIT_fp = join(data_dir, \"EMIT_L2A_RFL_001_20230401T203803_2309114_003.nc\") \n",
    "EMIT_fp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64a7e833",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "refl_ds = emit_xarray(EMIT_fp, ortho=True)\n",
    "refl_ds"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ddae16f",
   "metadata": {},
   "source": [
    "Here some functions are set up to "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed146f34",
   "metadata": {},
   "source": [
    "## 5.3 Color Hex-Codes for Raster Plots\n",
    "Here we have set up convenience functions for interpolating Color Hex-Codes for Raster Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25068015",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def interpolate_hex(hex1, hex2, ratio):\n",
    "    rgb1 = [int(hex1[i:i+2], 16) for i in (1, 3, 5)]\n",
    "    rgb2 = [int(hex2[i:i+2], 16) for i in (1, 3, 5)]\n",
    "    rgb = [int(rgb1[i] + (rgb2[i] - rgb1[i]) * ratio) for i in range(3)]\n",
    "    \n",
    "    return '#{:02x}{:02x}{:02x}'.format(*rgb)\n",
    "\n",
    "def create_gradient(colors, steps):\n",
    "    gradient = []\n",
    "    \n",
    "    for i in range(len(colors) - 1):\n",
    "        for j in range(steps):\n",
    "            ratio = j / float(steps)\n",
    "            gradient.append(interpolate_hex(colors[i], colors[i+1], ratio))\n",
    "    \n",
    "    gradient.append(colors[-1])\n",
    "    \n",
    "    return gradient\n",
    "\n",
    "def plot_cmap(cmap):\n",
    "    gradient = np.linspace(0, 1, 256)  # Gradient from 0 to 1\n",
    "    gradient = np.vstack((gradient, gradient))  # Make 2D image\n",
    "\n",
    "    # Display the colormap\n",
    "    plt.figure(figsize=(6, 2))\n",
    "    plt.imshow(gradient, aspect='auto', cmap=cmap)\n",
    "    plt.axis('off')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be1fd800",
   "metadata": {},
   "source": [
    "### Near-Infrared and Red Bands in EMIT Hyperspectral Cube"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77a64108",
   "metadata": {},
   "source": [
    "These are the `nearest` to 800 nm and 675 nm wavelengths, they can be used to calculate the NDVI using a ratio of the difference between the wavelengths to the sum of the wavelengths. NDVI is a metric by which we can estimate vegetation greenness.\n",
    "\n",
    "The `hvplot` package extends `xarray` to allow us to plot maps. We're reprojecting the raster geographic projection **EPSG 4326** to overlay on the basemap with a latitude and longitude graticule. We will be using hvplot a few more times to look at the data we are using."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0946eb26",
   "metadata": {},
   "source": [
    "#### Defining Color-Map for Near Infrared"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c78a83c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "NIR_colors = [\"#000000\", \"#FF0000\"]\n",
    "NIR_gradient = create_gradient(NIR_colors, 100)\n",
    "NIR_cmap = LinearSegmentedColormap.from_list(name=\"NIR\", colors=NIR_colors)\n",
    "plot_cmap(NIR_cmap)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b52747cc",
   "metadata": {},
   "source": [
    "#### Extracting 800 nm Near Infrared from EMIT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fef056a3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "NIR = refl_ds.sel(wavelengths=800, method=\"nearest\")\n",
    "NIR.rio.to_raster(join(data_dir,\"NIR_800.tif\"))\n",
    "NIR = rxr.open_rasterio(join(data_dir,\"NIR_800.tif\"), mask_and_scale=True).squeeze(\"band\", drop=True)\n",
    "NIR.data[NIR.data==-9999] = np.nan\n",
    "NIR_vmin, NIR_vmax = np.nanquantile(np.array(NIR), [0.02, 0.98])\n",
    "NIR_title = f\"{splitext(basename(EMIT_fp))[0]} ~800 nm\"\n",
    "\n",
    "NIR.rio.reproject(WEBMAP_PROJECTION).hvplot.image(\n",
    "    crs=WEBMAP_PROJECTION, \n",
    "    cmap=NIR_gradient, \n",
    "    alpha=FIG_ALPHA, \n",
    "    tiles=BASEMAP,\n",
    "    clim=(NIR_vmin, NIR_vmax),\n",
    "    title=NIR_title\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8743f891",
   "metadata": {},
   "source": [
    "### Defining Color-Map for Red"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c68c3954",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "red_colors = [\"#000000\", \"#FF0000\"]\n",
    "red_gradient = create_gradient(red_colors, 100)\n",
    "red_cmap = LinearSegmentedColormap.from_list(name=\"red\", colors=red_colors)\n",
    "plot_cmap(red_cmap)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7fa7d021",
   "metadata": {},
   "source": [
    "#### Extracting 675 nm Red from EMIT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a7a144b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "red = refl_ds.sel(wavelengths=675, method=\"nearest\")\n",
    "red.rio.to_raster(join(data_dir,\"red_675.tif\"))\n",
    "red = rxr.open_rasterio(join(data_dir,\"red_675.tif\")).squeeze(\"band\", drop=True)\n",
    "red.data[red.data==-9999] = np.nan\n",
    "red_vmin, red_vmax = np.nanquantile(np.array(red), [0.02, 0.98])\n",
    "red_title = f\"{splitext(basename(EMIT_fp))[0]} ~675 nm\"\n",
    "\n",
    "red.rio.reproject(WEBMAP_PROJECTION).hvplot.image(\n",
    "    crs=WEBMAP_PROJECTION, \n",
    "    cmap=red_gradient, \n",
    "    alpha=FIG_ALPHA, \n",
    "    tiles=BASEMAP,\n",
    "    clim=(red_vmin, red_vmax),\n",
    "    title=red_title\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5611f3cd",
   "metadata": {},
   "source": [
    "### Defining Color-Map for Vegetation Index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "025741aa",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "NDVI_colors=[\n",
    "    \"#0000ff\",\n",
    "    \"#000000\",\n",
    "    \"#745d1a\",\n",
    "    \"#e1dea2\",\n",
    "    \"#45ff01\",\n",
    "    \"#325e32\"\n",
    "]\n",
    "\n",
    "NDVI_gradient = create_gradient(NDVI_colors, 100)\n",
    "NDVI_cmap = LinearSegmentedColormap.from_list(name=\"NDVI\", colors=NDVI_colors)\n",
    "plot_cmap(NDVI_cmap)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e4ba0de",
   "metadata": {},
   "source": [
    "#### Calculating Vegetation Index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18976464",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "NDVI = (NIR - red)/(NIR + red)\n",
    "NDVI.rio.to_raster(join(data_dir,\"NDVI.tif\"))\n",
    "NDVI = rxr.open_rasterio(join(data_dir,\"NDVI.tif\")).squeeze(\"band\", drop=True)\n",
    "NDVI_vmin, NDVI_vmax = np.nanquantile(np.array(NDVI), [0.02, 0.98])\n",
    "NDVI_title = \"Normalized Difference Vegetation Index\"\n",
    "\n",
    "NDVI_scene_map = NDVI.rio.reproject(WEBMAP_PROJECTION).hvplot.image(\n",
    "    crs=WEBMAP_PROJECTION, \n",
    "    cmap=NDVI_gradient, \n",
    "    alpha=FIG_ALPHA, \n",
    "    tiles=BASEMAP,\n",
    "    clim=(NDVI_vmin, NDVI_vmax),\n",
    "    title=NDVI_title\n",
    ")\n",
    "\n",
    "NDVI_scene_map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a0f03de",
   "metadata": {},
   "source": [
    "### Loading in and visualizing our EWT data \n",
    "\n",
    "We can now open up the EWT data that you calculated earlier with the LP DAAC EWT. The dataset is now in .tif format which means we can just open this one up"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77e4bbe6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "EWT_filename = join(data_dir,'EMIT_L2A_RFL_001_20230401T203751_2309114_002_roi_bbox_cwc_merged.tif')\n",
    "EWT_filename"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9089514f",
   "metadata": {},
   "source": [
    "#### Defining Color-Map for EWT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6df691a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "EWT_colors = [\"#FFFFFF\", \"#0000FF\"]\n",
    "EWT_gradient = create_gradient(EWT_colors, 100)\n",
    "EWT_cmap = LinearSegmentedColormap.from_list(name=\"EWT\", colors=EWT_colors)\n",
    "plot_cmap(EWT_cmap)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7fec7b9b",
   "metadata": {},
   "source": [
    "#### Mapping EWT with color ramp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3db5968",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "EWT = rxr.open_rasterio(EWT_filename).squeeze(\"band\", drop=True)\n",
    "EWT_vmin, EWT_vmax = np.nanquantile(np.array(EWT), [0.02, 0.98])\n",
    "EWT_title = \"Equivalent Water Thickness \"\n",
    "\n",
    "EWT.rio.reproject(WEBMAP_PROJECTION,nodata=np.nan).hvplot.image(\n",
    "    crs=WEBMAP_PROJECTION, \n",
    "    cmap=EWT_gradient, \n",
    "    alpha=FIG_ALPHA, \n",
    "    tiles=BASEMAP,\n",
    "    clim=(EWT_vmin, EWT_vmax),\n",
    "    title=EWT_title\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05fbfd7c",
   "metadata": {},
   "source": [
    "### Loading an ECOSTRESS LST granule"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffa3237c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ST_filename = join(data_dir, \"ECOv002_L2T_LSTE_26860_001_10SGD_20230401T203733_0710_01_LST.tif\") \n",
    "ST_filename"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db9eadcd",
   "metadata": {},
   "source": [
    "The temperatures in the `L2T_LSTE` product are given in Kelvin. To convert them to Celsius, we subtract 273.15."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08b679f0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ST_K = rxr.open_rasterio(ST_filename).squeeze('band', drop=True)\n",
    "ST_C = ST_K - 273.15\n",
    "\n",
    "ST_C.rio.to_raster(join(data_dir,\"ST.tif\"))\n",
    "\n",
    "ST_vmin, ST_vmax = np.nanquantile(np.array(ST_C), [0.02, 0.98])\n",
    "ST_title = \"ECOSTRESS Surface Temperature (Celsius)\"\n",
    "\n",
    "ST_C.rio.reproject(WEBMAP_PROJECTION,nodata=np.nan).hvplot.image(\n",
    "    crs=WEBMAP_PROJECTION, \n",
    "    cmap=\"jet\", \n",
    "    alpha=FIG_ALPHA, \n",
    "    tiles=BASEMAP,\n",
    "    clim=(ST_vmin, ST_vmax),\n",
    "    title=ST_title\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7834b0c",
   "metadata": {},
   "source": [
    "### Defining Color Ramps \n",
    "\n",
    "Color ramps or maps, are used to map valued over a range of colors based on the data and goals of visualization. Here we show how to create a custom color ramp for the data using hex codes. The examples shown are red-green colorblind friendly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4faa5c43",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ST_colors = [\n",
    "    \"#054fb9\",\n",
    "    \"#0073e6\",\n",
    "    \"#8babf1\",\n",
    "    \"#cccccc\",\n",
    "    \"#e1ad01\",\n",
    "    \"#f57600\",\n",
    "    \"#c44601\"\n",
    "]\n",
    "\n",
    "ST_gradient = create_gradient(ST_colors, 100)\n",
    "ST_cmap = LinearSegmentedColormap.from_list(name=\"ST\", colors=ST_colors)\n",
    "plot_cmap(ST_cmap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1b314d4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ST_vmin, ST_vmax = np.nanquantile(np.array(ST_C), [0.02, 0.98])\n",
    "ST_title = \"ECOSTRESS Surface Temperature (Celsius)\"\n",
    "\n",
    "ST_C.rio.reproject(WEBMAP_PROJECTION).hvplot.image(\n",
    "    crs=WEBMAP_PROJECTION, \n",
    "    cmap=ST_gradient, \n",
    "    alpha=FIG_ALPHA, \n",
    "    tiles=BASEMAP,\n",
    "    clim=(ST_vmin, ST_vmax),\n",
    "    title=ST_title\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f393528",
   "metadata": {},
   "source": [
    "### Loading ECOSTRESS Evapotranspiration granule"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ad78d0c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ET_filename = join(data_dir, \"ECOv002_L3T_JET_26860_001_10SGD_20230401T203732_0700_01_ETdaily.tif\")\n",
    "ET_filename"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6f737ce",
   "metadata": {},
   "source": [
    "Now create a color ramp representing ECOSTRESS evapotranspiration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6338b083",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ET_colors = [\n",
    "    \"#f6e8c3\",\n",
    "    \"#d8b365\",\n",
    "    \"#99974a\",\n",
    "    \"#53792d\",\n",
    "    \"#6bdfd2\",\n",
    "    \"#1839c5\"\n",
    "]\n",
    "\n",
    "ET_gradient = create_gradient(ET_colors, 100)\n",
    "ET_cmap = LinearSegmentedColormap.from_list(name=\"ET\", colors=ET_colors)\n",
    "plot_cmap(ET_cmap)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7286942c",
   "metadata": {},
   "source": [
    "The daily evapotranspiration product from ECOSTRESS is given in millimeters per day."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf3d73dc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ET = rxr.open_rasterio(ET_filename).squeeze('band', drop=True)\n",
    "ET_vmin, ET_vmax = np.nanquantile(np.array(ET), [0.02, 0.98])\n",
    "ET_title = \"ECOSTRESS Evapotranspiration (mm / day)\"\n",
    "\n",
    "ET.rio.reproject(WEBMAP_PROJECTION).hvplot.image(\n",
    "    crs=WEBMAP_PROJECTION, \n",
    "    cmap=ET_gradient, \n",
    "    alpha=FIG_ALPHA, \n",
    "    tiles=BASEMAP,\n",
    "    clim=(ET_vmin, ET_vmax),\n",
    "    title=ET_title\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a767b923",
   "metadata": {},
   "source": [
    "## 5.4. Subsetting the data over our ROI\n",
    "\n",
    "To clip the raster image to the extent of the vector dataset, we want to subset the raster to the bounds of the vector dataset. This dataset is included here in GeoJSON format, which we'll load in as a geodataframe using the `geopandas` package.\n",
    "\n",
    "The GeoJSON is provided by CalAg and includes information on agricultural field boundaries in the Santa Barbara County. We can use this code to select the largest polygons that are classified as vineyards.\n",
    "\n",
    "Using the rioxarray `clip` function, we can subset the data using a geojson. to visualize the area that are vineyards are located, I have selected a boundary polygon. Overlaid are the outlines of some of the larger vineyards in the scene.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39cbf00c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "target_text = 'WINE'\n",
    "\n",
    "filtered_gdf = ag_latlon[ag_latlon['crop_list'].str.contains(target_text, case=False, na=False)]\n",
    "\n",
    "Q1 = filtered_gdf['size'].quantile(0.25)\n",
    "Q3 = filtered_gdf['size'].quantile(0.75)\n",
    "IQR = Q3 - Q1\n",
    "\n",
    "# Define the threshold as Q3 + 1.5 * IQR (adjust multiplier as needed)\n",
    "threshold_size = Q3 + 1.5 * IQR\n",
    "\n",
    "# Filter polygons based on size greater than the threshold\n",
    "vineyard_polygons = filtered_gdf[filtered_gdf['size'] > threshold_size]\n",
    "vineyard_polygons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41165ce9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "vineyard_plot = vineyard_polygons.to_crs(\"EPSG:4326\").hvplot.polygons(\n",
    "    crs=WEBMAP_PROJECTION,\n",
    "    tiles=BASEMAP,\n",
    "    alpha=FIG_ALPHA,\n",
    "    width=FIG_WIDTH_PX,\n",
    "    height=FIG_HEIGHT_PX,\n",
    "    fill_color = 'none',\n",
    "    line_color = 'r',\n",
    "    line_width = 3,\n",
    "    title=\"Vineyards of Acreage greater than IQR3\",\n",
    "    fontscale=1.5\n",
    "    )\n",
    "\n",
    "hv.extension('bokeh')\n",
    "vineyard_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82d310e8",
   "metadata": {},
   "source": [
    "## Interactively mapping our data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11f96d15",
   "metadata": {},
   "source": [
    "Now that we have selected a subset of the agriculture in the area we can subset our data using the `rioxarray` function `clip`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9b9436d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "NDVI_vineyard = NDVI.rio.clip(vineyard_polygons.geometry.values,vineyard_polygons.crs, all_touched=True)\n",
    "NDVI_vineyard.rio.to_raster(join(data_dir,\"NDVI_vineyard.tif\"))\n",
    "NDVI_vineyard = rxr.open_rasterio(join(data_dir,\"NDVI_vineyard.tif\")).squeeze(\"band\", drop=True)\n",
    "NDVI_vmin, NDVI_vmax = np.nanquantile(np.array(NDVI_vineyard), [0.02, 0.98])\n",
    "\n",
    "ST_vineyard = ST_C.rio.clip(vineyard_polygons.geometry.values,vineyard_polygons.crs, all_touched=True)\n",
    "ST_vineyard.rio.to_raster(join(data_dir,\"ST_vineyard.tif\"))\n",
    "ST_vineyard = rxr.open_rasterio(join(data_dir,\"ST_vineyard.tif\")).squeeze(\"band\", drop=True)\n",
    "ST_vmin, ST_vmax = np.nanquantile(np.array(ST_vineyard), [0.02, 0.98])\n",
    "\n",
    "EWT_vineyard = EWT.rio.clip(vineyard_polygons.geometry.values,vineyard_polygons.crs, all_touched=True)\n",
    "EWT_vineyard.rio.to_raster(join(data_dir,\"EWT_vineyard.tif\"))\n",
    "EWT_vineyard = rxr.open_rasterio(join(data_dir,\"EWT_vineyard.tif\")).squeeze(\"band\", drop=True)\n",
    "EWT_vmin, EWT_vmax = np.nanquantile(np.array(EWT_vineyard), [0.02, 0.98])\n",
    "\n",
    "ET_vineyard = ET.rio.clip(vineyard_polygons.geometry.values,vineyard_polygons.crs, all_touched=True)\n",
    "ET_vineyard.rio.to_raster(join(data_dir,\"ET_vineyard.tif\"))\n",
    "ET_vineyard = rxr.open_rasterio(join(data_dir,\"ET_vineyard.tif\")).squeeze(\"band\", drop=True)\n",
    "ET_vmin, ET_vmax = np.nanquantile(np.array(ET_vineyard), [0.02, 0.98])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d4e68a0",
   "metadata": {},
   "source": [
    "### 5.5 Interactive Maps"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5e68d13",
   "metadata": {},
   "source": [
    "We can use `holoviews` to plot an interactive clickable map of these subsets over our ROI. When you click a point on the map, you can see the coordinates and data for that pixel. **You can uncomment the crs and tiles option lines below to visualize the data on a basemap, but this breaks the interactive feature. Additionally, the interactive feature seems to also affect the color bar min and max limits (clim) whether applied manually or automatically. This may cause the color ramp to display incorrectly. Selecting and exporting values will work properly.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11d4a868",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "clicked_values = []\n",
    "\n",
    "# Reproject raster data\n",
    "NDVI_vineyard_projected = NDVI_vineyard.rio.reproject(WEBMAP_PROJECTION, nodata=np.nan)\n",
    "ST_vineyard_projected = ST_vineyard.rio.reproject(WEBMAP_PROJECTION, nodata=np.nan)\n",
    "EWT_vineyard_projected = EWT_vineyard.rio.reproject(WEBMAP_PROJECTION, nodata=np.nan)\n",
    "ET_vineyard_projected = ET_vineyard.rio.reproject(WEBMAP_PROJECTION, nodata=np.nan)\n",
    "\n",
    "# Load geospatial data\n",
    "vineyard_polygons = gpd.read_file(join(data_dir,\"vineyard_polygons_filtered.geojson\"))\n",
    "\n",
    "# Set up the SingleTap stream\n",
    "stream = SingleTap()\n",
    "\n",
    "# Define a function to process clicks\n",
    "def interactive_click(x, y):\n",
    "    if None not in [x, y]:\n",
    "        NDVI_value = NDVI_vineyard_projected.sel(x=x, y=y, method=\"nearest\").values\n",
    "        ST_value = ST_vineyard_projected.sel(x=x, y=y, method=\"nearest\").values\n",
    "        ET_value = ET_vineyard_projected.sel(x=x, y=y, method=\"nearest\").values\n",
    "        EWT_value = EWT_vineyard_projected.sel(x=x, y=y, method=\"nearest\").values\n",
    "        print(f\"Lon: {x}, Lat: {y}, NDVI: {NDVI_value}, ST: {ST_value}, EWT: {EWT_value}, ET: {ET_value}\")\n",
    "        clicked_values.append({ \"Lon\": x, \"Lat\": y, \"NDVI\": NDVI_value, \"ST\": ST_value, \"ET\": ET_value, \"EWT\": EWT_value})\n",
    "\n",
    "    return hv.Points((x, y)).opts(color='green', size=10)\n",
    "\n",
    "# Create two separate plots and share the stream between them\n",
    "NDVI_vineyard_map = NDVI_vineyard_projected.hvplot.image(\n",
    "    #crs=WEBMAP_PROJECTION,\n",
    "    cmap=NDVI_gradient,\n",
    "    clim=(NDVI_vmin, NDVI_vmax),\n",
    "    alpha=FIG_ALPHA, \n",
    "    #tiles=BASEMAP,\n",
    "    title=NDVI_title\n",
    ") * vineyard_polygons.to_crs(WEBMAP_PROJECTION).hvplot.polygons(\n",
    "    #crs=WEBMAP_PROJECTION,\n",
    "    line_color=LINE_COLOR,\n",
    "    line_width=LINE_WIDTH,\n",
    "    fill_color=FILL_COLOR\n",
    ").opts(width=600, height=400, fontscale=1.5)\n",
    "\n",
    "ST_vineyard_map = ST_vineyard_projected.hvplot.image(\n",
    "    #crs=WEBMAP_PROJECTION,\n",
    "    cmap=ST_gradient, \n",
    "    #clim=(ST_vmin, ST_vmax),\n",
    "    alpha=FIG_ALPHA, \n",
    "    #tiles=BASEMAP,\n",
    "    title=ST_title\n",
    ") * vineyard_polygons.to_crs(WEBMAP_PROJECTION).hvplot.polygons(\n",
    "    #crs=WEBMAP_PROJECTION,\n",
    "    line_color=LINE_COLOR,\n",
    "    line_width=LINE_WIDTH,\n",
    "    fill_color=FILL_COLOR\n",
    ").opts(width=600, height=400, fontscale=1.5)\n",
    "\n",
    "# Duplicate the plots for the second row\n",
    "EWT_vineyard_map = EWT_vineyard_projected.hvplot.image(\n",
    "    #crs=WEBMAP_PROJECTION,\n",
    "    cmap=EWT_gradient,\n",
    "    #clim=(EWT_vmin, EWT_vmax), \n",
    "    alpha=FIG_ALPHA, \n",
    "    #tiles=BASEMAP,\n",
    "    title=EWT_title\n",
    ") * vineyard_polygons.to_crs(WEBMAP_PROJECTION).hvplot.polygons(\n",
    "    #crs=WEBMAP_PROJECTION,\n",
    "    line_color=LINE_COLOR,\n",
    "    line_width=LINE_WIDTH,\n",
    "    fill_color=FILL_COLOR\n",
    ").opts(width=600, height=400, fontscale=1.5)\n",
    "\n",
    "ET_vineyard_map = ET_vineyard_projected.hvplot.image(\n",
    "    #crs=WEBMAP_PROJECTION,\n",
    "    rasterize=True, \n",
    "    cmap=ET_gradient, \n",
    "    #clim=(ET_vmin, ET_vmax),\n",
    "    alpha=FIG_ALPHA, \n",
    "    #tiles=BASEMAP,\n",
    "    title=ET_title\n",
    ") * vineyard_polygons.to_crs(WEBMAP_PROJECTION).hvplot.polygons(\n",
    "    #crs=WEBMAP_PROJECTION,\n",
    "    line_color=LINE_COLOR,\n",
    "    line_width=LINE_WIDTH,\n",
    "    fill_color=FILL_COLOR\n",
    ").opts(width=600, height=400, fontscale=1.5)\n",
    "\n",
    "# Create a DynamicMap for interactive clicks\n",
    "click = hv.DynamicMap(interactive_click, streams=[stream])\n",
    "\n",
    "# Display the plots in a 2x2 grid with two rows along with the interactive click behavior\n",
    "layout = (click * NDVI_vineyard_map * click + click * ST_vineyard_map * click +\n",
    "          click * EWT_vineyard_map * click + click * ET_vineyard_map * click).cols(2)\n",
    "layout\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21fe620b",
   "metadata": {},
   "source": [
    "Save the data from our clicked points into a `csv` that we can use for plotting.\n",
    "\n",
    "By default, saving is disabled below so this notebook can run unattended using a preselected `clicked_values.csv`. If you'd like to use your own clicked points instead, uncomment the cell below -- this will overwrite the tracked `clicked_values.csv` with whatever you've clicked above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd5320ed",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Save the clicked values to a CSV file after all clicks are processed\n",
    "# if clicked_values:\n",
    "#     df = pd.DataFrame(clicked_values)\n",
    "#     df.to_csv(join(data_dir,\"clicked_values.csv\"), index=False)\n",
    "#     print(\"Clicked values saved to 'clicked_values.csv'\")\n",
    "# else:\n",
    "#     print(\"No clicked values to save.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90eef3fb",
   "metadata": {},
   "source": [
    "### 5.6 Scatter Plots\n",
    "\n",
    "Now we can use the data saved into our csv to create scatter plots to visualize the correlations among ET, ST, EWT, and NDVI. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94749d11",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Replace \"your_data.csv\" with the actual CSV file containing your data\n",
    "df = pd.read_csv(join(data_dir,\"clicked_values.csv\"))\n",
    "\n",
    "# Filter out Lat/Lon Columns\n",
    "filtered_columns = df.columns[2:]\n",
    "\n",
    "# Get combinations of all pairs of columns\n",
    "combinations = [(col1, col2) for i, col1 in enumerate(filtered_columns) for col2 in filtered_columns[i + 1:]]\n",
    "\n",
    "# Set up subplots without equal aspect ratio and sharex/sharey set to False\n",
    "fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(10, 12), sharex=False, sharey=False)\n",
    "fig.subplots_adjust(hspace=0.5)\n",
    "\n",
    "# Iterate over combinations and create scatter plots with best-fit lines\n",
    "for (col1, col2), ax in zip(combinations, axes.flatten()):\n",
    "    sns.regplot(x=col1, y=col2, data=df, ax=ax, ci=None, line_kws={\"color\": \"red\"})\n",
    "\n",
    "    # Filter out rows with missing values in the selected columns\n",
    "    filtered_df = df[[col1, col2]].dropna()\n",
    "\n",
    "    # Calculate and plot the best-fit line using linear regression\n",
    "    slope, intercept, r_value, p_value, std_err = linregress(filtered_df[col1], filtered_df[col2])\n",
    "    x_range = np.linspace(filtered_df[col1].min(), filtered_df[col1].max(), 100)\n",
    "    y_fit = slope * x_range + intercept\n",
    "    ax.plot(x_range, y_fit, color='blue', linestyle='--', label=f'Best Fit (R={r_value:.2f})')\n",
    "\n",
    "    ax.legend()\n",
    "    ax.set_xlabel(col1)\n",
    "    ax.set_ylabel(col2)\n",
    "    ax.set_title(f\"{col1} vs {col2}\")\n",
    "\n",
    "    # Set x and y axis ticks with different values\n",
    "    ax.set_xticks(np.linspace(filtered_df[col1].min(), filtered_df[col1].max(), 5))\n",
    "    ax.set_yticks(np.linspace(filtered_df[col2].min(), filtered_df[col2].max(), 5))\n",
    "\n",
    "# Adjust layout\n",
    "plt.tight_layout()\n",
    "\n",
    "# Show the plot\n",
    "plt.show()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.13.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}