GPU_and_CPU_GEE-Data_Cube_Generator-F2A.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Li_TRYN9YJtq"
   },
   "source": [
    "# GEE - Data Cube Generator\n",
    "\n",
    "## Introduction\n",
    "\n",
    "This Colab notebook is designed to process and analyze geospatial data using Google Earth Engine (GEE). It focuses on deriving topographic, climatic, and spectral information for a specified region of interest (ROI). By leveraging GEE's powerful cloud computing capabilities, this workflow enables efficient data handling and export of key environmental variables for further analysis.\n",
    "\n",
    "\n",
    "### What Does This Notebook Do?:\n",
    "\n",
    "1) Loads a Region of Interest (ROI):\n",
    "\n",
    "- Accepts shapefiles or GeoPackages to define the area for analysis.\n",
    "- Applies spatial buffers to adapt the region for specific datasets.\n",
    "\n",
    "2) Processes Geospatial Variables:\n",
    "\n",
    "- Topographic Variables: Elevation, slope, aspect, landforms, and topographic position index (TPI).\n",
    "- Climate Variables: Total precipitation and average temperature over a defined period.\n",
    "- Spectral Indices: NDVI, NDBI, NDWI, and EVI derived from Sentinel-2 imagery.\n",
    "- Sentinel-2 Variables: Red, Green, Blue, R1, R2, R3, NIR, R4, SWIR1, and SWIR2 bands with cloud correction included.\n",
    "- Radar Data: VV and VH polarizations from Sentinel-1 with slope correction applied.\n",
    "- EVI Statistics: Summarizes vegetation dynamics for selected time periods.\n",
    "\n",
    "3) Defines and Analyzes Temporal Periods:\n",
    "\n",
    "- Period 1 (P1): 15 March (desired_year) - 30 May (desired_year).\n",
    "- Period 2 (P2): 01 June (desired_year) - 15 August (desired_year).\n",
    "- Period 3 (P3): 15 November (previous_year) - 15 February (desired_year).\n",
    "\n",
    "4) Exports Variables to Google Drive as a pickle file:\n",
    "\n",
    "- Variables are stored in Google Drive as a pickle file to change from GPU to CPU (GPU is needed for the preprocessing but when exporting the Data it's not being used, so if it's not change the free usage of GPU that Google Colab offer will be used and you will not be able to use after for a period of time )\n",
    "\n",
    "5) Exports Results to Google Drive:\n",
    "\n",
    "- The pickle file is loaded and the variables are organized into batches and exported as GeoTIFF files for easy integration with GIS tools like QGIS or ArcGIS.\n",
    "\n",
    "\n",
    "### Prerequisites:\n",
    "\n",
    "To use this notebook, you will need:\n",
    "\n",
    "- A Google Earth Engine account.\n",
    "- Basic familiarity with Python and geospatial concepts.\n",
    "- A GeoPackage (.gpkg) or Shapefile (.shp) defining your region of interest.\n",
    "\n",
    "## Instructions and considerations\n",
    "\n",
    "This are the steps to follow for a correct function of the Google Colab:\n",
    "\n",
    "### Part 1: Prepare the Data\n",
    "- Go to the [Google Earth Engine web site](https://code.earthengine.google.com/) and sing in to create a proyect.\n",
    "- Make sure that the GPU is selected, go to **\"Runtime\"** on the menu above and go to **\"Change runtime type\"** to check it.\n",
    "- Inside the cell **\"Initial Setup/Set-up and User Inputs\"** make the next changes:\n",
    "  - In the cell **\"/Initializing Google Earth Engine\"** insert your GEE proyect id.\n",
    "  - Execute the cell **Mounting Google Drive** and accept the terms\n",
    "  - Select the paths of the region file (.shp or .gpkg) stored on your drive by using the panel on the right by selecting the folder icon and then right-click the desired file and copy the path and paste it on the variable called **uploaded_file**. (If the region is too complex is recommended to simplify the region before).\n",
    "  - Insert the year that you want to extract the images from on the variable **season_year**.\n",
    "  \n",
    "   (Ex: If I select the year \"2021\" the images obtained will be for the 2021-2022 Fire Season)\n",
    "  - Define the name of the folder where the output will be stored in the variable **\"drive_folder\"**.\n",
    "  - Define de CRS which the tiff files will be downloaded in the variable **\"export_crs\"** (by default is set to EPSG:32718).\n",
    "- After doing the changes above you can execute the cell **Initial Setup** (accept the terms to use GEE and Drive) and then execute the cell **Preparing Data**.\n",
    "\n",
    "  (**WARNING: Do not execute the rest of the code yet, only that cell**)\n",
    "\n",
    "\n",
    "### Part 2: Export the Data\n",
    "- Change the **Runtime** to CPU.\n",
    "- Execute the cell Initial Setup (make sure that the variables are the same that you used in **Preparing Data** and also accept the terms to use GEE and Drive)\n",
    "- After this you can execute the cell **Exporting Data**\n",
    "\n",
    "### **Other Considerations**\n",
    "- The reason why it's important to change from GPU to CPU when starting to extract the Data it's because if you don't do it a message will appear that the GPU it's selected but it's not being used. So if you don't change and wait until the files are downloaded you'll consume the GPU usage that Google Colab gives to the free acounts. So to prevent this the code is divided so you can change it.\n",
    "- Be careful to not execute unwanted cells when trying to expand. Do not press the \"play\" button of the cells until it's required. To expand a cell click on the \">\" at the left of the title of the cell.\n",
    "- All the varaibles are processed in EPSG:4326 and then downloaded in the desired CRS (NASADEM variables can't be directly downloaded with some CRS so we donwload them in EPSG:4326 and then it's changed)\n",
    "- While a cell is executing you can press the execute button of other cells and they will be executed after the previous cell is done."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "NfAVArbkh3_r"
   },
   "source": [
    "# Intial Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Zb7mMmCeUcdq"
   },
   "source": [
    "## Installing and Importing Required Libraries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "kXjQGmeGfLb6"
   },
   "source": [
    "### Install Libraries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "0CKId_FreLe_"
   },
   "source": [
    "To begin, we need to install some Python libraries that are not pre-installed in Google Colab but are necessary for this project:\n",
    "\n",
    "- `xarray`: For handling multi-dimensional arrays.\n",
    "- `rioxarray`: For geospatial raster operations.\n",
    "- `rasterio`: For reading and writing geospatial raster data.\n",
    "- `numpy`: For numerical operations.\n",
    "- `matplotlib`: For visualizing data.\n",
    "- `pyproj`: For handling cartographic projections and CRS transformations.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "E5SgTIT8dUYN"
   },
   "outputs": [],
   "source": [
    "!pip install xarray rioxarray rasterio numpy matplotlib pyproj"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "VAnQN4azfPt5"
   },
   "source": [
    "### Import Libraries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WOEzG8WWeBRT"
   },
   "source": [
    "This section imports the following libraries:\n",
    "- `ee`: The Google Earth Engine (GEE) API for managing geospatial data.\n",
    "- `geopandas`: For working with vector geospatial data.\n",
    "- `google.colab.drive`: To access files stored in your Google Drive.\n",
    "- `os`: For directory and file path management.\n",
    "- `time`: For time measuring when downloading files from GEE.\n",
    "- `pickle`: For storing variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "nL9J6owSeESs"
   },
   "outputs": [],
   "source": [
    "import ee\n",
    "import geopandas as gpd\n",
    "from google.colab import drive\n",
    "import os\n",
    "import time\n",
    "import pickle\n",
    "import rasterio\n",
    "from rasterio.warp import calculate_default_transform, reproject, Resampling\n",
    "from rasterio.transform import Affine\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "hB4oty-qU-wf"
   },
   "source": [
    "## Setup and User Inputs\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "pLjQci1OeSVj"
   },
   "source": [
    "### Initializing Google Earth Engine (GEE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "1Nm9logeS-7-"
   },
   "source": [
    "How to Set Up a GEE Project\n",
    "- Visit the Google Earth Engine website.\n",
    "- Click on \"Sign In\" and use your Google account to authenticate.\n",
    "- Go to the \"Projects\" section in the Console.\n",
    "- Create a new project and note down the Project ID.\n",
    "- Replace 'your_project_id' in the ee.Initialize() call with your actual Project ID.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ALOU0_T27p5I"
   },
   "outputs": [],
   "source": [
    "# Authenticate GEE\n",
    "ee.Authenticate()\n",
    "\n",
    "# Initialize GEE with the user's project\n",
    "# IMPORTANT: Replace 'your_project_id' with your actual project ID from GEE\n",
    "ee.Initialize(project='ee-your_project_id')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "yM6M0xq3ekSK"
   },
   "source": [
    "### Mounting Google Drive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "qceIJKTBehM_"
   },
   "outputs": [],
   "source": [
    "# Mount Google Drive with your Google account to access files\n",
    "drive.mount('/content/drive')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "C9UlPJRtf_zR"
   },
   "source": [
    "### User Input: Define Paths and desired year to export the files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "PZKJdZ5UTSMH"
   },
   "source": [
    "How to Reference Your Shapefile/GeoPackage after you mount Google Drive\n",
    "- Open the left panel in Colab and click the \"Files\" tab.\n",
    "- Upload your .gpkg or .shp file to a folder in your Google Drive.\n",
    "- Copy the full path of your file (e.g., /content/drive/MyDrive/YourFolder/your_file.gpkg).\n",
    "- Enter this path in the uploaded_file variable at the start of the code.\n",
    "- (If you need to simplify the region do it before importing the file)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "-EgQJ5Rj7xCF"
   },
   "outputs": [],
   "source": [
    "# Prompt the user to provide the path of their shapefile (.shp) or GeoPackage (.gpkg)\n",
    "# Example: \"/content/drive/MyDrive/YourFolder/your_file.gpkg\"\n",
    "uploaded_file = \"/content/drive/MyDrive/YourFolder/your_file.gpkg\"\n",
    "\n",
    "# Prompt the user to define the year of the desired wildfire season (from november/december of the selected year to march of the next year)\n",
    "season_year = 2021\n",
    "print(f\"The Images will be obtained for the {season_year}-{season_year+1} wildfire season\")\n",
    "\n",
    "# Prompt the user to provide the folder in Google Drive where outputs will be saved (e.g., GEE_FuelM)\n",
    "drive_folder = \"Folder Name\"\n",
    "\n",
    "# Define de CRS which the tiff files will be downloaded\n",
    "export_crs = \"EPSG:32718\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "fwvJqGKfqCIA"
   },
   "source": [
    "# Preparing Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Fvm9Biy6VMD3"
   },
   "source": [
    "## Loading the Region of Interest (ROI) and Defining Buffers\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "tzKnAlD9kuES"
   },
   "source": [
    "### Region of Interest"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "F645y3YZdJWZ"
   },
   "source": [
    "This section:\n",
    "- Loads the shapefile or GeoPackage into a GeoDataFrame (`gpd`).\n",
    "- Converts the GeoDataFrame into a Google Earth Engine `FeatureCollection`.\n",
    "  \n",
    "This FeatureCollection defines the area of interest (AOI) for processing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "LZ0Gzo9gVPVH"
   },
   "outputs": [],
   "source": [
    "# Load the vector file into a GeoDataFrame\n",
    "gdf = gpd.read_file(uploaded_file)\n",
    "print(\"Loaded Region of Interest (ROI):\")\n",
    "print(gdf)\n",
    "\n",
    "# Ensure CRS is EPSG:4326\n",
    "if gdf.crs != \"EPSG:4326\":\n",
    "    gdf = gdf.to_crs(epsg=4326)\n",
    "\n",
    "# Convert GeoDataFrame to Earth Engine FeatureCollection\n",
    "region = ee.FeatureCollection(gdf.__geo_interface__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "UC0IUVFDVria"
   },
   "source": [
    "### Buffers\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "DxRuqIs-dF0Q"
   },
   "source": [
    "Buffers are applied to the region to extend the area of analysis based on specific requirements:\n",
    "- **100m Buffer:** General buffer for region processing.\n",
    "- **5000m Buffer:** For climate data variables.\n",
    "- **500m Buffer:** For TPI (Topographic Position Index).\n",
    "- **200m Buffer:** For landforms.\n",
    "\n",
    "Each buffer is created using the `buffer()` method in GEE.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "2qFNhDEbVkrb"
   },
   "outputs": [],
   "source": [
    "# Create a 100m buffer around the region\n",
    "region_buffered = region.map(lambda f: f.buffer(100))\n",
    "\n",
    "# Define specific buffers for different datasets\n",
    "climate_buffer = region.map(lambda f: f.buffer(5000))\n",
    "tpi_buffer = region.map(lambda f: f.buffer(500))\n",
    "landform_buffer = region.map(lambda f: f.buffer(200))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "MQq9lHhyV2Pp"
   },
   "source": [
    "## Defining Analysis Periods and Helper Functions\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "2BWc62M-lAWg"
   },
   "source": [
    "### Analysis Periods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "UpHnYk-vWXa9"
   },
   "source": [
    "The script defines three periods for temporal analysis:\n",
    "- **Period 1 (P1):** 15 March (selected year) - 30 May (selected year).\n",
    "- **Period 2 (P2):** 01 June (selected year) - 15 August (selected year).\n",
    "- **Period 3 (P3):** 15 November (previous year) - 15 February (selected year).\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "kF5Y0sGI-jbd"
   },
   "outputs": [],
   "source": [
    "periods = [\n",
    "    (f'{season_year}-03-15', f'{season_year}-05-30'),    # Period 1 (P1)\n",
    "    (f'{season_year}-06-01', f'{season_year}-08-15'),    # Period 2 (P2)\n",
    "    (f'{season_year-1}-11-15', f'{season_year}-02-15')   # Period 3 (P3)\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "zdm290_jlVDR"
   },
   "source": [
    "### Helper Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "rjq9eun5coI3"
   },
   "source": [
    "- `apply_region_mask`: Clips an image to the region of interest (ROI) and assigns a NoData value to pixels outside the ROI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "mZbDmAKo76a1"
   },
   "outputs": [],
   "source": [
    "def apply_region_mask(image, region, nodata_value=-9999):\n",
    "    \"\"\"\n",
    "    Clips an image to the region of interest (ROI) and sets NoData values outside the ROI.\n",
    "\n",
    "    Args:\n",
    "        image: The Earth Engine image to clip and mask.\n",
    "        region: The region of interest to clip the image to.\n",
    "        nodata_value: The value to assign to NoData pixels.\n",
    "\n",
    "    Returns:\n",
    "        An Earth Engine image clipped to the ROI with NoData values assigned.\n",
    "    \"\"\"\n",
    "    clipped = image.clip(region)\n",
    "    masked = clipped.updateMask(clipped.mask()).unmask(nodata_value)\n",
    "    return masked"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "trD44a5TlPDk"
   },
   "source": [
    "A helper function (`get_period_name`) maps each period to its corresponding name (P1, P2, or P3)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "smgHtyTMlQyB"
   },
   "outputs": [],
   "source": [
    "def get_period_name(period):\n",
    "    \"\"\"\n",
    "    Returns the name of the period based on the start and end dates.\n",
    "\n",
    "    Args:\n",
    "        period: Tuple with start_date and end_date.\n",
    "\n",
    "    Returns:\n",
    "        Name of the period (e.g., 'P1', 'P2', 'P3').\n",
    "    \"\"\"\n",
    "    if period == periods[0]:\n",
    "        return \"P1\"\n",
    "    elif period == periods[1]:\n",
    "        return \"P2\"\n",
    "    elif period == periods[2]:\n",
    "        return \"P3\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Z-bmXjnbWA0T"
   },
   "source": [
    "## Data Processing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "CaGdEn4KV9Wr"
   },
   "source": [
    "### Processing Topographic Variables\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "U8vdb8s9rIXY"
   },
   "source": [
    "This section calculates and combines the following topographic variables:\n",
    "- Elevation (from NASADEM).\n",
    "- Slope (derived from elevation).\n",
    "- Aspect (orientation of slope).\n",
    "- Landforms (specific topographic features).\n",
    "- Topographic Position Index (TPI).\n",
    "\n",
    "All these variables are combined into a single image for export."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "4Qds8xGW78Dx"
   },
   "outputs": [],
   "source": [
    "# Elevation, slope, aspect\n",
    "elevation = apply_region_mask(ee.Image(\"NASA/NASADEM_HGT/001\").select('elevation').toFloat(), region_buffered)\n",
    "slope = apply_region_mask(ee.Terrain.slope(elevation).rename('slope').toFloat(), region_buffered)\n",
    "aspect = apply_region_mask(ee.Terrain.aspect(elevation).rename('aspect').toFloat(), region_buffered)\n",
    "\n",
    "# Landforms (200m buffer)\n",
    "landform = apply_region_mask(ee.Image('CSP/ERGo/1_0/Global/SRTM_landforms').select('constant').rename('landform').toFloat(), landform_buffer)\n",
    "\n",
    "# TPI (500m buffer)\n",
    "tpi = apply_region_mask(ee.Image('CSP/ERGo/1_0/Global/SRTM_mTPI').rename('TPI').toFloat(), tpi_buffer)\n",
    "\n",
    "# Combine all NASADEM variables into a single image\n",
    "nasadem_variables = apply_region_mask(ee.Image([elevation, aspect, slope]), region)\n",
    "\n",
    "# Clip Landform and Tpi Variables to region size\n",
    "landform = apply_region_mask(landform, region)\n",
    "tpi = apply_region_mask(tpi, region)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "LE5tIy-9WSDq"
   },
   "source": [
    "### Processing Climate Variables\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "vlUtD9HnrSq2"
   },
   "source": [
    "This section calculates:\n",
    "- **Total Precipitation:** Summed over a one-year period.\n",
    "- **Average Temperature:** Mean temperature over the same period.\n",
    "\n",
    "These variables are combined into a single image for export."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "6A3H8x497-CH"
   },
   "outputs": [],
   "source": [
    "# Climate data: total precipitation and average temperature\n",
    "climate = ee.ImageCollection(\"ECMWF/ERA5_LAND/MONTHLY_AGGR\") \\\n",
    "    .filter(ee.Filter.date(f'{season_year-1}-07-01', f'{season_year}-07-01'))\n",
    "\n",
    "precipitation = apply_region_mask(climate.select('total_precipitation_sum').sum().rename('precipitation').toFloat(), climate_buffer)\n",
    "temperature = apply_region_mask(climate.select('temperature_2m').mean().rename('temperature').toFloat(), climate_buffer)\n",
    "\n",
    "# Combine all climate variables into a single image\n",
    "climate_variables = apply_region_mask(ee.Image([precipitation, temperature]), region)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "EcDMs62JWgIb"
   },
   "source": [
    "### Sentinel-2 Processing\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "JU1xDtNYrbJh"
   },
   "source": [
    "- Fetches Sentinel-2 surface reflectance data with cloud masking applied.\n",
    "- Calculates spectral indices:\n",
    "  - NDVI (Normalized Difference Vegetation Index).\n",
    "  - NDBI (Normalized Difference Built-Up Index).\n",
    "  - NDWI (Normalized Difference Water Index).\n",
    "  - EVI (Enhanced Vegetation Index).\n",
    "\n",
    "  These indices are added as additional bands for further analysis.\n",
    "\n",
    "- Get the main variables of Sentinel-2:\n",
    "  - RGB (Red-Green-Blue)\n",
    "  - R1, R2, R3 and R4 (Red Edges)\n",
    "  - NIR (Near InfraRed)\n",
    "  - SWIR 1 and 2 (Short Wave InfraRed)\n",
    "\n",
    "The cloud correction function retrieves Sentinel-2 images and their cloud probability data, filters them by area, date, and cloud thresholds, and masks pixels with high cloud probability. It joins the reflectance and cloud data, applies cloud masks, and computes spectral indices like NDVI, NDBI, NDWI, and EVI. The cleaned images are then sorted by date and aggregated for analysis, ensuring only cloud-free data is used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "koMnpO8-BCHV"
   },
   "outputs": [],
   "source": [
    "# Define cloud threshold and mask probability\n",
    "scene_cloud_threshold = 60\n",
    "cloud_mask_probability = 30\n",
    "\n",
    "def add_indices(image):\n",
    "    \"\"\"\n",
    "    Adds spectral indices (NDVI, NDBI, NDWI, EVI) to a Sentinel-2 image.\n",
    "    \"\"\"\n",
    "    ndbi = image.expression('(SWIR - NIR) / (SWIR + NIR)', {\n",
    "        'SWIR': image.select('swir1'),\n",
    "        'NIR': image.select('nir'),\n",
    "    }).multiply(100).rename('ndbi')\n",
    "\n",
    "    ndvi = image.normalizedDifference(['nir', 'red']).multiply(100).rename('ndvi')\n",
    "    ndwi = image.normalizedDifference(['green', 'nir']).multiply(100).rename('ndwi')\n",
    "    evi = image.expression(\n",
    "        '2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))', {\n",
    "            'NIR': image.select('nir'),\n",
    "            'RED': image.select('red'),\n",
    "            'BLUE': image.select('blue'),\n",
    "        }).rename('evi')\n",
    "\n",
    "    return image.addBands([ndvi, ndbi, ndwi, evi])\n",
    "\n",
    "def mask_clouds(cloud_probability_threshold):\n",
    "    \"\"\"\n",
    "    Masks clouds based on a given probability threshold.\n",
    "    \"\"\"\n",
    "    def _mask_image(img):\n",
    "        cloud_mask = img.select('probability').lt(cloud_probability_threshold)\n",
    "        return img.updateMask(cloud_mask)\n",
    "    return _mask_image\n",
    "\n",
    "def get_s2_sr_cloud_probability(aoi, start_date, end_date, scene_cloud_threshold, cloud_mask_probability):\n",
    "    \"\"\"\n",
    "    Fetches Sentinel-2 surface reflectance images with cloud masking.\n",
    "    \"\"\"\n",
    "    primary = ee.ImageCollection(\"COPERNICUS/S2_SR_HARMONIZED\") \\\n",
    "        .filterBounds(aoi) \\\n",
    "        .filterDate(start_date, end_date) \\\n",
    "        .filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', scene_cloud_threshold) \\\n",
    "        .select(['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B8A', 'B11', 'B12'],\n",
    "                ['blue', 'green', 'red', 'R1', 'R2', 'R3', 'nir', 'R4', 'swir1', 'swir2']) \\\n",
    "        .map(add_indices)\n",
    "\n",
    "    secondary = ee.ImageCollection(\"COPERNICUS/S2_CLOUD_PROBABILITY\") \\\n",
    "        .filterBounds(aoi) \\\n",
    "        .filterDate(start_date, end_date)\n",
    "\n",
    "    joined = ee.Join.inner().apply(primary=primary, secondary=secondary,\n",
    "                                   condition=ee.Filter.equals(leftField='system:index', rightField='system:index'))\n",
    "\n",
    "    def merge_image_bands(join_result):\n",
    "        return ee.Image(join_result.get('primary')).addBands(join_result.get('secondary'))\n",
    "\n",
    "    return ee.ImageCollection(joined.map(merge_image_bands)) \\\n",
    "        .map(mask_clouds(cloud_mask_probability)) \\\n",
    "        .sort('system:time_start')\n",
    "\n",
    "def process_sentinel2_index_variables(period, region):\n",
    "    start_date, end_date = period\n",
    "    sentinel2 = get_s2_sr_cloud_probability(\n",
    "        region, start_date, end_date, scene_cloud_threshold, cloud_mask_probability\n",
    "    )\n",
    "    vars = {\n",
    "        'ndvi': apply_region_mask(sentinel2.select('ndvi').median().rename(f'ndvi_{get_period_name(period)}').toFloat(), region),\n",
    "        'ndbi': apply_region_mask(sentinel2.select('ndbi').median().rename(f'ndbi_{get_period_name(period)}').toFloat(), region),\n",
    "        'ndwi': apply_region_mask(sentinel2.select('ndwi').median().rename(f'ndwi_{get_period_name(period)}').toFloat(), region),\n",
    "        'evi': apply_region_mask(sentinel2.select('evi').median().rename(f'evi_{get_period_name(period)}').toFloat(), region)\n",
    "    }\n",
    "\n",
    "    return ee.Image([vars['ndvi'], vars['ndbi'], vars['ndwi'], vars['evi']])\n",
    "\n",
    "def process_sentinel2(period, region):\n",
    "    start_date, end_date = period\n",
    "    sentinel2 = get_s2_sr_cloud_probability(\n",
    "        region, start_date, end_date, scene_cloud_threshold, cloud_mask_probability\n",
    "    )\n",
    "\n",
    "    vars = {\n",
    "        'red': apply_region_mask(sentinel2.select('red').median().rename(f'red_{get_period_name(period)}').toFloat(), region),\n",
    "        'green': apply_region_mask(sentinel2.select('green').median().rename(f'green_{get_period_name(period)}').toFloat(), region),\n",
    "        'blue': apply_region_mask(sentinel2.select('blue').median().rename(f'blue_{get_period_name(period)}').toFloat(), region),\n",
    "        'R1': apply_region_mask(sentinel2.select('R1').median().rename(f'R1_{get_period_name(period)}').toFloat(), region),\n",
    "        'R2': apply_region_mask(sentinel2.select('R2').median().rename(f'R2_{get_period_name(period)}').toFloat(), region),\n",
    "        'R3': apply_region_mask(sentinel2.select('R3').median().rename(f'R3_{get_period_name(period)}').toFloat(), region),\n",
    "        'nir': apply_region_mask(sentinel2.select('nir').median().rename(f'nir_{get_period_name(period)}').toFloat(), region),\n",
    "        'R4': apply_region_mask(sentinel2.select('R4').median().rename(f'R4_{get_period_name(period)}').toFloat(), region),\n",
    "        'swir1': apply_region_mask(sentinel2.select('swir1').median().rename(f'swir1_{get_period_name(period)}').toFloat(), region),\n",
    "        'swir2': apply_region_mask(sentinel2.select('swir2').median().rename(f'swir2_{get_period_name(period)}').toFloat(), region)\n",
    "    }\n",
    "\n",
    "    return ee.Image([vars['red'], vars['green'], vars['blue'], vars['R1'], vars['R2'], vars['R3'], vars['nir'], vars['R4'], vars['swir1'], vars['swir2']])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "PlKmP-1EWo4X"
   },
   "source": [
    "### Sentinel-1 Processing\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "vYDgNz7orgm1"
   },
   "source": [
    "This section processes Sentinel-1 radar data:\n",
    "- Applies radiometric slope correction to account for terrain effects.\n",
    "- Orbit propierties of the satelite are set to **\"Descending\"**\n",
    "- Generates two bands:\n",
    "  - **VV Polarization:** Vertical transmit and receive.\n",
    "  - **VH Polarization:** Vertical transmit, horizontal receive.\n",
    "\n",
    "The `process_sentinel1` function computes the mean value for each band, clips it to the region of interest, and names it based on the period.\n",
    "\n",
    "The `slope_correction` function adjusts Sentinel-1 radar images for terrain-induced distortions. It uses a digital elevation model (DEM) to calculate the terrain's slope and aspect, which are then used to correct for the angle of radar backscatter. The function applies a volume model correction to normalize radar signals and creates masks for layover (areas where slopes face the sensor, causing signal overlap) and shadow (areas blocked from radar signals). These corrections ensure the radar data reflects true ground conditions, reducing distortions caused by terrain. The processed bands (VV and VH polarizations) are averaged, clipped to the region of interest, and prepared for further analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "I7DD8TVRBE8_"
   },
   "outputs": [],
   "source": [
    "def slope_correction(collection, elevation=None, model='volume', buffer=10):\n",
    "    \"\"\"\n",
    "    Applies radiometric slope correction to a Sentinel-1 collection.\n",
    "\n",
    "    Args:\n",
    "        collection: ee.ImageCollection of Sentinel-1 images.\n",
    "        elevation: ee.Image of DEM (optional, default is NASADEM).\n",
    "        model: Correction model to apply ('volume' or 'surface').\n",
    "        buffer: Buffer distance in meters for layover/shadow mask (default is 10).\n",
    "\n",
    "    Returns:\n",
    "        ee.ImageCollection with corrected images and additional bands.\n",
    "    \"\"\"\n",
    "    if elevation is None:\n",
    "        elevation = ee.Image('NASA/NASADEM_HGT/001')\n",
    "\n",
    "    def _volumetric_model_SCF(theta_iRad, alpha_rRad):\n",
    "        \"\"\"\n",
    "        Calculates the volumetric model SCF.\n",
    "        \"\"\"\n",
    "        ninetyRad = ee.Image.constant(90).multiply(3.14159265359 / 180)\n",
    "        nominator = (ninetyRad.subtract(theta_iRad).add(alpha_rRad)).tan()\n",
    "        denominator = (ninetyRad.subtract(theta_iRad)).tan()\n",
    "        return nominator.divide(denominator)\n",
    "\n",
    "    def _masking(alpha_rRad, theta_iRad, buffer):\n",
    "        \"\"\"\n",
    "        Creates masks for layover and shadow.\n",
    "        \"\"\"\n",
    "        layover = alpha_rRad.lt(theta_iRad).rename('layover')\n",
    "        ninetyRad = ee.Image.constant(90).multiply(3.14159265359 / 180)\n",
    "        shadow = alpha_rRad.gt(ee.Image.constant(-1).multiply(ninetyRad.subtract(theta_iRad))).rename('shadow')\n",
    "\n",
    "        if buffer > 0:\n",
    "            layover = layover.Not().fastDistanceTransform(30).sqrt().multiply(ee.Image.pixelArea().sqrt()).gt(buffer).rename('layover')\n",
    "            shadow = shadow.Not().fastDistanceTransform(30).sqrt().multiply(ee.Image.pixelArea().sqrt()).gt(buffer).rename('shadow')\n",
    "\n",
    "        no_data_mask = layover.And(shadow).rename('no_data_mask')\n",
    "        return layover.addBands(shadow).addBands(no_data_mask)\n",
    "\n",
    "    def _correct(image):\n",
    "        \"\"\"\n",
    "        Applies slope correction to a single image and adds layover and shadow masks.\n",
    "        \"\"\"\n",
    "        theta_iRad = image.select('angle').multiply(3.14159265359 / 180)\n",
    "        alpha_sRad = ee.Terrain.slope(elevation).multiply(3.14159265359 / 180)\n",
    "        phi_sRad = ee.Terrain.aspect(elevation).multiply(3.14159265359 / 180)\n",
    "\n",
    "        phi_iRad = ee.Image.constant(0).multiply(3.14159265359 / 180)  # Assuming flat incidence direction\n",
    "\n",
    "        phi_rRad = phi_iRad.subtract(phi_sRad)\n",
    "        alpha_rRad = (alpha_sRad.tan().multiply(phi_rRad.cos())).atan()\n",
    "\n",
    "        gamma0 = image.divide(10.0).pow(10).divide(theta_iRad.cos())\n",
    "\n",
    "        if model == 'volume':\n",
    "            scf = _volumetric_model_SCF(theta_iRad, alpha_rRad)\n",
    "        else:\n",
    "            raise ValueError(\"Only 'volume' model is supported in this implementation.\")\n",
    "\n",
    "        gamma0_flat = gamma0.divide(scf)\n",
    "        gamma0_flatDB = ee.Image.constant(10).multiply(gamma0_flat.log10()).rename(image.bandNames())\n",
    "\n",
    "        masks = _masking(alpha_rRad, theta_iRad, buffer)\n",
    "\n",
    "        return gamma0_flatDB.addBands(masks).copyProperties(image, image.propertyNames())\n",
    "\n",
    "    return collection.map(_correct)\n",
    "\n",
    "\n",
    "def process_sentinel1_band(period, region, band_name):\n",
    "    \"\"\"\n",
    "    Processes a specific Sentinel-1 band, applying slope correction and clipping it to the region.\n",
    "\n",
    "    Args:\n",
    "        period: Tuple with the start and end dates of the period (start_date, end_date).\n",
    "        region: Region of interest as an `ee.FeatureCollection`.\n",
    "        band_name: Name of the band ('VV' or 'VH').\n",
    "\n",
    "    Returns:\n",
    "        Processed and clipped image of the selected band, or None if no valid data is available.\n",
    "    \"\"\"\n",
    "    start_date, end_date = period\n",
    "\n",
    "    # Filter the Sentinel-1 collection\n",
    "    collection = ee.ImageCollection(\"COPERNICUS/S1_GRD\") \\\n",
    "        .filterDate(start_date, end_date) \\\n",
    "        .filterBounds(region) \\\n",
    "        .filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING')) \\\n",
    "        .filter(ee.Filter.listContains('transmitterReceiverPolarisation', band_name)) \\\n",
    "        .filter(ee.Filter.eq('instrumentMode', 'IW'))\n",
    "\n",
    "    # Check if the collection is empty\n",
    "    collection_size = collection.size().getInfo()\n",
    "    if collection_size == 0:\n",
    "        print(f\"No images found for band {band_name} in period {start_date} to {end_date}\")\n",
    "        return None\n",
    "\n",
    "    # Apply slope correction\n",
    "    corrected = slope_correction(collection)\n",
    "\n",
    "    # Select the corrected band, calculate the mean, and clip it to the region\n",
    "    band_mean = corrected.select(band_name).mean()\n",
    "    band_clipped = band_mean.clip(region).rename(f'{band_name}_{get_period_name(period)}').toFloat()\n",
    "\n",
    "    return band_clipped\n",
    "\n",
    "\n",
    "def process_sentinel1(period, region):\n",
    "    \"\"\"\n",
    "    Processes Sentinel-1 for VV and VH bands, applying slope correction and clipping them to the region.\n",
    "\n",
    "    Args:\n",
    "        period: Tuple with the start and end dates of the period (start_date, end_date).\n",
    "        region: Region of interest as an `ee.FeatureCollection`.\n",
    "\n",
    "    Returns:\n",
    "        Image with processed VV and VH bands, or None if no valid data is available.\n",
    "    \"\"\"\n",
    "    vh = process_sentinel1_band(period, region, 'VH')\n",
    "    vv = process_sentinel1_band(period, region, 'VV')\n",
    "\n",
    "    # Check if any band is None\n",
    "    if vh is None or vv is None:\n",
    "        print(f\"Skipping Sentinel-1 processing for period {get_period_name(period)} due to missing data.\")\n",
    "        return None\n",
    "\n",
    "    return ee.Image([vh, vv])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "aaJ2IzR8WvIM"
   },
   "source": [
    "### EVI Statistics\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "VyoJ-Mrerkpz"
   },
   "source": [
    "Calculates summary statistics (sum, min, max, standard deviation) for the Enhanced Vegetation Index (EVI) over the region of interest for each period."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "DpeDTsnX-fQc"
   },
   "outputs": [],
   "source": [
    "def calculate_evi_statistics_by_period(period, region):\n",
    "    start_date, end_date = period\n",
    "    sentinel2 = get_s2_sr_cloud_probability(\n",
    "        region, start_date, end_date, scene_cloud_threshold, cloud_mask_probability\n",
    "    )\n",
    "\n",
    "    sentinel2_evi = sentinel2.map(\n",
    "    lambda img: img.expression(\n",
    "        '2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))',\n",
    "        {\n",
    "            'NIR': img.select('nir'),\n",
    "            'RED': img.select('red'),\n",
    "            'BLUE': img.select('blue'),\n",
    "        }\n",
    "    ).rename('evi')\n",
    "    )\n",
    "\n",
    "    vars = {\n",
    "        'evi_min': apply_region_mask(sentinel2_evi.reduce(ee.Reducer.min()).rename(f'evi_min_{get_period_name(period)}').toFloat(), region),\n",
    "        'evi_max': apply_region_mask(sentinel2_evi.reduce(ee.Reducer.max()).rename(f'evi_max_{get_period_name(period)}').toFloat(), region),\n",
    "        'evi_sd': apply_region_mask(sentinel2_evi.reduce(ee.Reducer.stdDev()).rename(f'evi_sd_{get_period_name(period)}').toFloat(), region)\n",
    "    }\n",
    "\n",
    "    return ee.Image([vars['evi_min'], vars['evi_max'], vars['evi_sd']])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "QkUo6gh3Wx2b"
   },
   "source": [
    "## Adding Data to Collections\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "E-8tsbtWc8Qy"
   },
   "source": [
    "For each period, the script:\n",
    "- Processes Sentinel-2 data and adds it to `sentinel2_collections`.\n",
    "- Processes Sentinel-2 indices and adds them to `sentinel2_index_variables_collections`.\n",
    "- Processes Sentinel-1 data and adds it to `sentinel1_collections`.\n",
    "- Calculates EVI statistics and adds them to `evi_statistics`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "4Woo2k8hBI9Q"
   },
   "outputs": [],
   "source": [
    "sentinel2_collections = []\n",
    "sentinel2_index_variables_collections = []\n",
    "sentinel1_collections = []\n",
    "evi_statistics = []\n",
    "\n",
    "for period in periods:\n",
    "    sentinel2_collections.append(apply_region_mask(process_sentinel2(period, region_buffered), region))\n",
    "    sentinel2_index_variables_collections.append(apply_region_mask(process_sentinel2_index_variables(period, region_buffered), region))\n",
    "    sentinel1_collections.append(apply_region_mask(process_sentinel1(period, region_buffered), region))\n",
    "    evi_statistics.append(apply_region_mask(calculate_evi_statistics_by_period(period, region_buffered), region))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "fOm-NtKq_tjm"
   },
   "source": [
    "## Save Data into Pickle file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "p0gZpst5ZRZR"
   },
   "outputs": [],
   "source": [
    "# Serialize the variables and lists\n",
    "data = {\n",
    "    'landform': landform,\n",
    "    'tpi': tpi,\n",
    "    'nasadem_variables': nasadem_variables,\n",
    "    'climate_variables': climate_variables,\n",
    "    'sentinel2_collections': sentinel2_collections,\n",
    "    'sentinel2_index_variables_collections': sentinel2_index_variables_collections,\n",
    "    'sentinel1_collections': sentinel1_collections,\n",
    "    'evi_statistics': evi_statistics,\n",
    "}\n",
    "\n",
    "# Save the Data into a file\n",
    "with open('/content/processed_data.pkl', 'wb') as f:\n",
    "    pickle.dump(data, f)\n",
    "\n",
    "print(\"Data Serilized and saved as processed_data.pkl\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "lrxcvUHlZYeM"
   },
   "outputs": [],
   "source": [
    "# Copiar el archivo serializado a Google Drive\n",
    "!cp /content/processed_data.pkl /content/drive/MyDrive/\n",
    "print(\"File processed_data.pkl uploaded to Google Drive.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "2eoMNfC5W413"
   },
   "source": [
    "# Exporting Data\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "EYWoqIgG7MGs"
   },
   "source": [
    "## Load Region"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "HeGEWJ-qjDG5"
   },
   "outputs": [],
   "source": [
    "# Load the vector file into a GeoDataFrame\n",
    "gdf = gpd.read_file(uploaded_file)\n",
    "print(\"Loaded Region of Interest (ROI):\")\n",
    "print(gdf)\n",
    "\n",
    "# Convert GeoDataFrame to Earth Engine FeatureCollection\n",
    "region = ee.FeatureCollection(gdf.__geo_interface__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "QXGCF6bH7VxW"
   },
   "source": [
    "## Load the variables and collections from the pickle file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ToCv0HxsgWMv"
   },
   "outputs": [],
   "source": [
    "# Load the serialized data from Google Drive\n",
    "with open('/content/drive/MyDrive/processed_data.pkl', 'rb') as f:\n",
    "    data = pickle.load(f)\n",
    "\n",
    "# Access the variables\n",
    "landform = data['landform']\n",
    "tpi = data['tpi']\n",
    "nasadem_variables = data['elevation_variables']\n",
    "climate_variables = data['climate_variables']\n",
    "sentinel2_collections = data['sentinel2_collections']\n",
    "sentinel2_index_variables_collections = data['sentinel2_index_variables_collections']\n",
    "sentinel1_collections = data['sentinel1_collections']\n",
    "evi_statistics = data['evi_statistics']\n",
    "\n",
    "print(\"Data loaded from the pickle file\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "8CvryX1vdAcC"
   },
   "source": [
    "Exports the processed data to Google Drive in batches:\n",
    "1. **Batch 1:** NASADEM Variables.\n",
    "2. **Batch 2:** TPI, Landform an Climate Variables.\n",
    "3. **Batch 3:** Sentinel-1, Sentinel-2 with EVI statistics and Index Variables for Period 1.\n",
    "4. **Batch 4:** Sentinel-1, Sentinel-2 with EVI statistics and Index Variables for Period 2.\n",
    "5. **Batch 5:** Sentinel-1, Sentinel-2 with EVI statistics and Index Variables for Period 3.\n",
    "Each batch waits for the previous one to complete before starting."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "8rV1DoYO7gWO"
   },
   "source": [
    "## Export Data to Drive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "-O6X7UiTVIZq"
   },
   "outputs": [],
   "source": [
    "# Define a function to export images and wait for the tasks to complete\n",
    "def export_batch(batch_tasks, defined_crs):\n",
    "    \"\"\"\n",
    "    Export a batch of tasks and wait for them to complete before proceeding.\n",
    "\n",
    "    Args:\n",
    "        batch_tasks: List of tuples, each containing (image, description, file_name, region).\n",
    "    \"\"\"\n",
    "    active_tasks = []\n",
    "\n",
    "    for image, description, file_name, region in batch_tasks:\n",
    "        print(f\"Starting export of {description} to Google Drive...\")\n",
    "        task = ee.batch.Export.image.toDrive(\n",
    "            image=image,\n",
    "            description=description,\n",
    "            folder=drive_folder,\n",
    "            fileNamePrefix=file_name,\n",
    "            region=region.geometry(),\n",
    "            scale=30,\n",
    "            crs=defined_crs,\n",
    "            maxPixels=1e13,\n",
    "            fileFormat='GeoTIFF',\n",
    "            shardSize=1024,\n",
    "            formatOptions={\n",
    "              'noData': -9999\n",
    "            }\n",
    "        )\n",
    "        task.start()\n",
    "        active_tasks.append((task, description))\n",
    "\n",
    "    # Monitor tasks until all are complete\n",
    "    for task, description in active_tasks:\n",
    "        print(f\"Waiting for task {description} to complete...\")\n",
    "        while task.active():\n",
    "            time.sleep(60)  # Check every minute\n",
    "        status = task.status()\n",
    "        if status['state'] == 'COMPLETED':\n",
    "            print(f\"Task {description} completed successfully.\")\n",
    "        else:\n",
    "            print(f\"Task {description} failed: {status.get('error_message', 'No error message provided')}.\")\n",
    "\n",
    "# Define the 5 batches for exports\n",
    "batch_1 = [\n",
    "    (nasadem_variables, 'NASADEM Variables', 'nasadem_variables', region)                                                 # 17 minutes aprox (+ 15 minutes to change crs)\n",
    "]\n",
    "\n",
    "batch_2 = [\n",
    "    (tpi, 'TPI', 'tpi', region),                                                                                          # 05 minutes aprox\n",
    "    (landform, 'Landform', 'landform', region),                                                                           # 03 minutes aprox\n",
    "    (climate_variables, 'Climate Variables', 'climate_variables', region),                                                # 06 minutes aprox\n",
    "]\n",
    "\n",
    "batch_3 = [\n",
    "    (sentinel2_collections[0], 'Sentinel-2 P1', 'sentinel-2_P1', region),                                                 # 27 minutes aprox\n",
    "    (sentinel2_index_variables_collections[0], 'Sentinel-2 Index Variables P1', 'sentinel-2-index-variables_P1', region), # 12 minutes aprox\n",
    "    (evi_statistics[0], 'EVI Stats P1', 'evi-stats_P1', region),                                                          # 11 minutes aprox\n",
    "    (sentinel1_collections[0], 'Sentinel-1 P1', 'sentinel-1_P1', region)                                                  # 11 minutes aprox\n",
    "]\n",
    "\n",
    "batch_4 = [\n",
    "    (sentinel2_collections[1], 'Sentinel-2 P2', 'sentinel-2_P2', region),                                                 # 20 minutes aprox\n",
    "    (sentinel2_index_variables_collections[1], 'Sentinel-2 Index Variables P2', 'sentinel-2-index-variables_P2', region), # 11 minutes aprox\n",
    "    (evi_statistics[1], 'EVI Stats P2', 'evi-stats_P2', region),                                                          # 07 minutes aprox\n",
    "    (sentinel1_collections[1], 'Sentinel-1 P2', 'sentinel-1_P2', region)                                                  # 13 minutes aprox\n",
    "]\n",
    "\n",
    "batch_5 = [\n",
    "    (sentinel2_collections[2], 'Sentinel-2 P3', 'sentinel-2_P3', region),                                                 # 17 minutes aprox\n",
    "    (sentinel2_index_variables_collections[2], 'Sentinel-2 Index Variables P3', 'sentinel-2-index-variables_P3', region), # 16 minutes aprox\n",
    "    (evi_statistics[2], 'EVI Stats P3', 'evi-stats_P3', region),                                                          # 11 minutes aprox\n",
    "    (sentinel1_collections[2], 'Sentinel-1 P3', 'sentinel-1_P3', region)                                                  # 11 minutes aprox\n",
    "]                                                                                                                  #Total: 213 minutes aprox (5,8 gb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "9YCCzN85Zhy3"
   },
   "source": [
    "### Download Batches"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "tfCu6YuAXzXf"
   },
   "outputs": [],
   "source": [
    "print(\"Starting Batch 1\")\n",
    "export_batch(batch_1, 'EPSG:4326') # Do not change this CRS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "eOq0BWEnXzvc"
   },
   "outputs": [],
   "source": [
    "print(\"Starting Batch 2\")\n",
    "export_batch(batch_2, export_crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "8mKCu3wXXz0l"
   },
   "outputs": [],
   "source": [
    "print(\"Starting Batch 3\")\n",
    "export_batch(batch_3, export_crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "QF4p6w7fXryk"
   },
   "outputs": [],
   "source": [
    "print(\"Starting Batch 4\")\n",
    "export_batch(batch_4, export_crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "068EOT8ZVQdg"
   },
   "outputs": [],
   "source": [
    "print(\"Starting Batch 5\")\n",
    "export_batch(batch_5, export_crs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "mqS2cdUvaZtC"
   },
   "source": [
    "#### Convert NASADEM variables to desired CRS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Hymz-QeYlz-g"
   },
   "outputs": [],
   "source": [
    "# Path to the original NASADEM file\n",
    "nasadem_tif = f'/content/drive/MyDrive/{drive_folder}/nasadem_variables.tif'  # Change this path as needed\n",
    "\n",
    "# Path to a reference raster (any correctly aligned file)\n",
    "reference_tif = f'/content/drive/MyDrive/{drive_folder}/evi-stats_P2.tif'  # Change this path to a correctly aligned file\n",
    "\n",
    "# Temporary file for reprojection\n",
    "temp_nasadem = nasadem_tif.replace('.tif', '_temp.tif')\n",
    "\n",
    "try:\n",
    "    # Open the reference raster to get the correct extent, width, height, and transform\n",
    "    with rasterio.open(reference_tif) as ref_src:\n",
    "        ref_transform = ref_src.transform\n",
    "        ref_width = ref_src.width\n",
    "        ref_height = ref_src.height\n",
    "        ref_crs = ref_src.crs\n",
    "        ref_bounds = ref_src.bounds\n",
    "\n",
    "    # Open the NASADEM raster\n",
    "    with rasterio.open(nasadem_tif) as src:\n",
    "        src_crs = src.crs\n",
    "\n",
    "        # Retrieve original band names\n",
    "        band_names = src.descriptions if src.descriptions else [f'Band {i}' for i in range(1, src.count + 1)]\n",
    "\n",
    "        # Update metadata with the correct width, height, transform, and CRS\n",
    "        kwargs = src.meta.copy()\n",
    "        kwargs.update({\n",
    "            'crs': ref_crs,\n",
    "            'transform': ref_transform,\n",
    "            'width': ref_width,\n",
    "            'height': ref_height,\n",
    "            'nodata': -9999,\n",
    "            'compress': 'DEFLATE',\n",
    "            'predictor': 2,\n",
    "            'zlevel': 9\n",
    "        })\n",
    "\n",
    "        # Create the reprojected NASADEM file\n",
    "        with rasterio.open(temp_nasadem, 'w', **kwargs) as dst:\n",
    "            for i in range(1, src.count + 1):\n",
    "                # Create an empty array for the reprojected data\n",
    "                dst_array = np.empty((ref_height, ref_width), dtype=src.dtypes[i - 1])\n",
    "\n",
    "                # Perform the reprojection\n",
    "                reproject(\n",
    "                    source=rasterio.band(src, i),\n",
    "                    destination=dst_array,\n",
    "                    src_transform=src.transform,\n",
    "                    src_crs=src_crs,\n",
    "                    dst_transform=ref_transform,\n",
    "                    dst_crs=ref_crs,\n",
    "                    resampling=Resampling.bilinear\n",
    "                )\n",
    "\n",
    "                # Write the reprojected data to the output file\n",
    "                dst.write(dst_array, i)\n",
    "\n",
    "                # Restore band names\n",
    "                dst.set_band_description(i, band_names[i - 1])\n",
    "\n",
    "    # Replace the original NASADEM file with the corrected one\n",
    "    os.remove(nasadem_tif)  # Remove the original file\n",
    "    os.rename(temp_nasadem, nasadem_tif)  # Rename the temporary file\n",
    "\n",
    "    print(f\"NASADEM file reprojected and aligned successfully: {nasadem_tif}\")\n",
    "\n",
    "except Exception as e:\n",
    "    print(f\"Error reprojecting NASADEM file: {e}\")"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [
    "NfAVArbkh3_r",
    "Zb7mMmCeUcdq",
    "hB4oty-qU-wf",
    "pLjQci1OeSVj",
    "yM6M0xq3ekSK",
    "C9UlPJRtf_zR",
    "fwvJqGKfqCIA",
    "Fvm9Biy6VMD3",
    "tzKnAlD9kuES",
    "UC0IUVFDVria",
    "MQq9lHhyV2Pp",
    "2BWc62M-lAWg",
    "zdm290_jlVDR",
    "Z-bmXjnbWA0T",
    "CaGdEn4KV9Wr",
    "LE5tIy-9WSDq",
    "EcDMs62JWgIb",
    "PlKmP-1EWo4X",
    "aaJ2IzR8WvIM",
    "QkUo6gh3Wx2b",
    "fOm-NtKq_tjm",
    "2eoMNfC5W413",
    "EYWoqIgG7MGs",
    "QXGCF6bH7VxW",
    "8rV1DoYO7gWO",
    "9YCCzN85Zhy3"
   ],
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}