mask tester.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "bc8ec45c-83cc-443d-bacc-a94aa372dfc3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import glob\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import json\n",
    "import rasterio\n",
    "from rasterio.features import rasterize\n",
    "import geopandas as gpd\n",
    "from shapely.geometry import box\n",
    "import matplotlib.pyplot as plt\n",
    "from rasterio.crs import CRS\n",
    "from rasterio.warp import calculate_default_transform, reproject, Resampling\n",
    "from spectral import open_image\n",
    "import warnings\n",
    "\n",
    "# Suppress NotGeoreferencedWarning for clarity\n",
    "warnings.filterwarnings(\"ignore\", category=rasterio.errors.NotGeoreferencedWarning)\n",
    "\n",
    "# ----- Existing ENVI Processing Classes and Functions -----\n",
    "\n",
    "class ENVIProcessor:\n",
    "    def __init__(self, file_path):\n",
    "        self.file_path = file_path\n",
    "        self.data = None  # This will hold the raster data array\n",
    "        self.file_type = \"envi\"\n",
    "\n",
    "    def load_data(self):\n",
    "        \"\"\"Loads the raster data from the file_path into self.data\"\"\"\n",
    "        with rasterio.open(self.file_path) as src:\n",
    "            self.data = src.read()  # Read all bands\n",
    "\n",
    "    def get_chunk_from_extent(self, corrections=[], resample=False):\n",
    "        self.load_data()  # Ensure data is loaded\n",
    "        return self.data\n",
    "\n",
    "def find_raster_files(directory):\n",
    "    \"\"\"\n",
    "    Searches for raster files in the given directory, capturing both original and corrected ENVI files,\n",
    "    plus resampled ones, while excluding .hdr, .json, .csv, and any files containing '_mask' or '_ancillary'.\n",
    "    We'll look for filenames containing '_reflectance' (original) or '_envi' (corrected/resampled).\n",
    "    \"\"\"\n",
    "    pattern = \"*\"\n",
    "    full_pattern = os.path.join(directory, pattern)\n",
    "    all_files = glob.glob(full_pattern)\n",
    "\n",
    "    filtered_files = [\n",
    "        file for file in all_files\n",
    "        if (\n",
    "            ('_reflectance' in os.path.basename(file) or '_envi' in os.path.basename(file)) and\n",
    "            '_mask' not in os.path.basename(file) and\n",
    "            '_ancillary' not in os.path.basename(file) and\n",
    "            not file.endswith('.hdr') and\n",
    "            not file.endswith('.json') and\n",
    "            not file.endswith('.csv')\n",
    "        )\n",
    "    ]\n",
    "\n",
    "    found_files_set = set(filtered_files)\n",
    "    found_files = list(found_files_set)\n",
    "    found_files.sort()\n",
    "\n",
    "    return found_files\n",
    "\n",
    "def load_and_combine_rasters(raster_paths):\n",
    "    \"\"\"\n",
    "    Loads and combines raster data from a list of file paths.\n",
    "    Assumes each raster has shape (bands, rows, cols) and that\n",
    "    all rasters can be concatenated along the band dimension.\n",
    "    \"\"\"\n",
    "    chunks = []\n",
    "    for path in raster_paths:\n",
    "        processor = ENVIProcessor(path)\n",
    "        chunk = processor.get_chunk_from_extent(corrections=['some_correction'], resample=False)\n",
    "        chunks.append(chunk)\n",
    "    combined_array = np.concatenate(chunks, axis=0)  # Combine along the first axis (bands)\n",
    "    return combined_array\n",
    "\n",
    "def process_and_flatten_array(array, json_dir='Resampling', original_bands=426, corrected_bands=426,\n",
    "                              original_wavelengths=None, corrected_wavelengths=None, folder_name=None,\n",
    "                              map_info=None):\n",
    "    \"\"\"\n",
    "    Processes a 3D numpy array to a DataFrame, adds metadata columns, \n",
    "    renames columns dynamically based on JSON configuration, and adds Pixel_id.\n",
    "    Uses provided wavelength lists to name original and corrected bands, and includes geocoordinates.\n",
    "\n",
    "    Parameters:\n",
    "    - array: A 3D numpy array of shape (bands, rows, cols).\n",
    "    - json_dir: Directory containing the landsat_band_parameters.json file.\n",
    "    - original_bands: Number of original bands expected.\n",
    "    - corrected_bands: Number of corrected bands expected.\n",
    "    - original_wavelengths: List of wavelengths for the original bands (floats).\n",
    "    - corrected_wavelengths: List of wavelengths for the corrected bands (floats).\n",
    "    - folder_name: Name of the subdirectory (flight line identifier).\n",
    "    - map_info: The map info array from the metadata for georeferencing.\n",
    "\n",
    "    Returns:\n",
    "    - A pandas DataFrame with additional metadata columns and renamed band columns.\n",
    "    \"\"\"\n",
    "    if len(array.shape) != 3:\n",
    "        raise ValueError(\"Input array must be 3-dimensional. Expected (bands, rows, cols).\")\n",
    "\n",
    "    json_file = os.path.join(json_dir, 'landsat_band_parameters.json')\n",
    "    if not os.path.isfile(json_file):\n",
    "        raise FileNotFoundError(f\"JSON file not found: {json_file}\")\n",
    "\n",
    "    with open(json_file, 'r') as f:\n",
    "        config = json.load(f)\n",
    "\n",
    "    bands, rows, cols = array.shape\n",
    "    print(f\"[DEBUG] array shape: bands={bands}, rows={rows}, cols={cols}\")\n",
    "\n",
    "    reshaped_array = array.reshape(bands, -1).T  # (pixels, bands)\n",
    "    pixel_indices = np.indices((rows, cols)).reshape(2, -1).T  # (pixels, 2)\n",
    "    df = pd.DataFrame(reshaped_array, columns=[f'Band_{i+1}' for i in range(bands)])\n",
    "\n",
    "    # Extract map info for georeferencing:\n",
    "    # Format: [projection, x_pixel_start, y_pixel_start, map_x, map_y, x_res, y_res, ...]\n",
    "    # Typically:\n",
    "    #   x_pixel_start, y_pixel_start = 1,1 for upper-left pixel\n",
    "    #   map_x, map_y = coordinates of that upper-left pixel\n",
    "    #   x_res, y_res = pixel sizes (y_res should be positive but we assume north-down in ENVI)\n",
    "    if map_info is not None and len(map_info) >= 7:\n",
    "        projection = map_info[0]\n",
    "        x_pixel_start = float(map_info[1])\n",
    "        y_pixel_start = float(map_info[2])\n",
    "        map_x = float(map_info[3])\n",
    "        map_y = float(map_info[4])\n",
    "        x_res = float(map_info[5])\n",
    "        y_res = float(map_info[6])\n",
    "    else:\n",
    "        # Fallback if map_info is not provided\n",
    "        projection = 'Unknown'\n",
    "        x_pixel_start, y_pixel_start = 1.0, 1.0\n",
    "        map_x, map_y = 0.0, 0.0\n",
    "        x_res, y_res = 1.0, 1.0\n",
    "\n",
    "    # Compute Easting, Northing\n",
    "    # Pixel_row and Pixel_col are zero-based. \n",
    "    # According to ENVI conventions:\n",
    "    # Easting = map_x + (pixel_col - (x_pixel_start - 1)) * x_res\n",
    "    # Northing = map_y - (pixel_row - (y_pixel_start - 1)) * y_res\n",
    "    pixel_row = pixel_indices[:, 0]\n",
    "    pixel_col = pixel_indices[:, 1]\n",
    "    Easting = map_x + (pixel_col - (x_pixel_start - 1)) * x_res\n",
    "    Northing = map_y - (pixel_row - (y_pixel_start - 1)) * y_res\n",
    "\n",
    "    # Insert Pixel info and coordinates\n",
    "    df.insert(0, 'Pixel_Col', pixel_col)\n",
    "    df.insert(0, 'Pixel_Row', pixel_row)\n",
    "    df.insert(0, 'Pixel_id', np.arange(len(df)))\n",
    "    df.insert(3, 'Easting', Easting)\n",
    "    df.insert(4, 'Northing', Northing)\n",
    "\n",
    "    # Check we have enough bands\n",
    "    if bands < (original_bands + corrected_bands):\n",
    "        raise ValueError(\n",
    "            f\"Not enough bands. Expected at least {original_bands + corrected_bands} (original+corrected), but got {bands}.\"\n",
    "        )\n",
    "\n",
    "    # Determine Corrected and Resampled flags\n",
    "    remaining_bands = bands - (original_bands + corrected_bands)\n",
    "    corrected_flag = \"Yes\" if corrected_bands > 0 else \"No\"\n",
    "    resampled_flag = \"Yes\" if remaining_bands > 0 else \"No\"\n",
    "\n",
    "    # Metadata columns: Subdirectory, Data_Source, Sensor_Type, Corrected, Resampled\n",
    "    # Insert these at the very front\n",
    "    df.insert(0, 'Resampled', resampled_flag)\n",
    "    df.insert(0, 'Corrected', corrected_flag)\n",
    "    df.insert(0, 'Sensor_Type', 'Hyperspectral')\n",
    "    df.insert(0, 'Data_Source', 'Flight line')\n",
    "    df.insert(0, 'Subdirectory', folder_name if folder_name else 'Unknown')\n",
    "\n",
    "    # Rename bands with wavelengths\n",
    "    band_names = []\n",
    "    # Original bands\n",
    "    if original_wavelengths is not None and len(original_wavelengths) >= original_bands:\n",
    "        for i in range(original_bands):\n",
    "            wl = original_wavelengths[i]\n",
    "            band_names.append(f\"Original_band_{i+1}_wl_{wl}nm\")\n",
    "    else:\n",
    "        for i in range(1, original_bands + 1):\n",
    "            band_names.append(f\"Original_band_{i}\")\n",
    "\n",
    "    # Corrected bands\n",
    "    if corrected_wavelengths is not None and len(corrected_wavelengths) >= corrected_bands:\n",
    "        for i in range(corrected_bands):\n",
    "            wl = corrected_wavelengths[i]\n",
    "            band_names.append(f\"Corrected_band_{i+1}_wl_{wl}nm\")\n",
    "    elif original_wavelengths is not None and len(original_wavelengths) >= corrected_bands:\n",
    "        for i in range(corrected_bands):\n",
    "            wl = original_wavelengths[i]\n",
    "            band_names.append(f\"Corrected_band_{i+1}_wl_{wl}nm\")\n",
    "    else:\n",
    "        for i in range(1, corrected_bands + 1):\n",
    "            band_names.append(f\"Corrected_band_{i}\")\n",
    "\n",
    "    print(f\"[DEBUG] remaining_bands for resampled sensors: {remaining_bands}\")\n",
    "\n",
    "    sensor_bands_assigned = 0\n",
    "    for sensor, details in config.items():\n",
    "        wavelengths = details.get('wavelengths', [])\n",
    "        for i, wl in enumerate(wavelengths, start=1):\n",
    "            if sensor_bands_assigned < remaining_bands:\n",
    "                band_names.append(f\"{sensor}_band_{i}_wl_{wl}nm\")\n",
    "                sensor_bands_assigned += 1\n",
    "            else:\n",
    "                break\n",
    "        if sensor_bands_assigned >= remaining_bands:\n",
    "            break\n",
    "\n",
    "    if sensor_bands_assigned < remaining_bands:\n",
    "        extra = remaining_bands - sensor_bands_assigned\n",
    "        print(f\"[DEBUG] {extra} leftover bands have no matching sensors/wavelengths in JSON. Naming them generically.\")\n",
    "        for i in range(1, extra + 1):\n",
    "            band_names.append(f\"Unassigned_band_{i}\")\n",
    "\n",
    "    # Now we have Pixel_id, Pixel_Row, Pixel_Col, Easting, Northing, and multiple metadata columns.\n",
    "    # Determine how many leading metadata columns we have before bands:\n",
    "    # Currently: Subdirectory, Data_Source, Sensor_Type, Corrected, Resampled, Pixel_id, Pixel_Row, Pixel_Col, Easting, Northing\n",
    "    # That's 10 columns before bands start.\n",
    "    metadata_count = 10\n",
    "\n",
    "    new_columns = list(df.columns[:metadata_count]) + band_names\n",
    "    if len(new_columns) != df.shape[1]:\n",
    "        raise ValueError(\n",
    "            f\"Band naming mismatch: {len(new_columns)} columns assigned vs {df.shape[1]} in df. Check indexing.\"\n",
    "        )\n",
    "\n",
    "    df.columns = new_columns\n",
    "\n",
    "    print(f\"[DEBUG] Final DataFrame shape: {df.shape}\")\n",
    "    print(\"[DEBUG] Columns assigned successfully.\")\n",
    "\n",
    "    return df\n",
    "\n",
    "def clean_data_and_write_to_csv(df, output_csv_path, chunk_size=100000):\n",
    "    \"\"\"\n",
    "    Cleans a large DataFrame by processing it in chunks and then writes it to a CSV file.\n",
    "    \"\"\"\n",
    "    total_rows = df.shape[0]\n",
    "    num_chunks = (total_rows // chunk_size) + (1 if total_rows % chunk_size else 0)\n",
    "\n",
    "    print(f\"Cleaning data and writing to CSV in {num_chunks} chunk(s).\")\n",
    "\n",
    "    first_chunk = True\n",
    "    for i, start_row in enumerate(range(0, total_rows, chunk_size)):\n",
    "        chunk = df.iloc[start_row:start_row + chunk_size].copy()\n",
    "        non_pixel_cols = [col for col in chunk.columns if not col.startswith('Pixel') and \n",
    "                          col not in ['Subdirectory','Data_Source','Sensor_Type','Corrected','Resampled',\n",
    "                                      'Easting','Northing']]\n",
    "\n",
    "        # Replace -9999 values with NaN\n",
    "        chunk[non_pixel_cols] = chunk[non_pixel_cols].apply(\n",
    "            lambda x: np.where(np.isclose(x, -9999, atol=1), np.nan, x)\n",
    "        )\n",
    "\n",
    "        # Drop rows with all NaNs in non-pixel columns (spectral data)\n",
    "        chunk.dropna(subset=non_pixel_cols, how='all', inplace=True)\n",
    "\n",
    "        mode = 'w' if first_chunk else 'a'\n",
    "        header = True if first_chunk else False\n",
    "        chunk.to_csv(output_csv_path, mode=mode, header=header, index=False)\n",
    "\n",
    "        print(f\"Chunk {i+1}/{num_chunks} processed and written.\")\n",
    "        first_chunk = False\n",
    "\n",
    "    print(f\"Data cleaning complete. Output written to: {output_csv_path}\")\n",
    "\n",
    "def control_function(directory):\n",
    "    \"\"\"\n",
    "    Orchestrates the finding, loading, processing of raster files found in a specified directory,\n",
    "    cleans the processed data, and saves it to a CSV file in the same directory.\n",
    "    \"\"\"\n",
    "    raster_paths = find_raster_files(directory)\n",
    "\n",
    "    if not raster_paths:\n",
    "        print(f\"No matching raster files found in {directory}.\")\n",
    "        return\n",
    "\n",
    "    # Assume original file name (without _envi etc.) is the directory name\n",
    "    base_name = os.path.basename(os.path.normpath(directory))\n",
    "    hdr_file = os.path.join(os.path.dirname(directory), base_name + '.hdr')\n",
    "    if not os.path.isfile(hdr_file):\n",
    "        hdr_file = os.path.join(directory, base_name + '.hdr')\n",
    "\n",
    "    original_wavelengths = None\n",
    "    map_info = None\n",
    "    if os.path.isfile(hdr_file):\n",
    "        img = open_image(hdr_file)\n",
    "        original_wavelengths = img.metadata.get('wavelength', [])\n",
    "        # Convert to float if they are strings\n",
    "        original_wavelengths = [float(w) for w in original_wavelengths]\n",
    "        map_info = img.metadata.get('map info', None)\n",
    "    else:\n",
    "        print(f\"No HDR file found at {hdr_file}. Will use generic band names and no geocoords.\")\n",
    "\n",
    "    corrected_wavelengths = original_wavelengths\n",
    "\n",
    "    # Load and combine raster data\n",
    "    combined_array = load_and_combine_rasters(raster_paths)  \n",
    "    print(f\"Combined array shape for directory {directory}: {combined_array.shape}\")\n",
    "\n",
    "    # Attempt to process and flatten the array into a DataFrame\n",
    "    try:\n",
    "        df_processed = process_and_flatten_array(\n",
    "            combined_array,\n",
    "            json_dir='Resampling',\n",
    "            original_bands=426,\n",
    "            corrected_bands=426,\n",
    "            original_wavelengths=original_wavelengths,\n",
    "            corrected_wavelengths=corrected_wavelengths,\n",
    "            folder_name=base_name,\n",
    "            map_info=map_info\n",
    "        )  \n",
    "        print(f\"DataFrame shape after flattening for directory {directory}: {df_processed.shape}\")\n",
    "    except ValueError as e:\n",
    "        print(f\"ValueError encountered during processing of {directory}: {e}\")\n",
    "        print(\"Check the number of bands vs. the expected column names in process_and_flatten_array().\")\n",
    "        return\n",
    "    except Exception as e:\n",
    "        print(f\"An unexpected error occurred while processing {directory}: {e}\")\n",
    "        return\n",
    "\n",
    "    # Extract the folder name from the directory path\n",
    "    folder_name = os.path.basename(os.path.normpath(directory))\n",
    "    output_csv_name = f\"{folder_name}_spectral_data_all_sensors.csv\"\n",
    "    output_csv_path = os.path.join(directory, output_csv_name)\n",
    "\n",
    "    # Always overwrite if CSV exists\n",
    "    if os.path.exists(output_csv_path):\n",
    "        print(f\"CSV {output_csv_path} already exists and will be overwritten.\")\n",
    "\n",
    "    # Clean data and write to CSV\n",
    "    clean_data_and_write_to_csv(df_processed, output_csv_path)  \n",
    "    print(f\"Processed and cleaned data saved to {output_csv_path}\")\n",
    "\n",
    "# ----- Masking Function Integration -----\n",
    "\n",
    "def mask_raster_with_polygons(\n",
    "    envi_path,\n",
    "    geojson_path,\n",
    "    raster_crs_override=None,\n",
    "    polygons_crs_override=None,\n",
    "    output_masked_suffix=\"_masked\",\n",
    "    plot_output=False,  # Disable plotting by default when processing multiple files\n",
    "    plot_filename=None,  # Not used when plot_output is False\n",
    "    dpi=300\n",
    "):\n",
    "    \"\"\"\n",
    "    Masks an ENVI raster using polygons from a GeoJSON file.\n",
    "\n",
    "    Parameters:\n",
    "    -----------\n",
    "    envi_path : str\n",
    "        Path to the ENVI raster file.\n",
    "    geojson_path : str\n",
    "        Path to the GeoJSON file containing polygons.\n",
    "    raster_crs_override : str, optional\n",
    "        CRS to assign to the raster if it's undefined (e.g., 'EPSG:4326').\n",
    "    polygons_crs_override : str or CRS, optional\n",
    "        CRS to assign to the polygons if they're undefined.\n",
    "    output_masked_suffix : str, optional\n",
    "        Suffix to append to the masked raster filename.\n",
    "    plot_output : bool, optional\n",
    "        Whether to generate and save plots of the results.\n",
    "    plot_filename : str, optional\n",
    "        Filename for the saved plot.\n",
    "    dpi : int, optional\n",
    "        Resolution of the saved plot in dots per inch.\n",
    "\n",
    "    Raises:\n",
    "    -------\n",
    "    FileNotFoundError\n",
    "        If the ENVI raster or GeoJSON file does not exist.\n",
    "    ValueError\n",
    "        If CRS assignments are needed but not provided.\n",
    "\n",
    "    Returns:\n",
    "    --------\n",
    "    masked_raster_path : str\n",
    "        Path to the saved masked raster file.\n",
    "    \"\"\"\n",
    "    def load_data(envi_path, geojson_path):\n",
    "        \"\"\"\n",
    "        Loads the ENVI raster and GeoJSON polygons.\n",
    "        \"\"\"\n",
    "        if not os.path.exists(envi_path):\n",
    "            raise FileNotFoundError(f\"ENVI raster file not found at: {envi_path}\")\n",
    "        if not os.path.exists(geojson_path):\n",
    "            raise FileNotFoundError(f\"GeoJSON file not found at: {geojson_path}\")\n",
    "\n",
    "        raster = rasterio.open(envi_path)\n",
    "        polygons = gpd.read_file(geojson_path)\n",
    "        return raster, polygons\n",
    "\n",
    "    def assign_crs(raster, polygons, raster_crs_override=None, polygons_crs_override=None):\n",
    "        \"\"\"\n",
    "        Determines the CRS for raster and polygons without modifying the raster object.\n",
    "        Returns the raster CRS and the aligned polygons.\n",
    "        \"\"\"\n",
    "        # Handle raster CRS\n",
    "        if raster.crs is None:\n",
    "            if raster_crs_override is not None:\n",
    "                raster_crs = CRS.from_string(raster_crs_override)\n",
    "                print(f\"Assigned CRS {raster_crs_override} to raster.\")\n",
    "            else:\n",
    "                raise ValueError(\"Raster CRS is undefined and no override provided.\")\n",
    "        else:\n",
    "            raster_crs = raster.crs\n",
    "\n",
    "        # Handle polygons CRS\n",
    "        if polygons.crs is None:\n",
    "            if polygons_crs_override is not None:\n",
    "                polygons = polygons.set_crs(polygons_crs_override)\n",
    "                print(f\"Assigned CRS {polygons_crs_override} to polygons.\")\n",
    "            else:\n",
    "                raise ValueError(\"Polygons CRS is undefined and no override provided.\")\n",
    "\n",
    "        return raster_crs, polygons\n",
    "\n",
    "    def align_crs(raster_crs, polygons):\n",
    "        \"\"\"\n",
    "        Reprojects polygons to match raster CRS.\n",
    "        \"\"\"\n",
    "        print(\"Raster CRS:\", raster_crs)\n",
    "        print(\"Polygons CRS:\", polygons.crs)\n",
    "\n",
    "        if polygons.crs != raster_crs:\n",
    "            print(\"Reprojecting polygons to match raster CRS...\")\n",
    "            polygons_aligned = polygons.to_crs(raster_crs)\n",
    "            print(\"Reprojection complete.\")\n",
    "        else:\n",
    "            print(\"CRS of raster and polygons already match.\")\n",
    "            polygons_aligned = polygons\n",
    "\n",
    "        return polygons_aligned\n",
    "\n",
    "    def clip_polygons(raster, polygons_aligned):\n",
    "        \"\"\"\n",
    "        Clips polygons to raster bounds.\n",
    "        \"\"\"\n",
    "        print(\"Clipping polygons to raster bounds...\")\n",
    "        raster_bounds_geom = gpd.GeoDataFrame({'geometry': [box(*raster.bounds)]}, crs=raster.crs)\n",
    "        clipped_polygons = gpd.overlay(polygons_aligned, raster_bounds_geom, how='intersection')\n",
    "        if clipped_polygons.empty:\n",
    "            print(\"No polygons overlap the raster extent.\")\n",
    "        else:\n",
    "            print(f\"Number of Clipped Polygons: {len(clipped_polygons)}\")\n",
    "            print(\"Clipped Polygons Bounds:\", clipped_polygons.total_bounds)\n",
    "        return clipped_polygons\n",
    "\n",
    "    def create_mask(raster, polygons):\n",
    "        \"\"\"\n",
    "        Creates a mask where pixels inside polygons are True and outside are False.\n",
    "        \"\"\"\n",
    "        print(\"Creating mask from polygons...\")\n",
    "        mask = rasterize(\n",
    "            [(geom, 1) for geom in polygons.geometry],\n",
    "            out_shape=(raster.height, raster.width),\n",
    "            transform=raster.transform,\n",
    "            fill=0,  # Background value\n",
    "            dtype='uint8',\n",
    "            all_touched=True  # Include all touched pixels\n",
    "        )\n",
    "        mask = mask.astype(bool)\n",
    "        print(f\"Mask created with shape {mask.shape}. Inside pixels: {np.sum(mask)}\")\n",
    "\n",
    "        # Additional debug: Check unique values\n",
    "        unique_values = np.unique(mask)\n",
    "        print(f\"Unique values in mask: {unique_values}\")\n",
    "\n",
    "        return mask\n",
    "\n",
    "    def apply_mask(raster, mask):\n",
    "        \"\"\"\n",
    "        Applies the mask to raster data, setting areas outside polygons to nodata.\n",
    "        \"\"\"\n",
    "        raster_data = raster.read()\n",
    "        nodata_value = raster.nodata if raster.nodata is not None else -9999\n",
    "        print(f\"Using nodata value: {nodata_value}\")\n",
    "\n",
    "        if raster.count > 1:  # For multi-band rasters\n",
    "            if mask.shape != raster.read(1).shape:\n",
    "                raise ValueError(\"Mask shape does not match raster band shape.\")\n",
    "            mask_expanded = np.repeat(mask[np.newaxis, :, :], raster.count, axis=0)\n",
    "            masked_data = np.where(mask_expanded, raster_data, nodata_value)\n",
    "        else:  # For single-band rasters\n",
    "            masked_data = np.where(mask, raster_data[0], nodata_value)\n",
    "\n",
    "        print(f\"Masked data shape: {masked_data.shape}\")\n",
    "        return masked_data\n",
    "\n",
    "    def save_masked_raster(envi_path, masked_data, nodata, raster, suffix):\n",
    "        \"\"\"\n",
    "        Saves the masked raster to a new file.\n",
    "        \"\"\"\n",
    "        base_name, ext = os.path.splitext(envi_path)\n",
    "        output_file = f\"{base_name}{suffix}{ext}\"\n",
    "        meta = raster.meta.copy()\n",
    "        meta.update({\n",
    "            'dtype': masked_data.dtype,\n",
    "            'nodata': nodata,\n",
    "            'count': raster.count if raster.count > 1 else 1\n",
    "        })\n",
    "\n",
    "        with rasterio.open(output_file, 'w', **meta) as dst:\n",
    "            if raster.count > 1:\n",
    "                for i in range(raster.count):\n",
    "                    dst.write(masked_data[i], i + 1)\n",
    "            else:\n",
    "                dst.write(masked_data, 1)\n",
    "        print(f\"Masked raster saved to: {output_file}\")\n",
    "        return output_file\n",
    "\n",
    "    def plot_results(raster, masked_data, nodata, polygons_aligned, clipped_polygons, plot_path):\n",
    "        \"\"\"\n",
    "        Plots the original raster, clipped polygons, clipped polygons on raster, and masked raster.\n",
    "        Saves the result as a high-resolution PNG and displays the plot.\n",
    "        \"\"\"\n",
    "        print(f\"Masked data stats: Min={masked_data.min()}, Max={masked_data.max()}, Unique={np.unique(masked_data)}\")\n",
    "\n",
    "        fig, axes = plt.subplots(2, 2, figsize=(18, 12))\n",
    "\n",
    "        # Original Raster\n",
    "        original = raster.read(1)\n",
    "        original = np.ma.masked_equal(original, nodata)\n",
    "        axes[0, 0].imshow(\n",
    "            original,\n",
    "            cmap='gray',\n",
    "            extent=(\n",
    "                raster.bounds.left,\n",
    "                raster.bounds.right,\n",
    "                raster.bounds.bottom,\n",
    "                raster.bounds.top\n",
    "            )\n",
    "        )\n",
    "        axes[0, 0].set_title('Original Raster')\n",
    "        axes[0, 0].axis('off')\n",
    "\n",
    "        # Clipped Polygons\n",
    "        clipped_polygons.plot(ax=axes[0, 1], facecolor='none', edgecolor='red')\n",
    "        axes[0, 1].set_title('Clipped Polygons')\n",
    "        axes[0, 1].set_xlim(raster.bounds.left, raster.bounds.right)\n",
    "        axes[0, 1].set_ylim(raster.bounds.bottom, raster.bounds.top)\n",
    "        axes[0, 1].axis('off')\n",
    "\n",
    "        # Clipped Polygons on Raster\n",
    "        axes[1, 0].imshow(\n",
    "            original,\n",
    "            cmap='gray',\n",
    "            extent=(\n",
    "                raster.bounds.left,\n",
    "                raster.bounds.right,\n",
    "                raster.bounds.bottom,\n",
    "                raster.bounds.top\n",
    "            )\n",
    "        )\n",
    "        clipped_polygons.plot(ax=axes[1, 0], facecolor='none', edgecolor='blue', linewidth=0.5)\n",
    "        axes[1, 0].set_title('Clipped Polygons on Raster')\n",
    "        axes[1, 0].set_xlim(raster.bounds.left, raster.bounds.right)\n",
    "        axes[1, 0].set_ylim(raster.bounds.bottom, raster.bounds.top)\n",
    "        axes[1, 0].axis('off')\n",
    "\n",
    "        # Masked Raster\n",
    "        # Find the first valid band to plot\n",
    "        valid_band_found = False\n",
    "        for band_index in range(masked_data.shape[0]):\n",
    "            band = np.ma.masked_equal(masked_data[band_index], nodata)\n",
    "            if np.any(band):  # Check if the band has non-nodata values\n",
    "                print(f\"Plotting band {band_index + 1} (valid band found).\")\n",
    "                valid_band_found = True\n",
    "                axes[1, 1].imshow(\n",
    "                    band,\n",
    "                    cmap='Reds',\n",
    "                    extent=(\n",
    "                        raster.bounds.left,\n",
    "                        raster.bounds.right,\n",
    "                        raster.bounds.bottom,\n",
    "                        raster.bounds.top\n",
    "                    )\n",
    "                )\n",
    "                axes[1, 1].set_title(f'Masked Raster (Band {band_index + 1})')\n",
    "                axes[1, 1].axis('off')\n",
    "                break\n",
    "\n",
    "        if not valid_band_found:\n",
    "            print(\"No valid band found to plot.\")\n",
    "            axes[1, 1].text(\n",
    "                0.5,\n",
    "                0.5,\n",
    "                \"No valid band found\",\n",
    "                ha='center',\n",
    "                va='center',\n",
    "                transform=axes[1, 1].transAxes\n",
    "            )\n",
    "\n",
    "        plt.tight_layout()\n",
    "\n",
    "        # Save the figure\n",
    "        print(f\"Saving plot to {plot_path}\")\n",
    "        plt.savefig(plot_path, dpi=dpi)  # Save as high-resolution PNG\n",
    "\n",
    "        if plot_output:\n",
    "            # Display the plot\n",
    "            plt.show()\n",
    "        else:\n",
    "            plt.close()\n",
    "\n",
    "    # Start of the masking function logic\n",
    "    try:\n",
    "        raster, polygons = load_data(envi_path, geojson_path)\n",
    "    except FileNotFoundError as e:\n",
    "        print(e)\n",
    "        raise\n",
    "\n",
    "    try:\n",
    "        raster_crs, polygons = assign_crs(\n",
    "            raster,\n",
    "            polygons,\n",
    "            raster_crs_override=raster_crs_override,\n",
    "            polygons_crs_override=polygons_crs_override\n",
    "        )\n",
    "    except ValueError as e:\n",
    "        print(e)\n",
    "        raster.close()\n",
    "        raise\n",
    "\n",
    "    polygons_aligned = align_crs(raster_crs, polygons)\n",
    "\n",
    "    clipped_polygons = clip_polygons(raster, polygons_aligned)\n",
    "\n",
    "    if clipped_polygons.empty:\n",
    "        print(\"No polygons overlap the raster extent. Skipping masking.\")\n",
    "        raster.close()\n",
    "        return None\n",
    "\n",
    "    # Check if raster has a valid transform\n",
    "    if not raster.transform.is_identity:\n",
    "        print(\"Raster has a valid geotransform.\")\n",
    "    else:\n",
    "        print(\"Raster has an identity transform. Geospatial operations may be invalid.\")\n",
    "        # Depending on your data, you might want to skip or handle differently\n",
    "        # For now, we'll proceed but be aware that spatial alignment may be incorrect\n",
    "\n",
    "    mask = create_mask(raster, clipped_polygons)\n",
    "\n",
    "    masked_data = apply_mask(raster, mask)\n",
    "\n",
    "    # Handle rasters with no geotransform by informing the user\n",
    "    if raster.transform.is_identity:\n",
    "        print(\"Warning: Raster has an identity transform. Masked data may not be georeferenced correctly.\")\n",
    "\n",
    "    masked_raster_path = save_masked_raster(\n",
    "        envi_path,\n",
    "        masked_data,\n",
    "        raster.nodata,\n",
    "        raster,\n",
    "        suffix=output_masked_suffix\n",
    "    )\n",
    "\n",
    "    # Plot results if enabled\n",
    "    if plot_output and plot_filename:\n",
    "        plot_results(\n",
    "            raster,\n",
    "            masked_data,\n",
    "            raster.nodata,\n",
    "            polygons_aligned,\n",
    "            clipped_polygons,\n",
    "            plot_path=plot_filename\n",
    "        )\n",
    "\n",
    "    # Close the raster dataset\n",
    "    raster.close()\n",
    "\n",
    "    return masked_raster_path\n",
    "\n",
    "# ----- New Function to Process All ENVI Files in a Directory -----\n",
    "\n",
    "def process_all_envi_files_with_mask(\n",
    "    directory,\n",
    "    geojson_path,\n",
    "    raster_crs_override=None,\n",
    "    polygons_crs_override=None,\n",
    "    output_masked_suffix=\"_masked\",\n",
    "    plot_output=False,  # Disable plotting by default when processing multiple files\n",
    "    plot_filename_template=\"plot_results_{envi_basename}.png\",\n",
    "    dpi=300\n",
    "):\n",
    "    \"\"\"\n",
    "    Processes all ENVI raster files in the specified directory by masking them with the provided GeoJSON polygons.\n",
    "\n",
    "    Parameters:\n",
    "    -----------\n",
    "    directory : str\n",
    "        Path to the directory containing ENVI raster files.\n",
    "    geojson_path : str\n",
    "        Path to the GeoJSON file containing polygons for masking.\n",
    "    raster_crs_override : str, optional\n",
    "        CRS to assign to the raster if it's undefined (e.g., 'EPSG:4326').\n",
    "    polygons_crs_override : str or CRS, optional\n",
    "        CRS to assign to the polygons if they're undefined.\n",
    "    output_masked_suffix : str, optional\n",
    "        Suffix to append to each masked raster filename.\n",
    "    plot_output : bool, optional\n",
    "        Whether to generate and save plots for each masked raster.\n",
    "    plot_filename_template : str, optional\n",
    "        Template for the plot filenames. Use '{envi_basename}' as a placeholder for the ENVI file basename.\n",
    "    dpi : int, optional\n",
    "        Resolution of the saved plots in dots per inch.\n",
    "\n",
    "    Returns:\n",
    "    --------\n",
    "    masked_raster_paths : list\n",
    "        List of paths to the saved masked raster files.\n",
    "    \"\"\"\n",
    "    # Find all relevant ENVI raster files in the directory\n",
    "    envi_files = find_raster_files(directory)\n",
    "    print(f\"Found {len(envi_files)} ENVI raster files in {directory}.\")\n",
    "\n",
    "    if not envi_files:\n",
    "        print(\"No ENVI raster files to process.\")\n",
    "        return []\n",
    "\n",
    "    masked_raster_paths = []\n",
    "\n",
    "    for envi_path in envi_files:\n",
    "        try:\n",
    "            print(f\"\\nProcessing ENVI file: {envi_path}\")\n",
    "\n",
    "            # Determine plot filename if plotting is enabled\n",
    "            if plot_output:\n",
    "                envi_basename = os.path.splitext(os.path.basename(envi_path))[0]\n",
    "                plot_filename = plot_filename_template.format(envi_basename=envi_basename)\n",
    "            else:\n",
    "                plot_filename = None\n",
    "\n",
    "            # Apply the masking function\n",
    "            masked_raster = mask_raster_with_polygons(\n",
    "                envi_path=envi_path,\n",
    "                geojson_path=geojson_path,\n",
    "                raster_crs_override=raster_crs_override,\n",
    "                polygons_crs_override=polygons_crs_override,\n",
    "                output_masked_suffix=output_masked_suffix,\n",
    "                plot_output=plot_output,\n",
    "                plot_filename=plot_filename,\n",
    "                dpi=dpi\n",
    "            )\n",
    "\n",
    "            if masked_raster:\n",
    "                masked_raster_paths.append(masked_raster)\n",
    "                print(f\"Successfully masked raster: {masked_raster}\")\n",
    "            else:\n",
    "                print(f\"Masking skipped for raster: {envi_path}\")\n",
    "\n",
    "        except Exception as e:\n",
    "            print(f\"Error processing {envi_path}: {e}\")\n",
    "            continue  # Continue with the next file\n",
    "\n",
    "    print(f\"\\nAll processing complete. {len(masked_raster_paths)} rasters masked successfully.\")\n",
    "    return masked_raster_paths\n",
    "\n",
    "# ----- Simplified Function: Only Folder Name and Polygon Layer as Required Parameters -----\n",
    "\n",
    "def process_folder_with_polygons(\n",
    "    folder_name,\n",
    "    polygon_layer,\n",
    "    raster_crs_override='EPSG:4326',\n",
    "    polygons_crs_override='EPSG:4326',\n",
    "    output_masked_suffix=\"_masked\",\n",
    "    plot_output=False,\n",
    "    plot_filename_template=\"plot_results_{envi_basename}.png\",\n",
    "    dpi=300\n",
    "):\n",
    "    \"\"\"\n",
    "    Simplified function to process all ENVI raster files in a folder by masking them with a polygon layer.\n",
    "    Only the folder name and polygon layer are required; all other parameters are optional with default values.\n",
    "\n",
    "    Parameters:\n",
    "    -----------\n",
    "    folder_name : str\n",
    "        Path to the directory containing ENVI raster files.\n",
    "    polygon_layer : str\n",
    "        Path to the GeoJSON file containing polygons for masking.\n",
    "    raster_crs_override : str, optional\n",
    "        CRS to assign to the raster if it's undefined (default: 'EPSG:4326').\n",
    "    polygons_crs_override : str or CRS, optional\n",
    "        CRS to assign to the polygons if they're undefined (default: 'EPSG:4326').\n",
    "    output_masked_suffix : str, optional\n",
    "        Suffix to append to each masked raster filename (default: \"_masked\").\n",
    "    plot_output : bool, optional\n",
    "        Whether to generate and save plots for each masked raster (default: False).\n",
    "    plot_filename_template : str, optional\n",
    "        Template for the plot filenames. Use '{envi_basename}' as a placeholder (default: \"plot_results_{envi_basename}.png\").\n",
    "    dpi : int, optional\n",
    "        Resolution of the saved plots in dots per inch (default: 300).\n",
    "\n",
    "    Returns:\n",
    "    --------\n",
    "    masked_raster_paths : list\n",
    "        List of paths to the saved masked raster files.\n",
    "    \"\"\"\n",
    "    return process_all_envi_files_with_mask(\n",
    "        directory=folder_name,\n",
    "        geojson_path=polygon_layer,\n",
    "        raster_crs_override=raster_crs_override,\n",
    "        polygons_crs_override=polygons_crs_override,\n",
    "        output_masked_suffix=output_masked_suffix,\n",
    "        plot_output=plot_output,\n",
    "        plot_filename_template=plot_filename_template,\n",
    "        dpi=dpi\n",
    "    )\n",
    "\n",
    "pass\n",
    "\n",
    "import os\n",
    "\n",
    "def process_base_folder(base_folder, polygon_layer, **kwargs):\n",
    "    \"\"\"\n",
    "    Applies a processing function to all subdirectories within the base folder.\n",
    "\n",
    "    Parameters:\n",
    "    -----------\n",
    "    base_folder : str\n",
    "        The path to the base folder containing subdirectories.\n",
    "    polygon_layer : str\n",
    "        Path to the GeoJSON file containing polygons for masking.\n",
    "    **kwargs : dict\n",
    "        Additional keyword arguments to pass to `process_folder_with_polygons`.\n",
    "    \"\"\"\n",
    "    # List all subdirectories\n",
    "    subdirectories = [\n",
    "        os.path.join(base_folder, sub_dir)\n",
    "        for sub_dir in os.listdir(base_folder)\n",
    "        if os.path.isdir(os.path.join(base_folder, sub_dir))\n",
    "    ]\n",
    "\n",
    "    print(f\"Found {len(subdirectories)} subdirectories in {base_folder}.\")\n",
    "\n",
    "    for sub_dir in subdirectories:\n",
    "        print(f\"\\nProcessing subdirectory: {sub_dir}\")\n",
    "        try:\n",
    "            process_folder_with_polygons(\n",
    "                folder_name=sub_dir,\n",
    "                polygon_layer=polygon_layer,\n",
    "                **kwargs  # Pass additional parameters to the processing function\n",
    "            )\n",
    "        except Exception as e:\n",
    "            print(f\"Error processing subdirectory {sub_dir}: {e}\")\n",
    "            continue  # Proceed to the next subdirectory\n",
    "\n",
    "    print(\"All subdirectories processed.\")\n",
    "\n",
    "# Example Usage:\n",
    "# Replace 'base_folder' with the path to your base folder containing subdirectories,\n",
    "# and 'polygon_layer' with the path to your GeoJSON polygon file.\n",
    "# Additional optional parameters can be passed as keyword arguments.\n",
    "\n",
    "pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ebcf170b-9e9c-4933-aec7-09ad87dd7ea7",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found 9 ENVI raster files in non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance.\n",
      "\n",
      "Processing ENVI file: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance\n",
      "Raster CRS: EPSG:32610\n",
      "Polygons CRS: EPSG:4326\n",
      "Reprojecting polygons to match raster CRS...\n",
      "Reprojection complete.\n",
      "Clipping polygons to raster bounds...\n",
      "Number of Clipped Polygons: 88\n",
      "Clipped Polygons Bounds: [ 577927.         5076096.65771862  578849.73218554 5077086.72333448]\n",
      "Raster has a valid geotransform.\n",
      "Creating mask from polygons...\n",
      "Mask created with shape (14396, 1098). Inside pixels: 1206\n",
      "Unique values in mask: [False  True]\n",
      "Using nodata value: -9999.0\n",
      "Masked data shape: (426, 14396, 1098)\n",
      "Masked raster saved to: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance_masked\n",
      "Successfully masked raster: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance_masked\n",
      "\n",
      "Processing ENVI file: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi\n",
      "Raster CRS: EPSG:32610\n",
      "Polygons CRS: EPSG:4326\n",
      "Reprojecting polygons to match raster CRS...\n",
      "Reprojection complete.\n",
      "Clipping polygons to raster bounds...\n",
      "Number of Clipped Polygons: 88\n",
      "Clipped Polygons Bounds: [ 577927.         5076096.65771862  578849.73218554 5077086.72333448]\n",
      "Raster has a valid geotransform.\n",
      "Creating mask from polygons...\n",
      "Mask created with shape (14396, 1098). Inside pixels: 1206\n",
      "Unique values in mask: [False  True]\n",
      "Using nodata value: -9999.0\n",
      "Masked data shape: (426, 14396, 1098)\n",
      "Masked raster saved to: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_masked\n",
      "Successfully masked raster: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_masked\n",
      "\n",
      "Processing ENVI file: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_Landsat_5_TM.img\n",
      "Raster CRS: EPSG:32610\n",
      "Polygons CRS: EPSG:4326\n",
      "Reprojecting polygons to match raster CRS...\n",
      "Reprojection complete.\n",
      "Clipping polygons to raster bounds...\n",
      "Number of Clipped Polygons: 88\n",
      "Clipped Polygons Bounds: [ 577927.         5076096.65771862  578849.73218554 5077086.72333448]\n",
      "Raster has a valid geotransform.\n",
      "Creating mask from polygons...\n",
      "Mask created with shape (14396, 1098). Inside pixels: 1206\n",
      "Unique values in mask: [False  True]\n",
      "Using nodata value: -9999\n",
      "Masked data shape: (6, 14396, 1098)\n",
      "Masked raster saved to: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_Landsat_5_TM_masked.img\n",
      "Successfully masked raster: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_Landsat_5_TM_masked.img\n",
      "\n",
      "Processing ENVI file: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_Landsat_7_ETMplus.img\n",
      "Raster CRS: EPSG:32610\n",
      "Polygons CRS: EPSG:4326\n",
      "Reprojecting polygons to match raster CRS...\n",
      "Reprojection complete.\n",
      "Clipping polygons to raster bounds...\n",
      "Number of Clipped Polygons: 88\n",
      "Clipped Polygons Bounds: [ 577927.         5076096.65771862  578849.73218554 5077086.72333448]\n",
      "Raster has a valid geotransform.\n",
      "Creating mask from polygons...\n",
      "Mask created with shape (14396, 1098). Inside pixels: 1206\n",
      "Unique values in mask: [False  True]\n",
      "Using nodata value: -9999\n",
      "Masked data shape: (6, 14396, 1098)\n",
      "Masked raster saved to: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_Landsat_7_ETMplus_masked.img\n",
      "Successfully masked raster: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_Landsat_7_ETMplus_masked.img\n",
      "\n",
      "Processing ENVI file: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_Landsat_8_OLI.img\n",
      "Raster CRS: EPSG:32610\n",
      "Polygons CRS: EPSG:4326\n",
      "Reprojecting polygons to match raster CRS...\n",
      "Reprojection complete.\n",
      "Clipping polygons to raster bounds...\n",
      "Number of Clipped Polygons: 88\n",
      "Clipped Polygons Bounds: [ 577927.         5076096.65771862  578849.73218554 5077086.72333448]\n",
      "Raster has a valid geotransform.\n",
      "Creating mask from polygons...\n",
      "Mask created with shape (14396, 1098). Inside pixels: 1206\n",
      "Unique values in mask: [False  True]\n",
      "Using nodata value: -9999\n",
      "Masked data shape: (7, 14396, 1098)\n",
      "Masked raster saved to: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_Landsat_8_OLI_masked.img\n",
      "Successfully masked raster: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_Landsat_8_OLI_masked.img\n",
      "\n",
      "Processing ENVI file: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_Landsat_9_OLI-2.img\n",
      "Raster CRS: EPSG:32610\n",
      "Polygons CRS: EPSG:4326\n",
      "Reprojecting polygons to match raster CRS...\n",
      "Reprojection complete.\n",
      "Clipping polygons to raster bounds...\n",
      "Number of Clipped Polygons: 88\n",
      "Clipped Polygons Bounds: [ 577927.         5076096.65771862  578849.73218554 5077086.72333448]\n",
      "Raster has a valid geotransform.\n",
      "Creating mask from polygons...\n",
      "Mask created with shape (14396, 1098). Inside pixels: 1206\n",
      "Unique values in mask: [False  True]\n",
      "Using nodata value: -9999\n",
      "Masked data shape: (7, 14396, 1098)\n",
      "Masked raster saved to: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_Landsat_9_OLI-2_masked.img\n",
      "Successfully masked raster: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_Landsat_9_OLI-2_masked.img\n",
      "\n",
      "Processing ENVI file: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_MicaSense-to-match_OLI_and_OLI-2.img\n",
      "Raster CRS: EPSG:32610\n",
      "Polygons CRS: EPSG:4326\n",
      "Reprojecting polygons to match raster CRS...\n",
      "Reprojection complete.\n",
      "Clipping polygons to raster bounds...\n",
      "Number of Clipped Polygons: 88\n",
      "Clipped Polygons Bounds: [ 577927.         5076096.65771862  578849.73218554 5077086.72333448]\n",
      "Raster has a valid geotransform.\n",
      "Creating mask from polygons...\n",
      "Mask created with shape (14396, 1098). Inside pixels: 1206\n",
      "Unique values in mask: [False  True]\n",
      "Using nodata value: -9999\n",
      "Masked data shape: (5, 14396, 1098)\n",
      "Masked raster saved to: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_MicaSense-to-match_OLI_and_OLI-2_masked.img\n",
      "Successfully masked raster: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_MicaSense-to-match_OLI_and_OLI-2_masked.img\n",
      "\n",
      "Processing ENVI file: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_MicaSense-to-match_TM_and_ETMplus.img\n",
      "Raster CRS: EPSG:32610\n",
      "Polygons CRS: EPSG:4326\n",
      "Reprojecting polygons to match raster CRS...\n",
      "Reprojection complete.\n",
      "Clipping polygons to raster bounds...\n",
      "Number of Clipped Polygons: 88\n",
      "Clipped Polygons Bounds: [ 577927.         5076096.65771862  578849.73218554 5077086.72333448]\n",
      "Raster has a valid geotransform.\n",
      "Creating mask from polygons...\n",
      "Mask created with shape (14396, 1098). Inside pixels: 1206\n",
      "Unique values in mask: [False  True]\n",
      "Using nodata value: -9999\n",
      "Masked data shape: (4, 14396, 1098)\n",
      "Masked raster saved to: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_MicaSense-to-match_TM_and_ETMplus_masked.img\n",
      "Successfully masked raster: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_MicaSense-to-match_TM_and_ETMplus_masked.img\n",
      "\n",
      "Processing ENVI file: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_MicaSense.img\n",
      "Raster CRS: EPSG:32610\n",
      "Polygons CRS: EPSG:4326\n",
      "Reprojecting polygons to match raster CRS...\n",
      "Reprojection complete.\n",
      "Clipping polygons to raster bounds...\n",
      "Number of Clipped Polygons: 88\n",
      "Clipped Polygons Bounds: [ 577927.         5076096.65771862  578849.73218554 5077086.72333448]\n",
      "Raster has a valid geotransform.\n",
      "Creating mask from polygons...\n",
      "Mask created with shape (14396, 1098). Inside pixels: 1206\n",
      "Unique values in mask: [False  True]\n",
      "Using nodata value: -9999\n",
      "Masked data shape: (10, 14396, 1098)\n",
      "Masked raster saved to: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_MicaSense_masked.img\n",
      "Successfully masked raster: non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_MicaSense_masked.img\n",
      "\n",
      "All processing complete. 9 rasters masked successfully.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "['non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance_masked',\n",
       " 'non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_masked',\n",
       " 'non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_Landsat_5_TM_masked.img',\n",
       " 'non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_Landsat_7_ETMplus_masked.img',\n",
       " 'non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_Landsat_8_OLI_masked.img',\n",
       " 'non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_Landsat_9_OLI-2_masked.img',\n",
       " 'non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_MicaSense-to-match_OLI_and_OLI-2_masked.img',\n",
       " 'non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_MicaSense-to-match_TM_and_ETMplus_masked.img',\n",
       " 'non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance__envi_resample_MicaSense_masked.img']"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "process_folder_with_polygons(\n",
    "     # Define the folder containing ENVI files\n",
    "    folder_name = 'non_provisional_flightlines_WREF/NEON_D16_WREF_DP1_L041-1_20230625_directional_reflectance',\n",
    "\n",
    "    # Define the path to the GeoJSON file for masking\n",
    "    polygon_layer = 'aop_macrosystems_data_1_7_25.geojson'\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ca31b9d5-d785-4000-bbfb-3de8602f711d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CRS: None\n",
      "Transform: | 1.00, 0.00, 0.00|\n",
      "| 0.00, 1.00, 0.00|\n",
      "| 0.00, 0.00, 1.00|\n",
      "Bounds: BoundingBox(left=0.0, bottom=7831.0, right=1102.0, top=0.0)\n",
      "Width: 1102, Height: 7831\n",
      "Number of Bands: 6\n",
      "Wavelengths: Not Found\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import rasterio\n",
    "\n",
    "def inspect_raster(raster_path):\n",
    "    with rasterio.open(raster_path) as src:\n",
    "        print(f\"CRS: {src.crs}\")\n",
    "        print(f\"Transform: {src.transform}\")\n",
    "        print(f\"Bounds: {src.bounds}\")\n",
    "        print(f\"Width: {src.width}, Height: {src.height}\")\n",
    "        print(f\"Number of Bands: {src.count}\")\n",
    "        wavelengths = src.tags().get('wavelength', 'Not Found')\n",
    "        print(f\"Wavelengths: {wavelengths}\\n\")\n",
    "\n",
    "# Example usage:\n",
    "resampled_raster = 'non_provisional_flightlines_YELL_2/NEON_D12_YELL_DP1_L005-1_20220626_directional_reflectance/NEON_D12_YELL_DP1_L005-1_20220626_directional_reflectance__envi_resample_Landsat_5_TM.img'\n",
    "inspect_raster(resampled_raster)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "86b30d1d-9460-453b-9f43-0793e1d1d3e3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found 2 subdirectories in full_test.\n",
      "\n",
      "Processing subdirectory: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance\n",
      "Found 9 ENVI raster files in full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance.\n",
      "\n",
      "Processing ENVI file: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance\n",
      "Raster CRS: EPSG:32613\n",
      "Polygons CRS: EPSG:4326\n",
      "Reprojecting polygons to match raster CRS...\n",
      "Reprojection complete.\n",
      "Clipping polygons to raster bounds...\n",
      "Number of Clipped Polygons: 531\n",
      "Clipped Polygons Bounds: [ 457121.94860384 4426034.20421104  457850.         4434663.4225484 ]\n",
      "Raster has a valid geotransform.\n",
      "Creating mask from polygons...\n",
      "Mask created with shape (11031, 839). Inside pixels: 18714\n",
      "Unique values in mask: [False  True]\n",
      "Using nodata value: -9999.0\n",
      "Masked data shape: (426, 11031, 839)\n",
      "Masked raster saved to: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance_processed\n",
      "Successfully masked raster: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance_processed\n",
      "\n",
      "Processing ENVI file: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance__envi\n",
      "Raster CRS: EPSG:32613\n",
      "Polygons CRS: EPSG:4326\n",
      "Reprojecting polygons to match raster CRS...\n",
      "Reprojection complete.\n",
      "Clipping polygons to raster bounds...\n",
      "Number of Clipped Polygons: 531\n",
      "Clipped Polygons Bounds: [ 457121.94860384 4426034.20421104  457850.         4434663.4225484 ]\n",
      "Raster has a valid geotransform.\n",
      "Creating mask from polygons...\n",
      "Mask created with shape (11031, 839). Inside pixels: 18714\n",
      "Unique values in mask: [False  True]\n",
      "Using nodata value: -9999.0\n",
      "Masked data shape: (426, 11031, 839)\n",
      "Masked raster saved to: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance__envi_processed\n",
      "Successfully masked raster: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance__envi_processed\n",
      "\n",
      "Processing ENVI file: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance__envi_resample_Landsat_5_TM.img\n",
      "Raster CRS: EPSG:32613\n",
      "Polygons CRS: EPSG:4326\n",
      "Reprojecting polygons to match raster CRS...\n",
      "Reprojection complete.\n",
      "Clipping polygons to raster bounds...\n",
      "Number of Clipped Polygons: 531\n",
      "Clipped Polygons Bounds: [ 457121.94860384 4426034.20421104  457850.         4434663.4225484 ]\n",
      "Raster has a valid geotransform.\n",
      "Creating mask from polygons...\n",
      "Mask created with shape (11031, 839). Inside pixels: 18714\n",
      "Unique values in mask: [False  True]\n",
      "Using nodata value: -9999\n",
      "Masked data shape: (6, 11031, 839)\n",
      "Masked raster saved to: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance__envi_resample_Landsat_5_TM_processed.img\n",
      "Successfully masked raster: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance__envi_resample_Landsat_5_TM_processed.img\n",
      "\n",
      "Processing ENVI file: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance__envi_resample_Landsat_7_ETMplus.img\n",
      "Raster CRS: EPSG:32613\n",
      "Polygons CRS: EPSG:4326\n",
      "Reprojecting polygons to match raster CRS...\n",
      "Reprojection complete.\n",
      "Clipping polygons to raster bounds...\n",
      "Number of Clipped Polygons: 531\n",
      "Clipped Polygons Bounds: [ 457121.94860384 4426034.20421104  457850.         4434663.4225484 ]\n",
      "Raster has a valid geotransform.\n",
      "Creating mask from polygons...\n",
      "Mask created with shape (11031, 839). Inside pixels: 18714\n",
      "Unique values in mask: [False  True]\n",
      "Using nodata value: -9999\n",
      "Masked data shape: (6, 11031, 839)\n",
      "Masked raster saved to: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance__envi_resample_Landsat_7_ETMplus_processed.img\n",
      "Successfully masked raster: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance__envi_resample_Landsat_7_ETMplus_processed.img\n",
      "\n",
      "Processing ENVI file: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance__envi_resample_Landsat_8_OLI.img\n",
      "Raster CRS: EPSG:32613\n",
      "Polygons CRS: EPSG:4326\n",
      "Reprojecting polygons to match raster CRS...\n",
      "Reprojection complete.\n",
      "Clipping polygons to raster bounds...\n",
      "Number of Clipped Polygons: 531\n",
      "Clipped Polygons Bounds: [ 457121.94860384 4426034.20421104  457850.         4434663.4225484 ]\n",
      "Raster has a valid geotransform.\n",
      "Creating mask from polygons...\n",
      "Mask created with shape (11031, 839). Inside pixels: 18714\n",
      "Unique values in mask: [False  True]\n",
      "Using nodata value: -9999\n",
      "Masked data shape: (7, 11031, 839)\n",
      "Masked raster saved to: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance__envi_resample_Landsat_8_OLI_processed.img\n",
      "Successfully masked raster: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance__envi_resample_Landsat_8_OLI_processed.img\n",
      "\n",
      "Processing ENVI file: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance__envi_resample_Landsat_9_OLI-2.img\n",
      "Raster CRS: EPSG:32613\n",
      "Polygons CRS: EPSG:4326\n",
      "Reprojecting polygons to match raster CRS...\n",
      "Reprojection complete.\n",
      "Clipping polygons to raster bounds...\n",
      "Number of Clipped Polygons: 531\n",
      "Clipped Polygons Bounds: [ 457121.94860384 4426034.20421104  457850.         4434663.4225484 ]\n",
      "Raster has a valid geotransform.\n",
      "Creating mask from polygons...\n",
      "Mask created with shape (11031, 839). Inside pixels: 18714\n",
      "Unique values in mask: [False  True]\n",
      "Using nodata value: -9999\n",
      "Masked data shape: (7, 11031, 839)\n",
      "Masked raster saved to: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance__envi_resample_Landsat_9_OLI-2_processed.img\n",
      "Successfully masked raster: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance__envi_resample_Landsat_9_OLI-2_processed.img\n",
      "\n",
      "Processing ENVI file: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance__envi_resample_MicaSense-to-match_OLI_and_OLI-2.img\n",
      "Raster CRS: EPSG:32613\n",
      "Polygons CRS: EPSG:4326\n",
      "Reprojecting polygons to match raster CRS...\n",
      "Reprojection complete.\n",
      "Clipping polygons to raster bounds...\n",
      "Number of Clipped Polygons: 531\n",
      "Clipped Polygons Bounds: [ 457121.94860384 4426034.20421104  457850.         4434663.4225484 ]\n",
      "Raster has a valid geotransform.\n",
      "Creating mask from polygons...\n",
      "Mask created with shape (11031, 839). Inside pixels: 18714\n",
      "Unique values in mask: [False  True]\n",
      "Using nodata value: -9999\n",
      "Masked data shape: (5, 11031, 839)\n",
      "Masked raster saved to: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance__envi_resample_MicaSense-to-match_OLI_and_OLI-2_processed.img\n",
      "Successfully masked raster: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance__envi_resample_MicaSense-to-match_OLI_and_OLI-2_processed.img\n",
      "\n",
      "Processing ENVI file: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance__envi_resample_MicaSense-to-match_TM_and_ETMplus.img\n",
      "Raster CRS: EPSG:32613\n",
      "Polygons CRS: EPSG:4326\n",
      "Reprojecting polygons to match raster CRS...\n",
      "Reprojection complete.\n",
      "Clipping polygons to raster bounds...\n",
      "Number of Clipped Polygons: 531\n",
      "Clipped Polygons Bounds: [ 457121.94860384 4426034.20421104  457850.         4434663.4225484 ]\n",
      "Raster has a valid geotransform.\n",
      "Creating mask from polygons...\n",
      "Mask created with shape (11031, 839). Inside pixels: 18714\n",
      "Unique values in mask: [False  True]\n",
      "Using nodata value: -9999\n",
      "Masked data shape: (4, 11031, 839)\n",
      "Masked raster saved to: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance__envi_resample_MicaSense-to-match_TM_and_ETMplus_processed.img\n",
      "Successfully masked raster: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance__envi_resample_MicaSense-to-match_TM_and_ETMplus_processed.img\n",
      "\n",
      "Processing ENVI file: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance__envi_resample_MicaSense.img\n",
      "Raster CRS: EPSG:32613\n",
      "Polygons CRS: EPSG:4326\n",
      "Reprojecting polygons to match raster CRS...\n",
      "Reprojection complete.\n",
      "Clipping polygons to raster bounds...\n",
      "Number of Clipped Polygons: 531\n",
      "Clipped Polygons Bounds: [ 457121.94860384 4426034.20421104  457850.         4434663.4225484 ]\n",
      "Raster has a valid geotransform.\n",
      "Creating mask from polygons...\n",
      "Mask created with shape (11031, 839). Inside pixels: 18714\n",
      "Unique values in mask: [False  True]\n",
      "Using nodata value: -9999\n",
      "Masked data shape: (10, 11031, 839)\n",
      "Masked raster saved to: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance__envi_resample_MicaSense_processed.img\n",
      "Successfully masked raster: full_test/NEON_D13_NIWO_DP1_20200731_155646_reflectance/NEON_D13_NIWO_DP1_20200731_155646_reflectance__envi_resample_MicaSense_processed.img\n",
      "\n",
      "All processing complete. 9 rasters masked successfully.\n",
      "\n",
      "Processing subdirectory: full_test/.ipynb_checkpoints\n",
      "Found 0 ENVI raster files in full_test/.ipynb_checkpoints.\n",
      "No ENVI raster files to process.\n",
      "All subdirectories processed.\n"
     ]
    }
   ],
   "source": [
    "# Define paths\n",
    "base_folder = \"full_test\"\n",
    "polygon_layer = \"aop_macrosystems_data_1_7_25.geojson\"\n",
    "\n",
    "# Process all subdirectories\n",
    "process_base_folder(\n",
    "    base_folder=base_folder,\n",
    "    polygon_layer=polygon_layer,\n",
    "    raster_crs_override=\"EPSG:4326\",  # Optional CRS override\n",
    "    polygons_crs_override=\"EPSG:4326\",  # Optional CRS override\n",
    "    output_masked_suffix=\"_processed\",  # Optional suffix for output\n",
    "    plot_output=False,  # Disable plotting\n",
    "    dpi=300  # Set plot resolution\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb665a2c-a31d-4e7e-8f54-e10a17d5946e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "macrosystems",
   "language": "python",
   "name": "macrosystems"
  },
  "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}