Sharpening_ECOSTRESS_S2_API.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sharpening ECOSTRESS LST using Sentinel-2 VSWIR products at High Resolution (<20m) including the downloading of ECOSTRESS granules"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notebook : Quentin Dehaene, mentored by Glynn Hulley    \n",
    "Original Python Implementation : [Radoslaw Guzinski](https://github.com/radosuav/pyDMS)  \n",
    "Original Implemenation : [Gao et al.](https://doi.org/10.3390/rs4113287) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Points of contact : quentindehaene@gmail.com and glynn.hulley@jpl.nasa.gov"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### How does the sharpening algorithm, pyDMS, work?\n",
    "\n",
    "We'll call the ECOSTRESS LST at native 70m **LST-LR** and we'll refer to the Sentinel-2 reflectance multiband raster at 20m as **S2-HR**.\n",
    "\n",
    "The first step is to **reproject and subset LST-LR to S2-HR Coordinates Reference System (CRS) and extent**. This means that the result, and all images in the process, will be in the S2-HR CRS and will share its extent, even if it implies padding a smaller LST-LR with no-data values to match the extents. Conversely, all the images larger than S2-HR would be clipped to its extent.\n",
    "\n",
    "The second step consists of **resampling S2-HR to the coarse resolution of LST-LR** (producing a temporary file). In the process, we compute the **homogeneity inside each coarse resolution pixel**.\n",
    "\n",
    "$$\n",
    "c_v = \\left( \\frac{1}{n} \\right) \\sum_{i=1}^{n} \\left( \\frac{\\sigma_{i}}{\\mu_{i}} \\right)\n",
    "$$\n",
    "\n",
    "where *i* represents the spectral band, *n* is the total number of bands, and $\\mu$ and $\\sigma$ are the mean and standard deviation of the fine resolution reflectances within a coarse resolution pixel.\n",
    "\n",
    "The $c_v$ is the **coefficient of variation** of the pixel; the closer to 0, the purer the pixel is. This will be used as **weights in our training**, with the most homogeneous pixels being weighted more and the least homogeneous being discarded (with a threshold set as $c_v = 0.2$), leaving around 80% of the pixels available for training usually.\n",
    "\n",
    "Using this **resampled Sentinel-2 reflectance (S2-LR)** and **LST-LR**, we can train the regression model.\n",
    "\n",
    "**Some more details about the model:**\n",
    "\n",
    "We are training a **bagging of an ensemble of regression trees**. Each of these trees is slightly more complicated than the usual random forest regression tree. The principle of the tree is classic: each node is built to minimize the **mean square error (MSE)** of the subsamples on each side of the node until we reach a satisfying tree depth or the minimum number of samples per leaf. The difference here lies in the way we determine the target value of each leaf: instead of assigning the average of the *y* values (i.e., the LST) to all features in the leaf, we apply a **Bayesian linear regression**. This means that for each feature in the same node, the target values obey the same linear function of the features (the reflectances). So each feature's target value depends on the actual values of the feature.\n",
    "\n",
    "The model is trained using the resampled **S2 imagery (S2-LR)** as $X_{train}$ and the reprojected **LST-LR** as $y_{train}$ (ground truth). Our forest is trained at low resolution because it is at this resolution that we have what we know to be true and can establish a link between reflectance and LST.\n",
    "\n",
    "To predict a high-resolution LST map, we'll use the newly trained model with the **S2-HR** as $X$. We thus obtain a first prediction of the LST ($y_{pred}$) at high resolution.\n",
    "\n",
    "Then comes the **residual analysis**. We resample this prediction to the coarse resolution and compute the residual ($y - y_{pred}$) (the only truth we have is at the low resolution). This compares our predicted LST, resampled at the coarse resolution, to the original **LST-LR**. For that, we don’t actually compare the LST directly but the $LST^4$, which is proportional to the exitance. It makes more sense since, physically, temperature doesn’t have to be conserved in the sharpening, but energy does. Sensors don’t measure temperature but radiance.\n",
    "\n",
    "The final step is to **smooth the residual**, resample it to the high resolution, and sum it with $y_{pred}$, our predicted LST. We now have a final LST prediction that verifies that the average radiance of all the high-resolution pixels composing a coarse pixel is equal to the radiance in the original LST.\n",
    "\n",
    "This is the downscaled image that we were looking for.\n",
    "\n",
    "---\n",
    "\n",
    "I would like to note that there are two other sharpening techniques implemented in the pyDMS code, developed (and still being developed as of October 2024) by R. Guzinski. All the steps are the same; the only difference is the regression model itself. There is, on one hand, the **Neural Network regressor**—instead of the trees presented here, we are using MLP trees from [scikit-neuralnetwork](https://github.com/aigamedev/scikit-neuralnetwork). On the other hand, the model closer to the original proposition of Gao is based on a [Cubist](https://www.rulequest.com/cubist-unix.html) regressor, recently made available in Python.\n",
    "\n",
    "As of October 2024, I haven't used the other models enough to express confidence in the results. This will require testing.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import cell\n",
    "from osgeo import gdal\n",
    "import rasterio\n",
    "import numpy as np\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "import getpass\n",
    "import requests as req\n",
    "from datetime import datetime \n",
    "import pandas as pd\n",
    "from shapely.geometry import Polygon\n",
    "import json\n",
    "import time\n",
    "import geopandas as gpd\n",
    "from sentinelhub import (SHConfig, DataCollection, SentinelHubCatalog, SentinelHubRequest, BBox, bbox_to_dimensions, CRS, MimeType, Geometry,MosaickingOrder)\n",
    "import rioxarray as rxr\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "import matplotlib.animation as animation\n",
    "from matplotlib import rc\n",
    "from pyDMS_main.run_pyDMS import *\n",
    "from dotenv import load_dotenv\n",
    "# If you recieve a No module named '...' error, it is likely because you haven't installed all the necessary packages (cf tutorial)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up and preprocessing\n",
    "This whole notebook is actually a wrap over the original run_pyDMS. It is intended to render the sharpening easy and automatic by dealing with an entire folder at a time. Indeed, when interested in an area, we are often looking at a series of images, a heatwave week, a summer month or more. The point here is to have all the files from your [AppEEARS](https://appeears.earthdatacloud.nasa.gov/) request in the same folders : This version is easier to use than the other one available, since you don't need to use [AppEEARS](https://appeears.earthdatacloud.nasa.gov/). Indeed, this code use the API, allowing to create a request directly from python and then to download the request files automatically too. This will simply require your [Earthdata](http://urs.earthdata.nasa.gov) logs and the inputs of your choosing. All the files downloaded will be downscaled using a single Sentinel 2 VSWIR image, downloaded here too, after being preprocessed ealier in the code. To avoid some errors due to seasonal effect (most notably changes in vegetation), I would advise to not sharpen different seasons with the same S2 image.\n",
    "The output files, will all be written in one folder in the GEOTIFF format, that you can use in any GIS software. For each image a residual GEOTIFF will also be produced, but in most cases you can ignore these files.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For convenience, all of the inputs requested from the user are grouped the next 4 cells."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up the dowloads and parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All the Sentinel data is free to access, but it requires to create an account to download data.\n",
    "First, if you don't have an account on [Copernicus Data Space](https://dataspace.copernicus.eu/), create one and log in.  \n",
    "Now access your User Settings : My Account > DashBoard > User Settings (bottom left)  \n",
    "Create a new OAuth client. And save your newly given token ID and password (called a secret here).  \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OAuth Copernicus Data Space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "config = SHConfig()\n",
    "config.sh_base_url = 'https://sh.dataspace.copernicus.eu'\n",
    "config.sh_token_url = 'https://identity.dataspace.copernicus.eu/auth/realms/CDSE/protocol/openid-connect/token'\n",
    "load_dotenv()\n",
    "config.sh_client_id = os.getenv('OAUTH_CLIENT_ID') # Alternatively you can type you client id\n",
    "config.sh_client_secret = os.getenv('OAUTH_CLIENT_SECRET') # Alternatively you can type your client secret"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Type your NASA Earthdata login and password. If you don't have an account you can create one on [Earthdata](http://urs.earthdata.nasa.gov). If you have any trouble refer to the relative tutorial available [here](https://github.com/ECOSTRESS-Tutorials/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "user = getpass.getpass(prompt = 'Enter NASA Earthdata Login Username: ')      # Input NASA Earthdata Login Username\n",
    "password = getpass.getpass(prompt = 'Enter NASA Earthdata Login Password: ')  # Input NASA Earthdata Login Password"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This cell is made for you to type the directories you want the dowloaded products and results to be written into. When you type a directory, it shouldn't include the final / or \\ (otherwise you will get a syntax error)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Choose your output folder for the downloaded Sentinel-2 products\n",
    "s2_output_folder = r''\n",
    "# The task name you enter here is the name of the AppEEARS request. Your S2 request will also be named accordingly for harmonization. \n",
    "task_name = input('Enter a Task Name: ')\n",
    "# Folder where all the ECOSTRESS products will be downloaded. This folder will then contain a subdirectory named after the task name above. \n",
    "eco_output_folder = r'' # In the newly created subdirectory, you'll see a QC file folder, an LST folder, an LST_scaled folder and a cloud_mask folder if you're using Collection 2.\n",
    "# Folder where all the sharpened ECOSTRESS LST files will be written for the scene of interest\n",
    "dst_dir = r''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the parameters for the products to be downloaded"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The coordinates of the bounding box of your choosing, in lat lon : (xmin,ymin,xmax,ymax)\n",
    "# Use the http://bboxfinder.com to find your box easily (already in the right order)\n",
    "aoi_coords_wgs_84 = (23.630219,37.909534,23.992081,38.076204) # example\n",
    "\n",
    "# Choose the resolution of the Sentinel data in meters 10 or 20, this will be the final resolution of the downscaled image. I would advise 20m at all times until further progress.\n",
    "s2_res = 20\n",
    "# Set the desired collection for the ECOSTRESS product to be downloaded (1 or 2). Pay attention Collection 2 hasn't been reprocessed for the full years of service of ECOSTRESS \n",
    "eco_collec = 1 \n",
    "\n",
    "# Choose your time interval (beginning, end) in the format (YYYY-MM-DD). You'll get a S2 image from the tile with the lowest cloud coverage available in the interval and all the ECOSTRESS scences overpassing the AOI during this interval : \n",
    "interval = (\"2023-07-01\", \"2023-08-31\") # example\n",
    "\n",
    "startYear,startMonth,startDay = interval[0][:4],interval[0][5:7],interval[0][8:]\n",
    "endYear,endMonth,endDay = interval[1][:4],interval[1][5:7],interval[1][8:]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Downloading Sentinel 2 product"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "If you encounter problem with the S2 imagery download, please refer to this [Documentation](https://sentinelhub-py.readthedocs.io/en/latest/index.html) or to this [Copernicus Web page](https://dataspace.copernicus.eu/news/2023-9-28-accessing-sentinel-mission-data-new-copernicus-data-space-ecosystem-apis)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Change the coordinates into a box for the Coipernicus API.  \n",
    "The bounding box is limited to 2500 pixels, you will recieve an error if the request is larger."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aoi_bbox = BBox(bbox = aoi_coords_wgs_84,crs=CRS.WGS84)\n",
    "aoi_size = bbox_to_dimensions(aoi_bbox,resolution = s2_res)\n",
    "print(f\"Image shape at {s2_res} m resolution: {aoi_size} pixels\") # The size of the box to 2500 pixels in each direction\n",
    "if (aoi_size[0]>2499 or aoi_size[1]>2499) : \n",
    "    raise(ValueError(\"The box is limited to 2500 pixels in each direction, try again with a smaller bounding box.\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Download the S2 image with the previously defined parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Request scripts, based off the sentinelhub documentation\n",
    "# This script will be used to download all the S2 whose resolution is 20m or below\n",
    "evalscript_all_bands_u20 = \"\"\"\n",
    "    //VERSION=3\n",
    "    function setup() {\n",
    "        return {\n",
    "            input: [{\n",
    "                bands: [\"B02\",\"B03\",\"B04\",\"B05\",\"B06\",\"B07\",\"B08\",\"B8A\",\"B11\",\"B12\"],\n",
    "                units: \"REFLECTANCE\"\n",
    "\n",
    "            }],\n",
    "            output: {\n",
    "                bands: 10,\n",
    "                sampleType: \"FLOAT32\"\n",
    "            }\n",
    "        };\n",
    "    }\n",
    "\n",
    "    function evaluatePixel(sample) {\n",
    "        return [sample.B02,\n",
    "                sample.B03,\n",
    "                sample.B04,\n",
    "                sample.B05,\n",
    "                sample.B06,\n",
    "                sample.B07,\n",
    "                sample.B08,\n",
    "                sample.B8A,\n",
    "                sample.B11,\n",
    "                sample.B12];\n",
    "    }\n",
    "\"\"\"\n",
    "# This script will be used to download all the S2 whose resolution is 10m\n",
    "evalscript_all_bands_u10 = \"\"\"\n",
    "    //VERSION=3\n",
    "    function setup() {\n",
    "        return {\n",
    "            input: [{\n",
    "                bands: [\"B02\",\"B03\",\"B04\",\"B08\"],\n",
    "                units: 'REFLECTANCE'\n",
    "            }],\n",
    "            output: {\n",
    "                bands: 4,\n",
    "                sampleType: \"FLOAT32\"\n",
    "            }\n",
    "        };\n",
    "    }\n",
    "\n",
    "    function evaluatePixel(sample) {\n",
    "        return [sample.B02,\n",
    "                sample.B03,\n",
    "                sample.B04,\n",
    "                sample.B08,];\n",
    "    }\n",
    "\"\"\"\n",
    "# Request the data at 20m\n",
    "if s2_res == 20 :\n",
    "    request_all_bands_u20 = SentinelHubRequest(\n",
    "        data_folder=s2_output_folder,\n",
    "        evalscript=evalscript_all_bands_u20,\n",
    "        input_data=[\n",
    "            SentinelHubRequest.input_data(\n",
    "                data_collection = DataCollection.SENTINEL2_L2A.define_from(\"s2l2a\", service_url=config.sh_base_url),\n",
    "                time_interval = interval,\n",
    "                mosaicking_order=MosaickingOrder.LEAST_CC, # selecting the tile in the interval with the least cloud coverage\n",
    "\n",
    "            )\n",
    "        ],\n",
    "        responses=[SentinelHubRequest.output_response(\"default\", MimeType.TIFF)],\n",
    "        bbox=aoi_bbox,\n",
    "        size=aoi_size,\n",
    "        config=config,\n",
    "    )\n",
    "    resp = request_all_bands_u20.save_data(show_progress = True)\n",
    "# Request the data at 10m\n",
    "elif s2_res == 10 :\n",
    "    request_all_bands_u10 = SentinelHubRequest(\n",
    "        data_folder=s2_output_folder,\n",
    "        evalscript=evalscript_all_bands_u10,\n",
    "        input_data=[\n",
    "            SentinelHubRequest.input_data(\n",
    "                data_collection = DataCollection.SENTINEL2_L2A.define_from(\"s2l2a\", service_url=config.sh_base_url),\n",
    "            time_interval = interval,\n",
    "                mosaicking_order = MosaickingOrder.LEAST_CC\n",
    "            )\n",
    "        ],\n",
    "        responses=[SentinelHubRequest.output_response(\"default\", MimeType.TIFF)],\n",
    "        bbox=aoi_bbox,\n",
    "        size=aoi_size,\n",
    "        config=config,\n",
    "    )\n",
    "\n",
    "    resp = request_all_bands_u10.save_data(show_progress= True)\n",
    "else :\n",
    "    raise(ValueError('Unexpected resolution please change it to 10 or 20'))\n",
    "\n",
    "# Rename the downloaded files for clarity\n",
    "dirs = [d for d in os.listdir(s2_output_folder)]\n",
    "\n",
    "most_recent_dir = max(dirs, key=lambda d: os.path.getmtime(os.path.join(s2_output_folder, d)))\n",
    "timestamp = datetime.now().strftime('%Y%m%d_%H%M%S')\n",
    "new_name = f\"S2_request_{timestamp}\"\n",
    "old_path = os.path.join(s2_output_folder, most_recent_dir)\n",
    "new_path = os.path.join(s2_output_folder, new_name)\n",
    "os.rename(old_path, new_path)\n",
    "beg = interval[0]\n",
    "end = interval [1]\n",
    "for files in os.listdir(new_path):\n",
    "    if files.__contains__('response') :\n",
    "        os.rename(os.path.join(new_path,files),os.path.join(new_path,f\"S2_{s2_res}m_{beg}_{end}.tif\"))\n",
    "        hr_img_path = os.path.join(new_path,f\"S2_{s2_res}m_{beg}_{end}.tif\") # This file will be used to train our model and downscale the ECOSTRESS image\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dowloading ECOSTRESS scenes\n",
    "\n",
    "This part is based off the ECOSTRESS API tutorials available [here](https://github.com/nasa/AppEEARS-Data-Resources/tree/main/Python/tutorials).\n",
    "\n",
    "Authentication and Token Retrieval for NASA AppEEARS API"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "api = 'https://appeears.earthdatacloud.nasa.gov/api/'  # Set the AρρEEARS API to a variable\n",
    "token_response = req.post('{}login'.format(api), auth=(user, password)).json() # Insert API URL, call login service, provide credentials & return json                                                        \n",
    "token = token_response['token']# Save login token to a variable\n",
    "head = {'Authorization': 'Bearer {}'.format(token)}  # Create a header to store token information, needed to submit a request\n",
    "print(token_response)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Above, you should see a Bearer token. Notice that this token will expire approximately 48 hours after being acquired. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Locate ECOSTRESS products and search for the layers of interest in our case : LST and QC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prods = ['ECO2LSTE.001','ECO_L2_LSTE.002']\n",
    "\n",
    "lst_response = req.get('{}product/{}'.format(api, prods[1])).json()  \n",
    "print(list(lst_response.keys()))\n",
    "if eco_collec == 1 :\n",
    "    layers = [(prods[0],'SDS_LST'),(prods[0],'SDS_QC')]\n",
    "elif eco_collec == 2 :\n",
    "    layers = [(prods[1],'LST'),(prods[1],'QC'),(prods[1],'cloud_mask')]\n",
    "\n",
    "# List of products and layers    \n",
    "prodLayer = []\n",
    "for l in layers: \n",
    "    prodLayer.append({\n",
    "            \"layer\": l[1],\n",
    "            \"product\": l[0]\n",
    "          })\n",
    "\n",
    "print(prodLayer)\n",
    "\n",
    "# Get the available projections list\n",
    "projections = req.get('{}spatial/proj'.format(api)).json()\n",
    "projs = {}                                 \n",
    "for p in projections: projs[p['Name']] = p  # Fill dictionary with `Name` of the projections available as keys\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Formulate the AppEEARS request"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Takes the bounding box entered above and creates a polygon for the ECOSTRESS request\n",
    "def polygon_from_bbox(bbox):\n",
    "    long0, lat0, long1, lat1 = bbox\n",
    "    return Polygon([[long0, lat0],\n",
    "                    [long1,lat0],\n",
    "                    [long1,lat1],\n",
    "                    [long0, lat1]])\n",
    "\n",
    "polyg = polygon_from_bbox(aoi_coords_wgs_84)\n",
    "\n",
    "aoi_df = gpd.GeoDataFrame(pd.DataFrame(columns = ['bbox']),\n",
    "        crs = 'epsg:4326',\n",
    "        geometry = [polyg])\n",
    "\n",
    "# Change the df to a json\n",
    "aoi_json = aoi_df.to_json()\n",
    "aoi_json = json.loads(aoi_json)\n",
    "\n",
    "# Create the actual task\n",
    "task_type = ['area']        # Type of task, area in our case\n",
    "proj = projs['geographic']['Name']  # Set output projection \n",
    "outFormat = ['geotiff', 'netcdf4']  # Set output file format type, we'll select geotiff here\n",
    "startDate = startMonth +'-'+startDay+'-'+startYear           # Start of the date range for which to extract data: MM-DD-YYYY\n",
    "endDate = endMonth +'-'+endDay+'-'+endYear               # End of the date range for which to extract data: MM-DD-YYYY\n",
    "\n",
    "\n",
    "task = {\n",
    "    'task_type': task_type[0],\n",
    "    'task_name': task_name,\n",
    "    'params': {\n",
    "         'dates': [\n",
    "         {\n",
    "             'startDate': startDate,\n",
    "             'endDate': endDate\n",
    "         }],\n",
    "         'layers': prodLayer,\n",
    "         'output': {\n",
    "                 'format': {\n",
    "                         'type': outFormat[0]}, \n",
    "                         'projection': proj},\n",
    "         'geo': aoi_json,\n",
    "    }\n",
    "}\n",
    "\n",
    "# Order the task\n",
    "task_response = req.post('{}task'.format(api), json=task, headers=head).json()                                                                  \n",
    "task_id = task_response['task_id']                                               \n",
    "status_response = req.get('{}status/{}'.format(api, task_id), headers=head).json() # Call status service with specific task ID & user credentials\n",
    "print(status_response)    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ping the API\n",
    "Ping the API every 30 seconds until the request is complete to display the status. It may occur that AppEEARS is slow to process data, I suggest to simply wait and come back later.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "starttime = time.time()\n",
    "while req.get('{}task/{}'.format(api, task_id), headers=head).json()['status'] != 'done':\n",
    "    print(req.get('{}task/{}'.format(api, task_id), headers=head).json()['status'])\n",
    "    time.sleep(30.0 - ((time.time() - starttime) % 30.0))\n",
    "print(req.get('{}task/{}'.format(api, task_id), headers=head).json()['status'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Download the ordered files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up output directory using input directory and task name\n",
    "eco_dest_dir = os.path.join(eco_output_folder,task_name)              \n",
    "if not os.path.exists(eco_dest_dir):\n",
    "    os.makedirs(eco_dest_dir) \n",
    "    \n",
    "bundle = req.get('{}bundle/{}'.format(api,task_id), headers=head).json()  # Call API and return bundle contents for the task_id as json\n",
    "\n",
    "files = {}                                                     \n",
    "for f in bundle['files']: files[f['file_id']] = f['file_name'] \n",
    "\n",
    "for f in files:\n",
    "    dl = req.get('{}bundle/{}/{}'.format(api, task_id, f), headers=head, stream=True, allow_redirects = 'True') # Get a stream to the bundle file\n",
    "    if files[f].endswith('.tif'):\n",
    "        filename = files[f].split('/')[1]\n",
    "    else:\n",
    "        filename = files[f] \n",
    "    filepath = os.path.join(eco_dest_dir, filename)                                                       \n",
    "    with open(filepath, 'wb') as f:                                                                  \n",
    "        for data in dl.iter_content(chunk_size=8192): f.write(data) \n",
    "print('Downloaded files can be found at: {}'.format(eco_dest_dir))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preprocessing ECOSTRESS Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sort the dowloaded files in appropriate subdirectories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For each directory, if it doesn't exist already, create it and place the dowloaded files in \n",
    "\n",
    "QC_dir = os.path.join(eco_dest_dir,'QC')\n",
    "if not os.path.exists(QC_dir) :\n",
    "    os.mkdir(QC_dir)\n",
    "lst_dir = os.path.join(eco_dest_dir,'LST')\n",
    "if not os.path.exists(lst_dir) :\n",
    "    os.mkdir(lst_dir)\n",
    "if eco_collec == 2 : \n",
    "    cld_dir = os.path.join(eco_dest_dir,'cloud_mask')\n",
    "    if not os.path.exists(cld_dir) :\n",
    "        os.mkdir(cld_dir)\n",
    "for file in os.listdir(eco_dest_dir) : \n",
    "    if file.__contains__('QC') and file.endswith('.tif') :\n",
    "        source_path = os.path.join(eco_dest_dir,file)\n",
    "        dst_path = os.path.join(QC_dir,file)\n",
    "        os.rename(source_path,dst_path)\n",
    "        \n",
    "    elif file.__contains__('LST_') and file.endswith('.tif') :\n",
    "        source_path = os.path.join(eco_dest_dir,file)\n",
    "        dst_path = os.path.join(lst_dir,file)\n",
    "        os.rename(source_path,dst_path)\n",
    "        \n",
    "    elif file.__contains__('cloud_mask') and file.endswith('.tif') :\n",
    "        source_path = os.path.join(eco_dest_dir,file)\n",
    "        dst_path = os.path.join(cld_dir,file)\n",
    "        os.rename(source_path,dst_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Preprocessing the QC files  \n",
    "\n",
    "The QC files are coded in 16 bits and thus can't be easily seen as a mask file. For convience, we write new Quality Flag (QF) files that respresent only the last two bits of the QC files. Then, there are only four possible values:  \n",
    "0 when the pixel is of best quality, 1 for nominal quality, 2 if a cloud is detected and 3 if the pixel is not produced. In the downscaling process, pixels with the last two values will be disregarded.   \n",
    "For more information on the QC files : https://lpdaac.usgs.gov/documents/423/ECO2_User_Guide_V1.pdf (section 2.4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for file in os.listdir(QC_dir) :\n",
    "    if not file.endswith('QF.tif') and not file.endswith('.xml') : \n",
    "        file_qc = os.path.join(QC_dir,file)\n",
    "        with rasterio.open(file_qc,'r') as f_qc :\n",
    "            # Read the QC file, coded in 16 bits\n",
    "            qc_img = f_qc.read((1)) \n",
    "            qc_img[qc_img==-99999] = -1  # Nodata values are read as -99999, we change it to -1 so that the last two bits appear as 11 (which means pixel not produced) and be masked out in the end\n",
    "            # Select only the last two bits\n",
    "            qc_img_2 = qc_img & 0b11 \n",
    "            out_meta = f_qc.meta.copy()\n",
    "        # Write the last two bits in a new file\n",
    "        file_qf = file_qc.replace('.tif','_QF.tif')\n",
    "        with rasterio.open(file_qf,'w',**out_meta) as dst : \n",
    "            dst.write(qc_img_2,1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scaling the ECOSTRESS LST to normal Kelvin scale.  \n",
    "  \n",
    "The LST product is actually scaled at 0.02, the GIS software takes that scale in account before display so you might not see it if you directly display on QGIS or ArcGIS. However, in Python it's easier to apply the scale rather than reading the metadata. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If the scaled subdirectory doesn't already exist, create it\n",
    "lst_dir_sc = os.path.join(eco_dest_dir,'LST_scaled')\n",
    "if not os.path.exists(lst_dir_sc) :\n",
    "        os.mkdir(lst_dir_sc)\n",
    "\n",
    "# Scale each file\n",
    "for file in os.listdir(lst_dir) : \n",
    "        if file.endswith('.tif') :\n",
    "                with rasterio.open(os.path.join(lst_dir,file),'r') as lr_img: \n",
    "                        out_image=lr_img.read().astype('float32')\n",
    "                        out_image[out_image==0]=np.nan\n",
    "                        out_meta = lr_img.meta\n",
    "\n",
    "                out_meta.update({\"driver\": \"GTiff\",\n",
    "                        \"height\": out_image.shape[1],\n",
    "                        \"width\": out_image.shape[2],\n",
    "                        'dtype' :'float32'})\n",
    "\n",
    "                dst_path = os.path.join(lst_dir_sc,file)\n",
    "                with rasterio.open(dst_path,'w',**out_meta) as dst :\n",
    "                        # Apply the scale \n",
    "                        dst.write(out_image*0.02 +0.49) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Upsampling using pyDMS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The preprocessing is now over. Let's sharpen using pyDMS. Use one of the following cells depending on the desired extent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you use this cell, then the output extent will be the extent of the orignal images downloaded."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "useDecisionTree = True # You could change this to False if you want to use the Neural Network intead of the Decision tree, not recommended\n",
    "files_sharpened = [] # list of the files sharpened\n",
    "dst_dir_f= os.path.join(dst_dir, f\"{task_name}-Results\")\n",
    "# Loop through the directory of LST images scaled\n",
    "for file in os.listdir(lst_dir_sc) :\n",
    "    if file.endswith('.tif') :\n",
    "        if not os.path.exists(dst_dir_f) : # create the output directory if it doesn't exist already\n",
    "                        os.mkdir(dst_dir_f)\n",
    "        outputFilename = os.path.join(dst_dir_f,file.replace('.tif','_sharp_S2.tif')) # destination path for the sharpened images\n",
    "        \n",
    "        highResFilename = hr_img_path # the high resolution file is the dowloaded s2 image\n",
    "        lowResFilename = os.path.join(lst_dir_sc,file) # the low resolution file is the scaled LST image\n",
    "        \n",
    "        valid = False # Bool that states if a sence is \"valid\", not too cloudy or not presenting too many unproduced pixels (limit at 25% by default)\n",
    "        thresh= 0.75 # you can modify this value between 0 (if you accept any file, the unusable pixels will be masked) and 1 (if you only want to sharpen files where every pixel is usable)\n",
    "        \n",
    "        # If the scene to be downscaled is a scene from Collection 2\n",
    "        if eco_collec == 2 : \n",
    "            # The QC cloud bit being unreliable in Collection 2, we use the cloud mask directly as a validity mask\n",
    "            file_cl = 'cloud_mask'.join(file.rsplit('LST', 1))\n",
    "            lowResMaskFilename = os.path.join(cld_dir,file_cl)\n",
    "            lowresflags = [0]\n",
    "            with rasterio.open(lowResMaskFilename,'r') as cld_msk : \n",
    "                cld_msk_arr = cld_msk.read(1)\n",
    "                mask_sz = cld_msk_arr.size\n",
    "                # Ensure that the scene to be sharpened isn't more than 25% cloudy\n",
    "            if np.count_nonzero(cld_msk_arr)<(1-thresh)*mask_sz :\n",
    "                valid = True\n",
    "        # If the scene to be downscaled belongs to Collection 1        \n",
    "        else : \n",
    "            # For Collection 1, we can rely on the QC (preprocessed earlier) file, no need to use the cloud mask\n",
    "            file_qc = 'QC'.join(file.rsplit('LST', 1))\n",
    "            file_qf = file_qc.replace('.tif','_QF.tif')\n",
    "            lowResMaskFilename = os.path.join(QC_dir,file_qf)        \n",
    "            with rasterio.open(lowResMaskFilename,'r') as mask :\n",
    "                mask_array = mask.read(1)\n",
    "                mask_sz = mask_array.size \n",
    "            lowresflags = [0,1]\n",
    "            # Ensure that the scene to be sharpened isn't more than 25% cloudy or invalid\n",
    "            if np.count_nonzero(mask_array==0) + np.count_nonzero(mask_array==1)>thresh*mask_sz :\n",
    "                valid = True\n",
    "        # Only downscale the files that are \"valid\" \n",
    "        if valid :    \n",
    "            commonOpts = {\"highResFiles\":               [highResFilename],\n",
    "                            \"lowResFiles\":                [lowResFilename],\n",
    "                            \"lowResQualityFiles\":         [lowResMaskFilename],\n",
    "                            \"lowResGoodQualityFlags\":     lowresflags, # flags for acceptable pixels\n",
    "                            \"cvHomogeneityThreshold\":     0, # the homogeneity threshold will be automatically computed\n",
    "                            \"movingWindowSize\":           0,\n",
    "                            \"disaggregatingTemperature\":  True} # we are dealing with LST\n",
    "            dtOpts =     {\"perLeafLinearRegression\":    True,\n",
    "                            \"linearRegressionExtrapolationRatio\": 0.25} # how much extrapolation we accept from the tree\n",
    "            sknnOpts =   {'hidden_layer_sizes':         (10,),\n",
    "                            'activation':                 'tanh'}\n",
    "            nnOpts =     {\"regressionType\":             REG_sklearn_ann,\n",
    "                            \"regressorOpt\":               sknnOpts}\n",
    "            \n",
    "            start_time = time.time()\n",
    "\n",
    "            if useDecisionTree: #  Regression tree\n",
    "                opts = commonOpts.copy()\n",
    "                opts.update(dtOpts)\n",
    "                disaggregator = DecisionTreeSharpener(**opts)\n",
    "            else: # Neural network\n",
    "                opts = commonOpts.copy()\n",
    "                opts.update(nnOpts)\n",
    "                disaggregator = NeuralNetworkSharpener(**opts)\n",
    "\n",
    "            # Training the tree\n",
    "            print(\"Training regressor...\")\n",
    "            disaggregator.trainSharpener()\n",
    "            # Applying the newly trained tree to the LR image\n",
    "            print(\"Sharpening...\")\n",
    "            downscaledFile = disaggregator.applySharpener(highResFilename, lowResFilename)\n",
    "            # Performing residual analysis, to ensure that the unsharpened images's exitance is conserved\n",
    "            print(\"Residual analysis...\")\n",
    "            residualImage, correctedImage = disaggregator.residualAnalysis(downscaledFile, lowResFilename,\n",
    "                                                                            lowResMaskFilename,\n",
    "                                                                            doCorrection=True)\n",
    "            # Saving the image and its residual \n",
    "            print(\"Saving output...\")\n",
    "            highResFile = gdal.Open(highResFilename)\n",
    "            if correctedImage is not None:\n",
    "                outImage = correctedImage\n",
    "            else:\n",
    "                outImage = downscaledFile\n",
    "\n",
    "            outFile = utils.saveImg(outImage.GetRasterBand(1).ReadAsArray(),\n",
    "                                    outImage.GetGeoTransform(),\n",
    "                                    outImage.GetProjection(),\n",
    "                                    outputFilename)\n",
    "            residualFile = utils.saveImg(residualImage.GetRasterBand(1).ReadAsArray(),\n",
    "                                            residualImage.GetGeoTransform(),\n",
    "                                            residualImage.GetProjection(),\n",
    "                                            os.path.splitext(outputFilename)[0] + \"_residual\" +\n",
    "                                            os.path.splitext(outputFilename)[1])\n",
    "            files_sharpened.append(file)\n",
    "            outFile = None\n",
    "            residualFile = None\n",
    "            downsaceldFile = None\n",
    "            highResFile = None\n",
    "\n",
    "            print(time.time() - start_time, \"seconds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you wish to sharpen only part of the image, for instance just the center of a city and not the suburbs, then use this cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "useDecisionTree = True # You could change this to False if you want to use the Neural Network intead of the Decision tree, not recommended\n",
    "\n",
    "files_sharpened = [] #list of the files sharpened\n",
    "dst_dir_f = os.path.join(dst_dir, f\"{task_name}-Results\")\n",
    "# Loop through the directory of LST images\n",
    "for file in os.listdir(lst_dir_sc) :\n",
    "    if file.endswith('.tif') :\n",
    "        if not os.path.exists(dst_dir_f) :  # create the output directory if it doesn't exist already\n",
    "                        os.mkdir(dst_dir_f)\n",
    "        outputFilename = os.path.join(dst_dir_f,file.replace('.tif','_sharp_S2_clipped.tif')) # destination path for the sharpened images\n",
    "        lowResFilename = os.path.join(lst_dir_sc,file) \n",
    "        lr_ds = rasterio.open(lowResFilename)\n",
    "        \n",
    "        # Make sure the coordinates here are in the correct order, otherwise you'll recieve an error. This extent has to be fully included in the bounding box.\n",
    "        # projwin =  [xmin,ymax,xmax,ymin]\n",
    "        projwin = [1.336555,43.650982,1.503754,43.549543] # example of a custom cutout  \n",
    "\n",
    "        hr_img_path_clipped = hr_img_path.replace('.tif',f'_clipped.tif') # also produces a clipped version of the HR image\n",
    "        ds = gdal.Open(hr_img_path)\n",
    "        ds = gdal.Translate(hr_img_path_clipped, ds, projWin = projwin)\n",
    "        ds = None\n",
    "        lr_ds = None\n",
    "        \n",
    "        highResFilename = hr_img_path_clipped \n",
    "        \n",
    "        valid = False # Bool that states if a sence is \"valid\", not too cloudy or not presenting too many unproduced pixels (limit at 25% by default)\n",
    "        thresh= 0.75 # you can modify this value between 0 (if you accept any file, the unusable pixels will be masked) and 1 (if you only want to sharpen files where every pixel is usable)\n",
    "        \n",
    "        # If the scene to be downscaled is a scene from Collection 2\n",
    "        if eco_collec == 2 : \n",
    "            # The QC cloud bit being unreliable in Collection 2, we use the cloud mask directly as a validity mask\n",
    "            file_cl = 'cloud_mask'.join(file.rsplit('LST', 1))\n",
    "            lowResMaskFilename = os.path.join(cld_dir,file_cl)\n",
    "            lowresflags = [0]\n",
    "            with rasterio.open(lowResMaskFilename,'r') as cld_msk : \n",
    "                cld_msk_arr = cld_msk.read(1)\n",
    "                mask_sz = cld_msk_arr.size\n",
    "                # Ensure that the scene to be sharpened isn't more than 25% cloudy\n",
    "            if np.count_nonzero(cld_msk_arr)<(1-thresh)*mask_sz :\n",
    "                valid = True\n",
    "        # If the scene to be downscaled belongs to Collection 1        \n",
    "        else : \n",
    "            # For Collection 1, we can rely on the QC (preprocessed earlier) file, no need to use the cloud mask\n",
    "            file_qc = 'QC'.join(file.rsplit('LST', 1))\n",
    "            file_qf = file_qc.replace('.tif','_QF.tif')\n",
    "            lowResMaskFilename = os.path.join(QC_dir,file_qf)        \n",
    "            with rasterio.open(lowResMaskFilename,'r') as mask :\n",
    "                mask_array = mask.read(1)\n",
    "                mask_sz = mask_array.size \n",
    "            lowresflags = [0,1]\n",
    "            # Ensure that the scene to be sharpened isn't more than 25% cloudy or invalid\n",
    "            if np.count_nonzero(mask_array==0) + np.count_nonzero(mask_array==1)>thresh*mask_sz :\n",
    "                valid = True\n",
    "        # Only downscale the files that are \"valid\" four our use case\n",
    "        if valid :\n",
    "            commonOpts = {\"highResFiles\":               [highResFilename],\n",
    "                            \"lowResFiles\":                [lowResFilename],\n",
    "                            \"lowResQualityFiles\":         [lowResMaskFilename],\n",
    "                            \"lowResGoodQualityFlags\":     lowresflags, # flags for acceptable pixels\n",
    "                            \"cvHomogeneityThreshold\":     0, # the homogeneity threshold will be automatically computed\n",
    "                            \"movingWindowSize\":           0,\n",
    "                            \"disaggregatingTemperature\":  True} # we are dealing with LST\n",
    "            dtOpts =     {\"perLeafLinearRegression\":    True,\n",
    "                            \"linearRegressionExtrapolationRatio\": 0.25} # how much extrapolation we accept from the tree\n",
    "            sknnOpts =   {'hidden_layer_sizes':         (10,),\n",
    "                            'activation':                 'tanh'}\n",
    "            nnOpts =     {\"regressionType\":             REG_sklearn_ann,\n",
    "                            \"regressorOpt\":               sknnOpts}\n",
    "            \n",
    "            commonOpts = {\"highResFiles\":               [highResFilename],\n",
    "                            \"lowResFiles\":                [lowResFilename],\n",
    "                            \"lowResQualityFiles\":         [lowResMaskFilename],\n",
    "                            \"lowResGoodQualityFlags\":     [0,1], # flags for acceptable pixels\n",
    "                            \"cvHomogeneityThreshold\":     0, # the homogeneity threshold will be automatically computed\n",
    "                            \"movingWindowSize\":           0,\n",
    "                            \"disaggregatingTemperature\":  True} # we are dealing with LST\n",
    "            dtOpts =     {\"perLeafLinearRegression\":    True,\n",
    "                            \"linearRegressionExtrapolationRatio\": 0.25} # how much extrapolation we accept from the tree\n",
    "            sknnOpts =   {'hidden_layer_sizes':         (10,),\n",
    "                            'activation':                 'tanh'}\n",
    "            nnOpts =     {\"regressionType\":             REG_sklearn_ann,\n",
    "                            \"regressorOpt\":               sknnOpts}\n",
    "            \n",
    "            start_time = time.time()\n",
    "\n",
    "            if useDecisionTree: #  Regression tree\n",
    "                opts = commonOpts.copy()\n",
    "                opts.update(dtOpts)\n",
    "                disaggregator = DecisionTreeSharpener(**opts)\n",
    "            else: # Neural network\n",
    "                opts = commonOpts.copy()\n",
    "                opts.update(nnOpts)\n",
    "                disaggregator = NeuralNetworkSharpener(**opts)\n",
    "\n",
    "            # Training the tree\n",
    "            print(\"Training regressor...\")\n",
    "            disaggregator.trainSharpener()\n",
    "            # Applying the newly trained tree to the LR image\n",
    "            print(\"Sharpening...\")\n",
    "            downscaledFile = disaggregator.applySharpener(highResFilename, lowResFilename)\n",
    "            # Performing residual analysis, to ensure that the unsharpened images's exitance is conserved\n",
    "            print(\"Residual analysis...\")\n",
    "            residualImage, correctedImage = disaggregator.residualAnalysis(downscaledFile, lowResFilename,\n",
    "                                                                            lowResMaskFilename,\n",
    "                                                                            doCorrection=True)\n",
    "            # Saving the image and its residual \n",
    "            print(\"Saving output...\")\n",
    "            highResFile = gdal.Open(highResFilename)\n",
    "            if correctedImage is not None:\n",
    "                outImage = correctedImage\n",
    "            else:\n",
    "                outImage = downscaledFile\n",
    "\n",
    "            outFile = utils.saveImg(outImage.GetRasterBand(1).ReadAsArray(),\n",
    "                                    outImage.GetGeoTransform(),\n",
    "                                    outImage.GetProjection(),\n",
    "                                    outputFilename)\n",
    "            residualFile = utils.saveImg(residualImage.GetRasterBand(1).ReadAsArray(),\n",
    "                                            residualImage.GetGeoTransform(),\n",
    "                                            residualImage.GetProjection(),\n",
    "                                            os.path.splitext(outputFilename)[0] + \"_residual\" +\n",
    "                                            os.path.splitext(outputFilename)[1])\n",
    "            files_sharpened.append(file)\n",
    "            outFile = None\n",
    "            residualFile = None\n",
    "            downsaceldFile = None\n",
    "            highResFile = None\n",
    "\n",
    "            print(time.time() - start_time, \"seconds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I adivse to check the residual bias and RMSD to analyze the perfomance of the downscaling. It performs better on images that are not nosiy, clear sky and ideally a quite homogenous envrionment.Any RMSD above 5 is really suspicious and usually means that the training was bad. From my observations, unsmasked clouds often explain a poor performance. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Display"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot one random sharpened image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pick one random sharpened file\n",
    "file = random.choice(files_sharpened)\n",
    "sharpened_file = os.path.join(dst_dir_f,file.replace('.tif','_sharp_S2.tif'))\n",
    "raw_file = os.path.join(lst_dir_sc,file)\n",
    "with rasterio.open(raw_file,'r') as lr_img : \n",
    "    raw_lst = lr_img.read(1)\n",
    "with rasterio.open(sharpened_file,'r') as shrp_img :\n",
    "    shrp_lst = shrp_img.read(1)\n",
    "    \n",
    "# Plot the ECOSTRESS original product\n",
    "plt.figure(1)\n",
    "plt.imshow(raw_lst,cmap='viridis')\n",
    "plt.axis('off')\n",
    "plt.colorbar(label ='LST(K)')\n",
    "vmin, vmax = plt.gci().get_clim() # save the limits to share them in the second figure\n",
    "plt.title(\"ECOSTRESS LST 70m\")\n",
    "plt.show()\n",
    "\n",
    "# Plot the ECOSTRESS sharpened product\n",
    "plt.figure(2)\n",
    "plt.imshow(shrp_lst,cmap='viridis',clim =(vmin,vmax))\n",
    "plt.axis('off')\n",
    "plt.colorbar(label ='LST(K)')\n",
    "plt.title(\"ECOSTRESS LST sharpened to 30m\")\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ECOSTRESS",
   "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.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}