full_pystac_workflow.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6da2ce6a-f70f-46d1-a1c8-689dd3dfd98c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The following codes to import, clean, and process Sentinel-2 data using VSI and AWS were written by Lilly Jones and Erick Verleye; edited by Ty Tuff, pseudocode by Cibele Amaral. \n",
    "\n",
    "#imports\n",
    "\n",
    "import os\n",
    "import glob\n",
    "from sentinelsat import SentinelAPI, read_geojson, geojson_to_wkt\n",
    "from datetime import date\n",
    "from osgeo import gdal, gdalconst\n",
    "import numpy as np\n",
    "import rasterio\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91f2f5db-d07d-49ad-8df1-aa98644a9f83",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# 1. Download images using Sentinel-2 API\n",
    "\n",
    "# download Sentinel-2 data by date and bounding box\n",
    "\n",
    "def download_sentinel_data(username, password, start_date, end_date, bounding_box, output_directory):\n",
    "    \"\"\"\n",
    "    Downloads Sentinel-2 data using the Copernicus API Hub.\n",
    "\n",
    "    Args:\n",
    "        username (str): Your Copernicus Open Access Hub username.\n",
    "        password (str): Your Copernicus Open Access Hub password.\n",
    "        start_date (str): Start date in 'YYYYMMDD' format (e.g., '20220101').\n",
    "        end_date (str): End date in 'YYYYMMDD' format (e.g., '20220131').\n",
    "        bounding_box (list): Bounding box coordinates in the format [lon_min, lat_min, lon_max, lat_max].\n",
    "        output_directory (str): Directory where downloaded data will be stored.\n",
    "\n",
    "    Returns:\n",
    "        None\n",
    "    \"\"\"\n",
    "    # Initialize the Sentinel API with your credentials\n",
    "    api = SentinelAPI(username, password, 'https://apihub.copernicus.eu/dhus')\n",
    "\n",
    "    # Convert the bounding box to Well-Known Text (WKT) format\n",
    "    footprint = geojson_to_wkt({\n",
    "        \"type\": \"Polygon\",\n",
    "        \"coordinates\": [bounding_box]\n",
    "    })\n",
    "\n",
    "    # Search for Sentinel-2 products within the specified date range and bounding box\n",
    "    products = api.query(footprint,\n",
    "                         date=(start_date, end_date),\n",
    "                         platformname='Sentinel-2',\n",
    "                         cloudcoverpercentage=(0, 30))  # You can adjust the cloud cover percentage as needed\n",
    "\n",
    "    # Download the products\n",
    "    for product_id, product_info in products.items():\n",
    "        print(f\"Downloading product {product_id}...\")\n",
    "        api.download(product_id, directory_path=output_directory)\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    # Set your Copernicus Open Access Hub credentials\n",
    "    username = \"your_username\"\n",
    "    password = \"your_password\"\n",
    "\n",
    "    # Define the date range and bounding box for your area of interest\n",
    "    start_date = \"20220101\"\n",
    "    end_date = \"20220131\"\n",
    "    bounding_box = [lon_min, lat_min, lon_max, lat_max]  # Replace with your bounding box coordinates\n",
    "    output_directory = \"path_to_output_directory\"  # Replace with your desired output directory\n",
    "\n",
    "    # Call the download_sentinel_data function\n",
    "    download_sentinel_data(username, password, start_date, end_date, bounding_box, output_directory)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "310127a5-b18d-4745-b171-ebd0507085a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2. Resample bands to 10-m resolution\n",
    "- Select bands 2 to 8A and 11 to 12 (10 bands total)\n",
    "- Resample the 20-m resolution bands to 10-m resolution \n",
    "- Stack the 10 bands into one raster\n",
    "- Gdal (https://gdal.org/api/python_bindings.html) does everything listed above\n",
    "\n",
    "import gdal\n",
    "from gdalconst import *\n",
    "\n",
    "# Input file paths\n",
    "input_file_20m = 'path_to_input_20m_image.tif'\n",
    "input_file_10m = 'path_to_input_10m_image.tif'\n",
    "\n",
    "# Output file path\n",
    "output_file = 'path_to_output_image.tif'\n",
    "\n",
    "# Open the 20m resolution raster\n",
    "ds_20m = gdal.Open(input_file_20m, GA_ReadOnly)\n",
    "\n",
    "# Resample the 20m resolution bands to 10m resolution\n",
    "ds_10m = gdal.Warp('', ds_20m, xRes=10, yRes=10)\n",
    "\n",
    "# Select bands 2 to 8A and 11 to 12 (10 bands total)\n",
    "band_selection = [2, 3, 4, 5, 6, 7, 8, 11, 12]\n",
    "\n",
    "# Stack the selected bands into one raster\n",
    "output_ds = gdal.BuildVRT('', [ds_10m.GetRasterBand(band) for band in band_selection])\n",
    "\n",
    "# Write the stacked raster to an output file\n",
    "gdal.Translate(output_file, output_ds, format='GTiff')\n",
    "\n",
    "# Clean up\n",
    "ds_20m = None\n",
    "ds_10m = None\n",
    "output_ds = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "746d4414-7342-4ae7-93f9-3493d05086f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Original AWS code by Erick Verleye, ESIIL Software Developer 2023-08-08, lightly edited by Lilly Jones for SCE project\n",
    "# Code downloads Sentinel Level-1 C data from S3 with Python for a given latitude, longitude, and date range.\n",
    "# An AWS account is needed and will be charged (a small amount) for any data downloaded.\n",
    "\n",
    "#imports\n",
    "\n",
    "import json\n",
    "\n",
    "import multiprocessing as mp\n",
    "import os\n",
    "\n",
    "from argparse import Namespace\n",
    "from datetime import datetime, timedelta\n",
    "from typing import List, Tuple, Union\n",
    "\n",
    "import boto3\n",
    "import geojson\n",
    "import tqdm\n",
    "from shapely.geometry import Polygon\n",
    "\n",
    "# Use boto3 to create a connection to the AWS account. Initialization of the boto3.session object changes depending on the environment you are running this code in:\n",
    "\n",
    "# Using the AWS CLI from a personal computer to log in, the profile name will typically be default.\n",
    "# Using a federated access login, the profile name will typically be the name of the software.\n",
    "# Using S3 configured IAM profile, no profile name is required.\n",
    "# Each of these arguments, if applicable, can be found in the aws credentials file typically found at ~/.aws/credentials\n",
    "\n",
    "\n",
    "# Replace with your AWS profile name if applicable\n",
    "profile = 'default'\n",
    "\n",
    "session = boto3.Session(profile_name=profile)\n",
    "\n",
    "# AWS CLI\n",
    "session = boto3.Session(profile_name='default')\n",
    "\n",
    "# Federated login (saml2aws example)\n",
    "session = boto3.Session(profile_name='saml')\n",
    "\n",
    "# S3 configured IAM profile on EC2 instance\n",
    "session = boto3.Session()\n",
    "\n",
    "# Define constants\n",
    "\n",
    "SENTINEL_2_BUCKET_NAME = 'insert_name_of_S3_bucket'  # Name of the s3 bucket on AWS hosting the sentinel-2 data\n",
    "\n",
    "\n",
    "# find_overlapping_mgrs Sentinel-2 data is organized depending on which Military Grid Reference System (MGRS) square that it belongs to. This function converts a bounding box defined in    # lat/lon as [min_lon, min_lat, max_lon, max_lat] to the military grid squares that is overlaps. More information on # # the MGRS can be found at                                                       # https://www.maptools.com/tutorials/mgrs/quick_guide. NOTE: You will have to download the mgrs_lookup.geojson file from # # # # # https://github.com/CU-ESIIL/data-                        # library/blob/main/docs/remote_sensing/sentinel2_aws/mgrs_lookup.geojson and place it into the working directory that # # this code is being run from.\n",
    "\n",
    "def find_overlapping_mgrs(bounds: List[float]) -> List[str]:\n",
    "    \"\"\"\n",
    "    Files in the Sinergise Sentinel2 S3 bucket are organized by which military grid they overlap. Thus, the\n",
    "    military grids that overlap the input bounding box defined in lat / lon must be found. A lookup table that\n",
    "    includes each grid name and its geometry is used to find the overlapping grids.\n",
    "    Args:\n",
    "        bounds (list): Bounding box definition as [min_lon, min_lat, max_lon, max_lat]\n",
    "    \"\"\"\n",
    "    print('Finding overlapping tiles... ')\n",
    "    input_bounds = Polygon([(bounds[0], bounds[1]), (bounds[2], bounds[1]), (bounds[2], bounds[3]),\n",
    "                            (bounds[0], bounds[3]), (bounds[0], bounds[1])])\n",
    "    with open('mgrs_lookup.geojson', 'r') as f:\n",
    "        ft = geojson.load(f)\n",
    "        return [i['properties']['mgs'] for i in ft[1:] if\n",
    "                input_bounds.intersects(Polygon(i['geometry']['coordinates'][0]))]\n",
    "\n",
    "# find_available_files finds the set of available files in the s3 bucket given a date range, lat/lon bounds, and list of bands.\n",
    "\n",
    "def find_available_files(s3_client, bounds: List[float], start_date: datetime, end_date: datetime,\n",
    "                         bands: List[str]) -> List[Tuple[str, str]]:\n",
    "    \"\"\"\n",
    "    Given a bounding box and start / end date, finds which files are available on the bucket and\n",
    "    meet the search criteria\n",
    "    Args:\n",
    "        bounds (list): Lower left and top right corner of bounding box defined in\n",
    "        lat / lon [min_lon, min_lat, max_lon, max_lat]\n",
    "        start_date (str): Beginning of requested data creation date YYYY-MM-DD\n",
    "        end_date (str): End of requested data creation date YYYY-MM-DD\n",
    "    \"\"\"\n",
    "    date_paths = []\n",
    "    ref_date = start_date\n",
    "    while ref_date <= end_date:\n",
    "        tt = ref_date.timetuple()\n",
    "        date_paths.append(f'/{tt.tm_year}/{tt.tm_mon}/{tt.tm_mday}/')\n",
    "        ref_date = ref_date + timedelta(days=1)\n",
    "\n",
    "    info = []\n",
    "    mgrs_grids = find_overlapping_mgrs(bounds)\n",
    "    for grid_string in mgrs_grids:\n",
    "        utm_code = grid_string[:2]\n",
    "        latitude_band = grid_string[2]\n",
    "        square = grid_string[3:5]\n",
    "        grid = f'tiles/{utm_code}/{latitude_band}/{square}'\n",
    "        response = s3_client.list_objects_v2(\n",
    "            Bucket=SENTINEL_2_BUCKET_NAME,\n",
    "            Prefix=grid + '/',\n",
    "            MaxKeys=300,\n",
    "            RequestPayer='requester'\n",
    "        )\n",
    "        if 'Contents' not in list(response.keys()):\n",
    "            continue\n",
    "\n",
    "        for date in date_paths:\n",
    "            response = s3_client.list_objects_v2(\n",
    "                Bucket=SENTINEL_2_BUCKET_NAME,\n",
    "                Prefix=grid + date + '0/',  # '0/' is for the sequence, which in most cases will be 0\n",
    "                MaxKeys=100,\n",
    "                RequestPayer='requester'\n",
    "            )\n",
    "            if 'Contents' in list(response.keys()):\n",
    "                info += [\n",
    "                    (v['Key'], v['Size']) for v in response['Contents'] if\n",
    "                    any([band + '.jp2' in v['Key'] for band in bands]) or 'MSK_CLOUDS_B00.gml' in v['Key']\n",
    "                ]\n",
    "\n",
    "    return info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0add4cd6-246c-4e83-bcff-4f83e1a8d098",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This function is designed to download the files in parallel, so a download_task function is defined as well. \n",
    "# Note that multiprocessing will not work as implemented here in iPython (Jupyter # # notebooks, etc.) and so that block of code is commented out. \n",
    "# If you are running this in a different environment and would like to download in parallel, un-comment this block and comment the sequential block below.\n",
    "\n",
    "\n",
    "def download_task(namespace: Namespace) -> None:\n",
    "    \"\"\"\n",
    "    Downloads a single file from the indicated s3 bucket. This function is intended to be spawned in parallel from the\n",
    "    parent process.\n",
    "    Args:\n",
    "        namespace (Namespace): Contains the bucket name, s3 file name, and destination required for s3 download request.\n",
    "        Each value in the namespace must be a pickle-izable type (i.e. str).\n",
    "    \"\"\"\n",
    "    session = boto3.Session(profile_name=namespace.profile_name)\n",
    "    s3 = session.client('s3')\n",
    "    s3.download_file(namespace.bucket_name, namespace.available_file,\n",
    "                     namespace.dest,\n",
    "                     ExtraArgs={'RequestPayer': 'requester'}\n",
    "                     )\n",
    "\n",
    "\n",
    "def download(session, bounds: List[float], start_date: datetime, end_date: datetime,\n",
    "             bands: List[float], out_dir: str, buffer: float = None) -> None:\n",
    "    \"\"\"\n",
    "    Downloads a list of .jp2 files from the Sinergise Sentinel2 LC1 bucket given a bounding box defined in lat/long, a buffer in meters, and a start and end date. Only Bands 2-4 are requested.\n",
    "     Args:\n",
    "         bounds (list): Bounding box defined in lat / lon [min_lon, min_lat, max_lon, max_lat]\n",
    "         buffer (float): Amount by which to extend the bounding box by, in meters\n",
    "         start_date (str): Beginning of requested data creation date YYYY-MM-DD\n",
    "         end_date (str): End of requested data creation date YYYY-MM-DD\n",
    "         bands (list): The bands to download for each file. Ex. ['B02', 'B03', 'B04', 'B08'] for R, G, B, and\n",
    "         near wave IR, respectively\n",
    "         out_dir (str): Path to directory where downloaded files will be written to\n",
    "    \"\"\"\n",
    "    # Convert the buffer from meters to degrees lat/long at the equator\n",
    "    if buffer is not None:\n",
    "        buffer /= 111000\n",
    "\n",
    "        # Adjust the bounding box to include the buffer (subtract from min lat/long values, add to max lat/long values)\n",
    "        bounds[0] -= buffer\n",
    "        bounds[1] -= buffer\n",
    "        bounds[2] += buffer\n",
    "        bounds[3] += buffer\n",
    "\n",
    "    s3_client = session.client('s3')\n",
    "    available_files = find_available_files(s3_client, bounds, start_date, end_date, bands)\n",
    "\n",
    "    args = []\n",
    "    total_data = 0\n",
    "    for file_info in available_files:\n",
    "        file_path = file_info[0]\n",
    "        if '/preview/' in file_path:\n",
    "            continue\n",
    "\n",
    "        created_file_path = os.path.join(out_dir, file_path.replace('_qi_', '').replace('/', '_').replace('tiles_', ''))\n",
    "\n",
    "        # Skip if file is already local\n",
    "        if os.path.exists(created_file_path):\n",
    "            continue\n",
    "\n",
    "        total_data += file_info[1]\n",
    "\n",
    "        args.append(Namespace(available_file=file_path, bucket_name=SENTINEL_2_BUCKET_NAME, profile_name=session.profile_name,\n",
    "                              dest=created_file_path))\n",
    "\n",
    "    total_data /= 1E9\n",
    "    print(f'Found {len(args)} files for download. Total size of files is'\n",
    "          f' {round(total_data, 2)}GB and estimated cost will be ${round(0.09 * total_data, 2)}'\n",
    "          )\n",
    "\n",
    "    # For multiprocessing when being run in iPython (Jupyter notebook, etc.)\n",
    "    # with mp.Pool(mp.cpu_count() - 1) as pool:\n",
    "        # for _ in tqdm.tqdm(pool.imap_unordered(download_task, args), total=len(args)):\n",
    "            # pass\n",
    "\n",
    "    # Sequential download for use in iPython (Jupyter notebooks, etc.)\n",
    "    for _ in tqdm.tqdm(args, total=len(args)):\n",
    "        download_task(args)\n",
    "\n",
    "# Once all of the above functions and variables are defined, download for a specific time range and region. \n",
    "# Make sure the directory corresponding to the out_dir parameter you pass in already exists. If not, create the directory.\n",
    "\n",
    "# Define the bounding box for California [min_lon, min_lat, max_lon, max_lat]\n",
    "sce_bounds = [-124.482003, 32.528832, -114.130775, 42.009518]\n",
    "\n",
    "# Define the output directory\n",
    "out_dir = 'SCE'\n",
    "os.makedirs(out_dir, exist_ok=True)\n",
    "\n",
    "\n",
    "# Download RGB band data for California from a specific date range\n",
    "download(session=session, bounds=california_bounds, start_date=datetime(YYYY, MM, DD),  # Replace with your date range\n",
    "         end_date=datetime(YYYY, MM, DD),  # Replace with your date range\n",
    "         bands=['B02', 'B03', 'B04'], out_dir=out_dir)\n",
    "\n",
    "\n",
    "# Cloud files will be downloaded alongside the band data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4680d0e-fc94-46c9-87bb-e9595b47a02f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The following cloud probability mask code is by Erick Verleye\n",
    "\n",
    "# Define constants\n",
    "MAX_BAND_VAL = 4000  # counts\n",
    "\n",
    "# Helper function for opening Sentinel-2 jp2 files, with optional slicing\n",
    "\n",
    "def get_img_from_file(img_path, g_ncols, dtype, row_bound=None):\n",
    "    img = rasterio.open(img_path, driver='JP2OpenJPEG')\n",
    "    ncols, nrows = img.meta['width'], img.meta['height']\n",
    "    assert g_ncols == ncols, f'Imgs have different size ncols: {ncols} neq {g_ncols}'\n",
    "    if row_bound is None:\n",
    "        pixels = img.read(1).astype(np.float32)\n",
    "    else:\n",
    "        pixels = img.read(\n",
    "            1,\n",
    "            window=Window.from_slices(\n",
    "                slice(row_bound[0], row_bound[1]),\n",
    "                slice(0, ncols)\n",
    "            )\n",
    "        ).astype(dtype)\n",
    "    return pixels\n",
    "\n",
    "OR: def get_img_from_file(img_path, g_ncols, dtype, row_bound=None):\n",
    "    img = rasterio.open(img_path, driver='JP2OpenJPEG')\n",
    "    ncols, nrows = img.meta['width'], img.meta['height']\n",
    "    assert g_ncols == ncols, f'Imgs have different size ncols: {ncols} neq {g_ncols}'\n",
    "    if row_bound is None:\n",
    "        pixels = img.read(1).astype(np.float32)\n",
    "    else:\n",
    "        pixels = img.read(\n",
    "            1,\n",
    "            window=Window.from_slices(\n",
    "                slice(row_bound[0], row_bound[1]),\n",
    "                slice(0, ncols)\n",
    "            )\n",
    "        ).astype(dtype)\n",
    "    return pixels\n",
    "\n",
    "# Helper function for reading in cloud file array, with optional slicing\n",
    "def get_cloud_mask_from_file(cloud_path, crs, transform, shape, row_bound=None):\n",
    "    # filter out RuntimeWarnings, due to geopandas/fiona read file spam\n",
    "    # https://stackoverflow.com/questions/64995369/geopandas-warning-on-read-file\n",
    "    warnings.filterwarnings(\"ignore\", category=RuntimeWarning)\n",
    "    try:\n",
    "        cloud_file = gpd.read_file(cloud_path)\n",
    "        cloud_file.crs = (str(crs))\n",
    "        # convert the cloud mask data to a raster that has the same shape and transformation as the\n",
    "        # img raster data\n",
    "        cloud_img = features.rasterize(\n",
    "            (\n",
    "                (g['geometry'], 1) for v, g in cloud_file.iterrows()\n",
    "            ),\n",
    "            out_shape=shape,\n",
    "            transform=transform,\n",
    "            all_touched=True\n",
    "        )\n",
    "        if row_bound is None:\n",
    "            return np.where(cloud_img == 0, 1, 0)\n",
    "        return np.where(cloud_img[row_bound[0]:row_bound[1], :] == 0, 1, 0)\n",
    "    except Exception as e:\n",
    "        return None\n",
    "\n",
    "# Function for filtering out cloud pixels\n",
    "\n",
    "def nan_clouds(pixels, cloud_channels, max_pixel_val: float = MAX_BAND_VAL):\n",
    "    cp = pixels * cloud_channels\n",
    "    mask = np.where(np.logical_or(cp == 0, cp > max_pixel_val))\n",
    "    cp[mask] = np.nan\n",
    "    return cp\n",
    "\n",
    "OR: def nan_clouds(pixels, cloud_channels, max_pixel_val: float = MAX_BAND_VAL):\n",
    "    cp = pixels * cloud_channels\n",
    "    mask = np.where(np.logical_or(cp == 0, cp > max_pixel_val))\n",
    "    cp[mask] = np.nan\n",
    "    return cp\n",
    "\n",
    "# Function to create cloud-cleaned composite for a single military tile and optical band\n",
    "\n",
    "def create_cloud_cleaned_composite(in_dir: str, mgrs_tile: str, band: str, out_file: str, num_slices: int = 12) -> None:\n",
    "    \"\"\"\n",
    "    Creates a cloud cleaned composite tif file from a set of sentinel 2 files\n",
    "    Args:\n",
    "        in_dir (str): Directory containing all of the date directories which contain the sentinel 2 data\n",
    "        mgrs_tile (str): The military grid coordinate (35MGR, 36MTV, etc.) to create the composite for\n",
    "        band (str): The optical band for which to create the composite for (B02, B03, B08) etc.\n",
    "        num_slices (int): The amount of slices to split the composite up into while building it.\n",
    "                          More slices will use less RAM\n",
    "    \"\"\"\n",
    "    # Loop through each band, getting a median estimate for each\n",
    "    crs = None\n",
    "    transform = None\n",
    "\n",
    "# Find each optical and cloud file in the input directory for the mgrs_tile and band combination\n",
    "\n",
    "    mgrs_str = f'{mgrs_tile[:2]}_{mgrs_tile[2]}_{mgrs_tile[3:]}'.upper()\n",
    "    band_str = band.upper()\n",
    "    dates = os.listdir(in_dir)\n",
    "    file_sets = []\n",
    "    for date_dir in dates:\n",
    "        optical_file = None\n",
    "        cloud_file = None\n",
    "        for file in os.listdir(os.path.join(in_dir, date_dir)):\n",
    "            f_up = file.upper()\n",
    "            file_path = os.path.join(in_dir, date_dir, file)\n",
    "            if mgrs_str in file and band_str in file:\n",
    "                optical_file = file_path\n",
    "            elif mgrs_str in file and 'MSK_CLOUDS_B00' in file:\n",
    "                cloud_file = file_path\n",
    "\n",
    "        if optical_file is None or cloud_file is None:\n",
    "            continue\n",
    "\n",
    "        file_sets.append((optical_file, cloud_file))  \n",
    "\n",
    " # Resolve the crs and other optical file attributes (loop through all until found)\n",
    "    crs = None\n",
    "    transform = None\n",
    "    g_nrows = None\n",
    "    g_ncols = None\n",
    "\n",
    "    for file_set in file_sets:\n",
    "        with rasterio.open(file_set[0], 'r', driver='JP2OpenJPEG') as rf:\n",
    "            g_nrows = rf.meta['height'] if g_nrows is None else g_nrows\n",
    "            g_ncols = rf.meta['width'] if g_ncols is None else g_ncols\n",
    "            crs = rf.crs if crs is None else crs\n",
    "            transform = rf.transform if transform is None else transform\n",
    "            if crs is not None and transform is not None and g_nrows is not None and g_ncols is not None:\n",
    "                break\n",
    "\n",
    "    if crs is None or transform is None or g_nrows is None or g_ncols is None:\n",
    "        raise LookupError(f'Could not determine the following projection attributes from the available '\n",
    "                          f'sentinel2 files in {in_dir}: \\n' \\\n",
    "                          f'{\"CRS\" if crs is None else \"\"} '\n",
    "                          f'{\"Transform\" if transform is None else \"\"} '\n",
    "                          f'{\"Number of rows\" if g_nrows is None else \"\"} '\n",
    "                          f'{\"Number of columns\" if g_ncols is None else \"\"}')\n",
    "\n",
    "    # Determine the slicing bounds to save memory as we process\n",
    "    slice_height = g_nrows / num_slices\n",
    "    slice_end_pts = [int(i) for i in np.arange(0, g_nrows + slice_height, slice_height)]\n",
    "    slice_bounds = [(slice_end_pts[i], slice_end_pts[i + 1] - 1) for i in range(num_slices - 1)]\n",
    "    slice_bounds.append((slice_end_pts[-2], slice_end_pts[-1]))  \n",
    "\n",
    "# Correct the images one slice at a time, and then combine the slices. \n",
    "# must create a out_dir to store the slices\n",
    "\n",
    "    slice_dir = os.path.join(os.path.dirname(out_file), 'slices')\n",
    "    os.makedirs(slice_dir, exist_ok=True)\n",
    "    for k, row_bound in tqdm(enumerate(slice_bounds), desc=f'band={band}', total=num_slices, position=2):\n",
    "        slice_file_path = os.path. join(slice_dir, f'{row_bound[0]}_{row_bound[1]}.tif')\n",
    "        cloud_correct_imgs = []\n",
    "        for file_set in tqdm(file_sets, desc=f'slice {k + 1}', leave=False, position=3):\n",
    "            # Get data from files\n",
    "            optical_file = file_set[0]\n",
    "            cloud_file = file_set[1]\n",
    "            pixels = get_img_from_file(optical_file, g_ncols, np.float32, row_bound)\n",
    "            cloud_channels = get_cloud_mask_from_file(cloud_file, crs, transform, (g_nrows, g_ncols), row_bound)\n",
    "            if cloud_channels is None:\n",
    "                continue\n",
    "            # add to list to do median filter later\n",
    "            cloud_correct_imgs.append(nan_clouds(pixels, cloud_channels))\n",
    "            del pixels\n",
    "        corrected_stack = np.vstack([img.ravel() for img in cloud_correct_imgs])\n",
    "        median_corrected = np.nanmedian(corrected_stack, axis=0, overwrite_input=True)\n",
    "        median_corrected = median_corrected.reshape(cloud_correct_imgs[0].shape)\n",
    "        with rasterio.open(slice_file_path, 'w', driver='GTiff', width=g_ncols, height=g_nrows,\n",
    "                           count=1, crs=crs, transform=transform, dtype=np.float32) as wf:\n",
    "            wf.write(median_corrected.astype(np.float32), 1)\n",
    "\n",
    "        # release mem\n",
    "        median_corrected = []\n",
    "        del median_corrected\n",
    "        corrected_stack = []\n",
    "        del corrected_stack\n",
    "\n",
    "# Combine slices\n",
    "    with rasterio.open(out_file, 'w', driver='GTiff', width=g_ncols, height=g_nrows,\n",
    "                       count=1, crs=crs, transform=transform, dtype=np.float32) as wf:\n",
    "        for slice_file in os.listdir(slice_dir):\n",
    "            bound_split = slice_file.split('.')[0].split('_')\n",
    "            top_bound = int(bound_split[0])\n",
    "            bottom_bound = int(bound_split[1])\n",
    "            with rasterio.open(os.path.join(slice_dir, slice_file), 'r', driver='GTiff') as rf:\n",
    "                wf.write(\n",
    "                    rf.read(1),\n",
    "                    window=Window.from_slices(\n",
    "                        slice(top_bound, bottom_bound),\n",
    "                        slice(0, g_ncols)\n",
    "                    ),\n",
    "                    indexes=1\n",
    "                )\n",
    "\n",
    "    shutil.rmtree(slice_dir)\n",
    "    print(f'Wrote file to {out_file}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13282828-d527-4d75-b052-a0fe0d071888",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The following codes by Lilly Jones. Alternate function to mask out clouds and water using the SCL (Scene Classification Layer) attribute of Sentinel-2 data.\n",
    "\n",
    "# Read the SCL band, apply a mask to identify cloud and water areas, and then apply this mask to the other bands. \n",
    "# Choice depends on the quality and accuracy of the SCL band in your data, otherwise use the cloud probability method (cloud mask) by Erick Verleye above.\n",
    "\n",
    "def mask_clouds_and_water(input_directory, output_directory):\n",
    "    \"\"\"\n",
    "    Masks out clouds and water in Sentinel-2 data using the SCL band.\n",
    "\n",
    "    Args:\n",
    "        input_directory (str): Directory containing Sentinel-2 data.\n",
    "        output_directory (str): Directory where masked data will be stored.\n",
    "\n",
    "    Returns:\n",
    "        None\n",
    "    \"\"\"\n",
    "\n",
    "    # List all Sentinel-2 band files in the input directory\n",
    "    band_files = glob.glob(os.path.join(input_directory, 'B*.jp2'))\n",
    "\n",
    "    # Create output directory\n",
    "    os.makedirs(output_directory, exist_ok=True)\n",
    "\n",
    "    # Determine the SCL file path (assuming it's in the same directory)\n",
    "    scl_file = os.path.join(input_directory, 'SCL.jp2')\n",
    "\n",
    "    for band_file in band_files:\n",
    "        # Define the output file path\n",
    "        output_file = os.path.join(output_directory, os.path.basename(band_file))\n",
    "\n",
    "        # Open the input band and SCL files\n",
    "        ds_band = gdal.Open(band_file, gdalconst.GA_ReadOnly)\n",
    "        ds_scl = gdal.Open(scl_file, gdalconst.GA_ReadOnly)\n",
    "\n",
    "        # Read band data as a numpy array\n",
    "        band_data = ds_band.ReadAsArray()\n",
    "\n",
    "        # Read SCL data as a numpy array\n",
    "        scl_data = ds_scl.ReadAsArray()\n",
    "\n",
    "        # Create a mask for cloudy and water pixels\n",
    "        # In the SCL band, value 3 corresponds to cloud and value 6 corresponds to water\n",
    "        cloud_water_mask = np.logical_or(scl_data == 3, scl_data == 6)\n",
    "\n",
    "        # Apply the mask to the band data by setting masked values to a nodata value (e.g., 0)\n",
    "        band_data[cloud_water_mask] = 0\n",
    "\n",
    "        # Create an output GeoTIFF file with the same metadata as the input band\n",
    "        driver = gdal.GetDriverByName(\"GTiff\")\n",
    "        output_ds = driver.Create(output_file, ds_band.RasterXSize, ds_band.RasterYSize, 1, gdalconst.GDT_Float32)\n",
    "        output_ds.SetGeoTransform(ds_band.GetGeoTransform())\n",
    "        output_ds.SetProjection(ds_band.GetProjection())\n",
    "        output_ds.GetRasterBand(1).WriteArray(band_data)\n",
    "\n",
    "        # Close the datasets\n",
    "        ds_band = None\n",
    "        ds_scl = None\n",
    "        output_ds = None\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    # Define input and output directories\n",
    "    input_directory = \"path_to_input_directory\"  # Replace with the directory containing Sentinel-2 data\n",
    "    output_directory = \"path_to_output_directory\"  # Replace with your desired output directory\n",
    "\n",
    "    # Call the mask_clouds_and_water function\n",
    "    mask_clouds_and_water(input_directory, output_directory)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3ae725a-a1a4-4ad4-9204-87d5c271b943",
   "metadata": {},
   "outputs": [],
   "source": [
    "# VSI cloud cover\n",
    "\n",
    "# Imports\n",
    "\n",
    "import gdal\n",
    "import os\n",
    "import numpy as np\n",
    "import rasterio\n",
    "\n",
    "# Open Sentinel scene using gdal\n",
    "# Set GDAL environment variables\n",
    "\n",
    "os.environ['CPL_TMPDIR'] = '/vsimem'  # Use /vsimem for in-memory VSI\n",
    "gdal.UseExceptions()\n",
    "\n",
    "# Input Sentinel-2 image URL\n",
    "input_image_path = '/vsicurl/https://your-sentinel-2-image-url'\n",
    "\n",
    "# Open the Sentinel-2 image using GDAL\n",
    "dataset = gdal.Open(input_image_path)\n",
    "\n",
    "# Define a cloud detection function\n",
    "def cloud_detection(band_red, band_swir, threshold=0.2):\n",
    "    # Calculate a cloud probability map using custom logic\n",
    "    cloud_prob = band_swir / (band_red + band_swir + 1e-10)\n",
    "\n",
    "    # Set cloudy pixels to 1 and non-cloudy pixels to 0\n",
    "    cloud_mask = cloud_prob > threshold\n",
    "\n",
    "    return cloud_mask\n",
    "\n",
    "# Read red and SWIR bands using rasterio\n",
    "with rasterio.open(input_image_path) as src:\n",
    "    red_band = src.read(3)\n",
    "    swir_band = src.read(11)\n",
    "\n",
    "# Simple algorithm for cloud detection using a threshold\n",
    "cloud_mask = cloud_detection(red_band, swir_band)\n",
    "\n",
    "# Apply the cloud mask to the scene (sets the pixel values to 0 for cloudy areas based on the cloud mask)\n",
    "\n",
    "for i in range(1, dataset.RasterCount + 1):\n",
    "    band = dataset.GetRasterBand(i)\n",
    "    data = band.ReadAsArray()\n",
    "    data[cloud_mask] = 0  # Set cloudy pixels to 0\n",
    "    band.WriteArray(data)\n",
    "\n",
    "# Save the output\n",
    "\n",
    "output_image_path = '/vsimem/cloud_corrected_sentinel2.tif'  # Output VSI path\n",
    "dataset.FlushCache()\n",
    "dataset = None\n",
    "\n",
    "# Create a copy of the corrected image with the VSI path\n",
    "\n",
    "gdal.Translate(output_image_path, input_image_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f239039e-655a-4330-9a19-292f909fc154",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to extract spectra from Sentinel-2 data and resample relevant bands\n",
    "\n",
    "def extract_spectra_and_resample(xml_file, geotiff_file, output_csv_file):\n",
    "    try:\n",
    "        # Parse the Sentinel-2 XML file to extract band information\n",
    "        tree = ET.parse(xml_file)\n",
    "        root = tree.getroot()\n",
    "\n",
    "        bands_info = []\n",
    "        for elem in root.iterfind(\".//BANDS/BAND\"):\n",
    "            band_name = elem.find(\"BAND_NAME\").text\n",
    "            band_center_wavelength = float(elem.find(\"CENTRAL_WAVELENGTH\").text)\n",
    "            bands_info.append((band_name, band_center_wavelength))\n",
    "\n",
    "        # Open the geotiff file using rasterio\n",
    "        with rasterio.open(geotiff_file) as src:\n",
    "            band_count = src.count\n",
    "            band_data = []\n",
    "\n",
    "            # Iterate over each band and read data\n",
    "            for band_index in range(1, band_count + 1):\n",
    "                band_data.append(src.read(band_index))\n",
    "\n",
    "        # Resample the 20m resolution bands to 10m resolution using the average method\n",
    "        resampled_band_data = []\n",
    "        resample_indices = [1, 2, 3, 4, 8, 11, 12]  # Indices of 20m bands to be resampled\n",
    "        for i in range(len(band_data)):\n",
    "            if i in resample_indices:\n",
    "                resampled_band = (band_data[i] + band_data[i + 1]) / 2.0\n",
    "                resampled_band_data.append(resampled_band)\n",
    "            elif i == 0 or i == 5:  # Bands 2 and 8 are already 10m resolution\n",
    "                resampled_band_data.append(band_data[i])\n",
    "\n",
    "        # Extract reflectance data from each band and create a dataframe\n",
    "        spectra_data = {band_name: resampled_band_data[i] for i, (band_name, _) in enumerate(bands_info)}\n",
    "        df = pd.DataFrame(spectra_data)\n",
    "\n",
    "        # Add the center wavelength information to the dataframe\n",
    "        df['Wavelength (micrometers)'] = [wavelength for (_, wavelength) in bands_info]\n",
    "\n",
    "        # Rename columns for samples\n",
    "        df.columns = [f'Sample_{i}' for i in range(df.shape[1] - 1)] + ['Wavelength (micrometers)']\n",
    "\n",
    "        # Display the dataframe (Optional: comment out if you have many samples)\n",
    "        print(df)\n",
    "\n",
    "        # Save the dataframe to a CSV file\n",
    "        df.to_csv(output_csv_file, index=False)\n",
    "\n",
    "    except Exception as e:\n",
    "        print(f\"Error: {e}\")\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    sentinel2_xml_file = \"path/to/sentinel2.xml\"\n",
    "    sentinel2_geotiff_file = \"path/to/sentinel2_geotiff.tif\"\n",
    "    output_csv_file = \"spectra_table.csv\"\n",
    "    \n",
    "    # Extract spectra and resample relevant bands\n",
    "    extract_spectra_and_resample(sentinel2_xml_file, sentinel2_geotiff_file, output_csv_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bacf10e8-a4c4-4961-9f48-9373b557f799",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the Normalized Difference Vegetation Index (NDVI) from Sentinel-2 bands 4 (red) and 8 (nir). NDVI is calculated as (NIR - Red) / (NIR + Red).\n",
    "\n",
    "def calculate_ndvi(input_directory, output_directory):\n",
    "    \"\"\"\n",
    "    Calculates NDVI from Sentinel-2 bands 4 (red) and 8 (NIR).\n",
    "\n",
    "    Args:\n",
    "        input_directory (str): Directory containing Sentinel-2 data.\n",
    "        output_directory (str): Directory where NDVI data will be stored.\n",
    "\n",
    "    Returns:\n",
    "        None\n",
    "    \"\"\"\n",
    "\n",
    "    # List all Sentinel-2 band files in the input directory\n",
    "    band4_file = glob.glob(os.path.join(input_directory, 'B04.jp2'))[0]  # Red band\n",
    "    band8_file = glob.glob(os.path.join(input_directory, 'B08.jp2'))[0]  # NIR band\n",
    "\n",
    "    # Create output directory\n",
    "    os.makedirs(output_directory, exist_ok=True)\n",
    "\n",
    "    # Initialize variables to store band data\n",
    "    red_band = None\n",
    "    nir_band = None\n",
    "\n",
    "    for band_file in [band4_file, band8_file]:\n",
    "        # Define the output file path\n",
    "        output_file = os.path.join(output_directory, os.path.basename(band_file).replace('B0', 'NDVI'))\n",
    "\n",
    "        # Open the input band file\n",
    "        ds_band = gdal.Open(band_file, gdalconst.GA_ReadOnly)\n",
    "\n",
    "        # Read band data as a numpy array\n",
    "        band_data = ds_band.ReadAsArray()\n",
    "\n",
    "        # Assign band data to respective variables\n",
    "        if 'B04' in band_file:\n",
    "            red_band = band_data.astype(np.float32)\n",
    "        elif 'B08' in band_file:\n",
    "            nir_band = band_data.astype(np.float32)\n",
    "\n",
    "        # Close the dataset\n",
    "        ds_band = None\n",
    "\n",
    "    # Check if both bands are available\n",
    "    if red_band is not None and nir_band is not None:\n",
    "        # Calculate NDVI\n",
    "        ndvi = (nir_band - red_band) / (nir_band + red_band)\n",
    "\n",
    "        # Create an output GeoTIFF file with the same metadata as the input band\n",
    "        driver = gdal.GetDriverByName(\"GTiff\")\n",
    "        output_ds = driver.Create(output_file, ds_band.RasterXSize, ds_band.RasterYSize, 1, gdalconst.GDT_Float32)\n",
    "        output_ds.SetGeoTransform(ds_band.GetGeoTransform())\n",
    "        output_ds.SetProjection(ds_band.GetProjection())\n",
    "        output_ds.GetRasterBand(1).WriteArray(ndvi)\n",
    "\n",
    "        # Close the output dataset\n",
    "        output_ds = None\n",
    "    else:\n",
    "        print(\"Error: Both red and NIR bands are required to calculate NDVI.\")\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    # Define input and output directories\n",
    "    input_directory = \"path_to_input_directory\"  # Replace with the directory containing Sentinel-2 data\n",
    "    output_directory = \"path_to_output_directory\"  # Replace with your desired output directory\n",
    "\n",
    "    # Call the calculate_ndvi function\n",
    "    calculate_ndvi(input_directory, output_directory)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c6af4c1-5e29-4b67-a477-bdfaf5197d1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate mean NDVI for all time-steps\n",
    "\n",
    "\n",
    "def calculate_ndvi(input_directory, output_directory):\n",
    "    \"\"\"\n",
    "    Calculates NDVI from Sentinel-2 bands 4 (red) and 8 (NIR).\n",
    "\n",
    "    Args:\n",
    "        input_directory (str): Directory containing Sentinel-2 data for multiple time-steps.\n",
    "        output_directory (str): Directory where mean NDVI data will be stored.\n",
    "\n",
    "    Returns:\n",
    "        None\n",
    "    \"\"\"\n",
    "    # List all time-steps (subdirectories) in the input directory\n",
    "    time_steps = [d for d in os.listdir(input_directory) if os.path.isdir(os.path.join(input_directory, d))]\n",
    "\n",
    "    # Create the output directory if it doesn't exist\n",
    "    os.makedirs(output_directory, exist_ok=True)\n",
    "\n",
    "    # Initialize arrays to store NDVI values and count of valid pixels\n",
    "    ndvi_sum = None\n",
    "    valid_pixel_count = None\n",
    "\n",
    "    for time_step in time_steps:\n",
    "        # Get the directory containing Sentinel-2 data for the current time-step\n",
    "        time_step_directory = os.path.join(input_directory, time_step)\n",
    "\n",
    "        # List all Sentinel-2 band files in the current time-step directory\n",
    "        band4_file = glob.glob(os.path.join(time_step_directory, 'B04.jp2'))[0]  # Red band\n",
    "        band8_file = glob.glob(os.path.join(time_step_directory, 'B08.jp2'))[0]  # NIR band\n",
    "\n",
    "        for band_file in [band4_file, band8_file]:\n",
    "            # Open the input band file\n",
    "            ds_band = gdal.Open(band_file, gdalconst.GA_ReadOnly)\n",
    "\n",
    "            # Read band data as a numpy array\n",
    "            band_data = ds_band.ReadAsArray()\n",
    "\n",
    "            # Calculate NDVI\n",
    "            if 'B04' in band_file:\n",
    "                red_band = band_data.astype(np.float32)\n",
    "            elif 'B08' in band_file:\n",
    "                nir_band = band_data.astype(np.float32)\n",
    "\n",
    "        # Calculate NDVI for the current time-step\n",
    "        ndvi = (nir_band - red_band) / (nir_band + red_band)\n",
    "\n",
    "        # Initialize the arrays for NDVI sum and valid pixel count\n",
    "        if ndvi_sum is None:\n",
    "            ndvi_sum = np.zeros_like(ndvi, dtype=np.float64)\n",
    "            valid_pixel_count = np.zeros_like(ndvi, dtype=np.int64)\n",
    "\n",
    "        # Update the sum of NDVI values and valid pixel count\n",
    "        ndvi_sum += ndvi\n",
    "        valid_pixel_count += 1\n",
    "\n",
    "    # Calculate mean NDVI by dividing the sum by the valid pixel count\n",
    "    mean_ndvi = np.divide(ndvi_sum, valid_pixel_count, where=valid_pixel_count != 0)\n",
    "\n",
    "    # Get the metadata (geotransform and projection) from one of the input bands\n",
    "    ds_band = gdal.Open(band4_file, gdalconst.GA_ReadOnly)\n",
    "    geotransform = ds_band.GetGeoTransform()\n",
    "    projection = ds_band.GetProjection()\n",
    "    ds_band = None\n",
    "\n",
    "    # Create an output GeoTIFF file with the same metadata as the input band\n",
    "    output_file = os.path.join(output_directory, 'mean_ndvi.tif')\n",
    "    driver = gdal.GetDriverByName(\"GTiff\")\n",
    "    output_ds = driver.Create(output_file, mean_ndvi.shape[1], mean_ndvi.shape[0], 1, gdalconst.GDT_Float32)\n",
    "    output_ds.SetGeoTransform(geotransform)\n",
    "    output_ds.SetProjection(projection)\n",
    "    output_ds.GetRasterBand(1).WriteArray(mean_ndvi)\n",
    "\n",
    "    # Close the output dataset\n",
    "    output_ds = None\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    # Define input and output directories\n",
    "    input_directory = \"path_to_input_directory\"  # Replace with the directory containing Sentinel-2 data for multiple time-steps\n",
    "    output_directory = \"path_to_output_directory\"  # Replace with your desired output directory\n",
    "\n",
    "    # Call the calculate_ndvi function\n",
    "    calculate_ndvi(input_directory, output_directory)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8900e4ed-5fa7-43a7-ba9c-f00f6d3f4cea",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# plot mean NDVI\n",
    "\n",
    "\n",
    "def plot_mean_ndvi(ndvi_file):\n",
    "    \"\"\"\n",
    "    Plots the mean NDVI from a GeoTIFF file.\n",
    "\n",
    "    Args:\n",
    "        ndvi_file (str): Path to the GeoTIFF file containing mean NDVI data.\n",
    "\n",
    "    Returns:\n",
    "        None\n",
    "    \"\"\"\n",
    "    # Open the GeoTIFF file using rasterio\n",
    "    with rasterio.open(ndvi_file) as src:\n",
    "        # Read the NDVI data as a numpy array\n",
    "        ndvi_data = src.read(1)\n",
    "\n",
    "        # Create a colormap for NDVI (green to red)\n",
    "        cmap = plt.cm.RdYlGn\n",
    "\n",
    "        # Set nodata values to transparent\n",
    "        cmap.set_bad(alpha=0)\n",
    "\n",
    "        # Create a plot\n",
    "        plt.figure(figsize=(10, 10))\n",
    "        plt.imshow(ndvi_data, cmap=cmap, vmin=-1, vmax=1)\n",
    "        plt.colorbar(label='NDVI')\n",
    "        plt.title('Mean NDVI')\n",
    "        plt.axis('off')\n",
    "\n",
    "        # Show the plot\n",
    "        plt.show()\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    # Define the path to the mean NDVI GeoTIFF file\n",
    "    ndvi_file = \"path_to_mean_ndvi.tif\"  # Replace with the actual file path\n",
    "\n",
    "    # Call the plot_mean_ndvi function\n",
    "    plot_mean_ndvi(ndvi_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b408655-16da-482a-afcc-04e1df0b1b03",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot spectral signature for all bands reads and plots the pixel values for each band\n",
    "\n",
    "\n",
    "def plot_spectral_signature(input_directory):\n",
    "    \"\"\"\n",
    "    Plots the spectral signature from all bands of Sentinel-2 data.\n",
    "\n",
    "    Args:\n",
    "        input_directory (str): Directory containing Sentinel-2 data.\n",
    "\n",
    "    Returns:\n",
    "        None\n",
    "    \"\"\"\n",
    "    # List all Sentinel-2 band files in the input directory\n",
    "    band_files = sorted(glob.glob(os.path.join(input_directory, 'B*.jp2')))\n",
    "\n",
    "    # Initialize an empty list to store band names\n",
    "    band_names = []\n",
    "\n",
    "    # Initialize empty arrays to store spectral data for each band\n",
    "    spectral_data = []\n",
    "\n",
    "    # Open and read data for each band\n",
    "    for band_file in band_files:\n",
    "        # Get the band name from the file name (e.g., B02.jp2)\n",
    "        band_name = os.path.splitext(os.path.basename(band_file))[0]\n",
    "        band_names.append(band_name)\n",
    "\n",
    "        # Open the band file using GDAL\n",
    "        ds_band = gdal.Open(band_file, gdal.GA_ReadOnly)\n",
    "\n",
    "        # Read the band data as a numpy array\n",
    "        band_data = ds_band.ReadAsArray()\n",
    "\n",
    "        # Append the band data to the spectral_data list\n",
    "        spectral_data.append(band_data)\n",
    "\n",
    "    # Create a plot for the spectral signature\n",
    "    plt.figure(figsize=(12, 6))\n",
    "    for band_name, band_data in zip(band_names, spectral_data):\n",
    "        plt.plot(band_data, label=band_name)\n",
    "\n",
    "    # Set plot labels and title\n",
    "    plt.xlabel('Pixel Value')\n",
    "    plt.ylabel('Reflectance')\n",
    "    plt.title('Spectral Signature of Sentinel-2 Bands')\n",
    "    plt.legend()\n",
    "    plt.grid(True)\n",
    "\n",
    "    # Show the plot\n",
    "    plt.show()\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    # Define input directory containing Sentinel-2 data\n",
    "    input_directory = \"path_to_input_directory\"  # Replace with the directory containing Sentinel-2 data\n",
    "\n",
    "    # Call the plot_spectral_signature function\n",
    "    plot_spectral_signature(input_directory)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c02959eb-ccaf-427a-9aa9-4b27e276deda",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create spectral libraries\n",
    "# Collect spectra across seasons and years (> 50 per class)\n",
    "# Classes: each vegetation type's green vegetation and dry vegetation, bare soil, water and shade\n",
    "# Collect spectra at the pixel level (capture the variability within classes, sunlight pixels preferentially - we will have the shade class)\n",
    "# NEON R code to be translated to Python: https://github.com/earthlab/neonhs\n",
    "# May use Python functions from this: https://www.spectralpython.net/\n",
    "\n",
    "# 4. Create metadata\n",
    "# A final, merged spectral library must have a metadata (.csv) with at least spectrum ID (column 1), class (2), and sub-class (3) description. \n",
    "# For the vegetation type classes the order is: 1. class: green vegetation or non-photosynthetic vegetation (i.e. dry matter); 2. sub-class: vegetation type. \n",
    "# For the other classes (bare # # soil, water..), repeat the class content in the sub-class column.\n",
    "\n",
    "# 5. Endmember Selection - WIP\n",
    "\n",
    "# 6. Multiple Endmber Spectral Mixture Analysis Emulation WIP\n",
    "\n",
    "# 7. Dry matter Index Calculation WIP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c78751dc-2e24-456c-9319-0112c36dc0cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# export output to PetaLibrary \n",
    "\n",
    "#!/bin/bash\n",
    "\n",
    "# Variables\n",
    "local_file=\"/path/to/local/file.zip\"\n",
    "remote_username=\"your_username\"\n",
    "remote_host=\"petalibrary.colorado.edu\"\n",
    "remote_directory=\"/path/to/remote/directory/\"\n",
    "\n",
    "# Copy file using SCP\n",
    "scp \"$local_file\" \"$remote_username@$remote_host:$remote_directory\"\n",
    "\n",
    "echo \"File transferred successfully.\"\n",
    "\n",
    "bash transfer_script.sh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de1aa568-06c4-4169-b74a-573a16fd8a97",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Attempts to use vsicurl and STAC\n",
    "\n",
    "import geopandas as gpd\n",
    "import osmdata\n",
    "from rasterio.plot import show\n",
    "from rasterio.windows import from_bounds\n",
    "import boto3\n",
    "from io import BytesIO\n",
    "import subprocess\n",
    "import json\n",
    "from zipfile import ZipFile\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# AWS credentials and bucket name\n",
    "aws_access_key_id = 'YOUR_ACCESS_KEY_ID'\n",
    "aws_secret_access_key = 'YOUR_SECRET_ACCESS_KEY'\n",
    "bucket_name = 'your-bucket-name'\n",
    "\n",
    "# Create an S3 client\n",
    "s3 = boto3.client('s3', aws_access_key_id=aws_access_key_id, aws_secret_access_key=aws_secret_access_key)\n",
    "\n",
    "# Load and filter shapefiles\n",
    "world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))\n",
    "usa = world[world['name'] == 'United States']\n",
    "california = usa[usa['name'] == 'California']  # Filter for California\n",
    "\n",
    "# Read and filter GMW shapefiles\n",
    "base_url = \"https://datadownload-production.s3.amazonaws.com\"\n",
    "y2017 = \"GMW_v3_2017.zip\"\n",
    "y2020 = \"GMW_v3_2020.zip\"\n",
    "\n",
    "gmw2017_california = gpd.read_file(f\"{base_url}/{y2017}\", vfs=\"zip://\")\n",
    "gmw2017_california = gmw2017_california[gmw2017_california.intersects(california.unary_union)]\n",
    "\n",
    "gmw2020_california = gpd.read_file(f\"{base_url}/{y2020}\", vfs=\"zip://\")\n",
    "gmw2020_california = gmw2020_california[gmw2020_california.intersects(california.unary_union)]\n",
    "\n",
    "# Plot the filtered shapefiles\n",
    "gmw2017_california.plot()\n",
    "gmw2020_california.plot()\n",
    "\n",
    "# Define a bounding box for California and transform it to different coordinate systems\n",
    "bbox = osmdata.getbb(\"California, USA\")\n",
    "bbox_4326 = bbox.to_crs(epsg=4326)\n",
    "bbox_32618 = bbox.to_crs(epsg=32618)\n",
    "\n",
    "def lambda_handler(event, context):\n",
    "    # Replace 'YOUR_SENTINEL_API_KEY' with your actual Sentinel Hub API key\n",
    "    sentinel_api_key = 'YOUR_SENTINEL_API_KEY'\n",
    "    \n",
    "    # Get the input parameters from the Step Function event\n",
    "    bbox = event['bbox']\n",
    "    time_range = event['time_range']\n",
    "    \n",
    "    # Call Sentinel Hub API to download Sentinel-2 data\n",
    "    download_url = f'https://services.sentinel-hub.com/api/v1/process/download'\n",
    "    request_payload = {\n",
    "        'input': {\n",
    "            'bounds': {\n",
    "                'geometry': {\n",
    "                    'type': 'Polygon',\n",
    "                    'coordinates': [bbox]\n",
    "                }\n",
    "            },\n",
    "            'data': [{\n",
    "                'type': 'S2L1C',\n",
    "                'dataFilter': {\n",
    "                    'timeRange': time_range,\n",
    "                    'maxCloudCoverage': 10\n",
    "                }\n",
    "            }]\n",
    "        },\n",
    "        'output': {\n",
    "            'width': 512,\n",
    "            'height': 512,\n",
    "            'responses': [{\n",
    "                'identifier': 'default',\n",
    "                'format': {\n",
    "                    'type': 'image/tiff'\n",
    "                }\n",
    "            }]\n",
    "        }\n",
    "    }\n",
    "    \n",
    "    response = requests.post(download_url, json=request_payload, headers={'Authorization': f'Bearer {sentinel_api_key}'})\n",
    "    \n",
    "    # Process the downloaded data \n",
    "    with ZipFile(BytesIO(response.content)) as zip_file:\n",
    "        zip_file.extractall('/tmp')  # Extract the contents to /tmp directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a7c9794-6157-4d3a-b570-1b7ca372318f",
   "metadata": {},
   "outputs": [],
   "source": [
    "    \n",
    "# VSI processing steps\n",
    "\n",
    "import subprocess\n",
    "import json\n",
    "from io import BytesIO\n",
    "from zipfile import ZipFile\n",
    "\n",
    "def lambda_handler(event, context):\n",
    "    # Replace 'YOUR_STAC_API_URL' with the actual STAC API endpoint URL\n",
    "    stac_api_url = 'YOUR_STAC_API_URL'\n",
    "    \n",
    "    # Get the input parameters from the Step Function event\n",
    "    # For example, you might pass the bounding box coordinates, time range, etc.\n",
    "    bbox = event['bbox']\n",
    "    time_range = event['time_range']\n",
    "    \n",
    "    # Use STAC API to search for Sentinel-2 data\n",
    "    search_payload = {\n",
    "        'bbox': bbox,\n",
    "        'time': time_range,\n",
    "        'collections': ['sentinel-2-l1c']\n",
    "    }\n",
    "    \n",
    "    search_url = f'{stac_api_url}/search'\n",
    "    search_response = subprocess.run(['vsi', 'curl', '-X', 'POST', '-d', json.dumps(search_payload), search_url], capture_output=True, text=True)\n",
    "    \n",
    "    # Parse the search response to get the download URL\n",
    "    download_url = json.loads(search_response.stdout)['features'][0]['assets']['data']['href']\n",
    "    \n",
    "    # Download Sentinel-2 data using vsi curl\n",
    "    download_response = subprocess.run(['vsi', 'curl', download_url], capture_output=True, text=True)\n",
    "    \n",
    "    # Process the downloaded data (you might want to use a more sophisticated method)\n",
    "    with ZipFile(BytesIO(download_response.stdout.encode())) as zip_file:\n",
    "        zip_file.extractall('/tmp')  # Extract the contents to /tmp directory\n",
    "    \n",
    "    # Your VSI processing steps here (replace with actual processing logic)\n",
    "    processed_data = b'Your processed data'\n",
    "    \n",
    "    return {\n",
    "        'statusCode': 200,\n",
    "        'body': processed_data\n",
    "    }\n",
    "   \n",
    "\n",
    "# Step Function definition\n",
    "{\n",
    "  \"Comment\": \"Sentinel-2 Data Processing\",\n",
    "  \"StartAt\": \"DownloadAndProcess\",\n",
    "  \"States\": {\n",
    "    \"DownloadAndProcess\": {\n",
    "      \"Type\": \"Task\",\n",
    "      \"Resource\": \"arn:aws:lambda:REGION:ACCOUNT_ID:function:YourLambdaFunctionName\",\n",
    "      \"End\": true\n",
    "    }\n",
    "  }\n",
    "}\n",
    "\n",
    "# Save the processed data back to S3\n",
    "processed_data = b'Your processed data in the cloud'  # Replace this with your actual processed data\n",
    "with BytesIO(processed_data) as processed_data_buffer:\n",
    "    s3.upload_fileobj(processed_data_buffer, bucket_name, 'path/to/processed_data.tif')\n",
    "\n",
    "# Data Visualization \n",
    "data = pd.read_csv('path/to/processed_data.csv') \n",
    "plt.figure(figsize=(10, 5))\n",
    "plt.plot(data['wavelength_nm'], data['reflectance'])\n",
    "plt.xlabel('Wavelength (nm)')\n",
    "plt.ylabel('Reflectance')\n",
    "plt.title('Reflectance vs. Wavelength')\n",
    "plt.savefig('reflectance_plot.png')\n",
    "\n",
    "# Upload the visualization to S3\n",
    "with open('reflectance_plot.png', 'rb') as plot_file:\n",
    "    s3.upload_fileobj(plot_file, bucket_name, 'path/to/reflectance_plot.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1613db7a-af08-4c50-b96c-e762969a350e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# use output shapefile in ArcGIS with arcpy\n",
    "\n",
    "import arcpy\n",
    "\n",
    "# Set the workspace where you want to save the output\n",
    "arcpy.env.workspace = r'C:\\Path\\To\\Output\\Directory'\n",
    "\n",
    "# Create a feature class to save your data\n",
    "output_feature_class = \"OutputData.shp\"  # Change the name as needed\n",
    "\n",
    "# Create a new feature class (shapefile)\n",
    "arcpy.CreateFeatureclass_management(arcpy.env.workspace, output_feature_class, \"POINT\")\n",
    "\n",
    "# Open an insert cursor to add points to the feature class\n",
    "cursor = arcpy.da.InsertCursor(output_feature_class, [\"SHAPE@XY\"])\n",
    "\n",
    "# Loop through your data and add points to the feature class\n",
    "for point in data:\n",
    "    cursor.insertRow([point])\n",
    "\n",
    "# Clean up\n",
    "del cursor\n",
    "\n",
    "print(f\"Data has been saved to {output_feature_class}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:SCE]",
   "language": "python",
   "name": "conda-env-SCE-py"
  },
  "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.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}