Drone_correction_and_pixel_extraction_2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5a993e8c-ea40-4a9c-863d-79470d9c2ff0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: spectral in /opt/conda/envs/macrosystems/lib/python3.10/site-packages (0.23.1)\n",
      "Requirement already satisfied: numpy in /opt/conda/envs/macrosystems/lib/python3.10/site-packages (from spectral) (1.26.3)\n",
      "Note: you may need to restart the kernel to use updated packages.\n"
     ]
    }
   ],
   "source": [
    "%pip install spectral"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "792dfdb9-2f39-41eb-a890-623e66ff3abd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import hytools.masks.cloud\n",
    "\n",
    "# Define a function to override the original `zhai_cloud`\n",
    "def patched_zhai_cloud(hy_obj, cloud, shadow, T1=0.01, t2=.1, t3=.25, t4=.5, T7=9, T8=9):\n",
    "    blue = hy_obj.get_wave(444)   # New wavelength\n",
    "    green = hy_obj.get_wave(560)  # New wavelength\n",
    "    red = hy_obj.get_wave(650)    # New wavelength\n",
    "    nir = hy_obj.get_wave(862)    # New wavelength\n",
    "    \n",
    "    # Call the original function but with modified parameters\n",
    "    return hytools.masks.cloud.zhai_cloud(hy_obj, cloud, shadow, T1, t2, t3, t4, T7, T8)\n",
    "\n",
    "# Patch the function in the `hytools.masks.cloud` module\n",
    "hytools.masks.cloud.zhai_cloud = patched_zhai_cloud\n",
    "\n",
    "# Now, when `zhai_cloud()` is called from `hytools`, it will use the modified version\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2d673681-65e3-40af-8cd0-b6627a9ccc2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import hytools as ht\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import random\n",
    "import nbconvert\n",
    "import time\n",
    "import os\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d5e0f78d-5405-4f8b-94e4-2468a943dc3f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/home/jovyan/data-store/cross-sensor-cal'"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "os.getcwd()"
   ]
  },
  {
   "cell_type": "raw",
   "id": "d916d094-4ea2-4748-af18-c6fd2b53a940",
   "metadata": {},
   "source": [
    "import rasterio\n",
    "from rasterio.transform import Affine\n",
    "import numpy as np\n",
    "\n",
    "def modify_tiff_to_match_envi(tiff_path, output_path):\n",
    "    # Open the original GeoTIFF\n",
    "    with rasterio.open(tiff_path) as src:\n",
    "        # Read metadata and data\n",
    "        meta = src.meta.copy()\n",
    "        data = src.read()  # Read all bands\n",
    "\n",
    "        # Flip the data vertically to match positive y-resolution\n",
    "        flipped_data = np.flip(data, axis=1)  # Flip along the height (rows)\n",
    "\n",
    "        # Update the transform to use a positive y-resolution\n",
    "        transform = src.transform\n",
    "        updated_transform = Affine(transform.a, transform.b, transform.c, transform.d, -transform.e, transform.f)\n",
    "\n",
    "        # Update the metadata with the new transform\n",
    "        meta.update({\"transform\": updated_transform})\n",
    "\n",
    "    # Save the modified GeoTIFF\n",
    "    with rasterio.open(output_path, \"w\", **meta) as dst:\n",
    "        dst.write(flipped_data)\n",
    "\n",
    "    print(f\"Modified GeoTIFF saved to {output_path}\")\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "raw",
   "id": "60eb5a83-9b7a-4db0-9293-2490fecb2191",
   "metadata": {},
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import drone_spectral_unmixing_tools_original as el_spectral\n",
    "\n",
    "#print(dir(el_spectral))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "eba67b7b-da27-4b34-962a-fe20a20caef2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/home/jovyan/data-store/cross-sensor-cal/data_set2/aligned_chm.tif'"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_dir = '/home/jovyan/data-store/cross-sensor-cal/data_set2'\n",
    "file_path = os.path.join(data_dir,'aligned_orthomosaic.tif')\n",
    "slope_path = os.path.join(data_dir,'aligned_slope.tif')\n",
    "aspect_path = os.path.join(data_dir,'aligned_aspect.tif')\n",
    "chm_path = os.path.join(data_dir,'aligned_chm.tif')\n",
    "\n",
    "date_time_str = '2023-08-16 18:53:18'\n",
    "chm_path"
   ]
  },
  {
   "cell_type": "raw",
   "id": "907a4d64-c3ce-4f54-89db-1ba503b67b44",
   "metadata": {},
   "source": [
    "\n",
    "# chnage here based on AOP underflights or not. If underflights, the polygons are at the AOP geojson file. If not within the AOP, the polygons are at uas geojson file\n",
    "# geojson_path = '/home/jovyan/data-store/cross-sensor-cal/uas_polygons_1_24_2025_analysis_ready.geojson'\n",
    "# geojson_path_half = '/home/jovyan/data-store/cross-sensor-cal/uas_polygons_1_24_2025_analysis_ready_half_diam.geojson'\n",
    "\n",
    "geojson_path = '/home/jovyan/data-store/cross-sensor-cal/aop_macrosystems_data_1_7_25.geojson'\n",
    "geojson_path_half = '/home/jovyan/data-store/cross-sensor-cal/aop_polygons_1_24_2025_analysis_ready_half_diam.geojson'\n",
    "\n",
    "\n",
    "output_csv = \"/home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/NEON_D13_NIWO_reflectance_per_pixel.csv\"\n",
    "output_csv_half = \"/home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/NEON_D13_NIWO_reflectance_per_pixel_half.csv\"\n",
    "\n",
    "raster_path = '/home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/NEON_D13_NIWO_aligned_orthomosaic'\n",
    "# main_image_name = os.path.basename(directory)\n",
    "# raster_path = os.path.join(directory, main_image_name + '__reflectance_corrected' )\n",
    "\n",
    "# raster_path_2 = os.path.join(directory, main_image_name + '__reflectance_corrected_flipped' )\n",
    "\n",
    "\n",
    "raster_path"
   ]
  },
  {
   "cell_type": "raw",
   "id": "d539d4ec-e66a-4707-828b-47fd77cdd816",
   "metadata": {},
   "source": [
    "# Input and output paths\n",
    "tiff_path = file_path  # Path to the original GeoTIFF\n",
    "output_path = os.path.join(data_dir,\"modified_to_match_envi.tif\")  # Path to the modified GeoTIFF\n",
    "\n",
    "# Run the function\n",
    "modify_tiff_to_match_envi(tiff_path, output_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "54647211-8400-4fa0-a6a9-e6f9cc8bb9b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from osgeo import gdal, osr\n",
    "import re\n",
    "from typing import Union\n",
    "import h5py\n",
    "import numpy as np\n",
    "import os\n",
    "\n",
    "# Define Wavelengths and Full Width at Half Maximum (FWHM) for spectral bands\n",
    "wavelengths = [444, 475, 531, 560, 650, 668, 705, 717, 740, 862]  # In nm\n",
    "fwhm = [28, 32, 14, 27, 16, 14, 10, 12, 18, 57]\n",
    "\n",
    "def create_h5_file_from_dict(data, h5_file, group_name=\"/\"):\n",
    "    \"\"\"Recursively create groups and datasets in an HDF5 file from a dictionary.\"\"\"\n",
    "    for key, value in data.items():\n",
    "        if isinstance(value, dict):\n",
    "            subgroup_name = f\"{group_name}/{key}\"\n",
    "            _ = h5_file.create_group(subgroup_name)\n",
    "            create_h5_file_from_dict(value, h5_file, subgroup_name)\n",
    "        else:\n",
    "            dataset_name = f\"{group_name}/{key}\"\n",
    "            if os.path.basename(group_name) == 'Coordinate_System':\n",
    "                dataset = h5_file.create_dataset(dataset_name, data=str(value), dtype=h5py.string_dtype())\n",
    "            else:\n",
    "                dataset = h5_file.create_dataset(dataset_name, data=value)\n",
    "            if key in ['Wavelength', 'FWHM']:\n",
    "                dataset.attrs['Units'] = 'nanometers'\n",
    "\n",
    "def remove_non_numeric(input_string):\n",
    "    \"\"\"Remove non-numeric characters from a string.\"\"\"\n",
    "    return re.sub(r'\\D', '', input_string)\n",
    "\n",
    "def load_data(data: Union[str, np.array], expected_shape):\n",
    "    \"\"\"Loads data from file or ensures an array has the expected shape.\"\"\"\n",
    "    if isinstance(data, str):\n",
    "        data_array = gdal.Open(data).ReadAsArray()\n",
    "    else:\n",
    "        data_array = np.array(data)\n",
    "    \n",
    "    if data_array.shape != expected_shape:\n",
    "        raise ValueError(f\"Data shape mismatch: Expected {expected_shape}, but got {data_array.shape}\")\n",
    "    \n",
    "    return data_array\n",
    "\n",
    "def tiff_to_h5(reflectance_tiff_path: str, slope_data: Union[str, np.array],\n",
    "               aspect_data: Union[str, np.array], path_length_data: Union[str, np.array],\n",
    "               solar_zenith_data: Union[str, np.array], solar_azimuth_data: Union[str, np.array]):\n",
    "    \"\"\"Converts a TIFF reflectance file into an HDF5 format with metadata and ancillary data.\"\"\"\n",
    "\n",
    "    # Open reflectance TIFF\n",
    "    reflectance_ds = gdal.Open(reflectance_tiff_path)\n",
    "    reflectance_data = reflectance_ds.ReadAsArray()[:10, :, :]  # Use first 10 bands\n",
    "    img_shape = reflectance_data.shape[1:]  # (rows, cols)\n",
    "       \n",
    "\n",
    "    # Load ancillary data\n",
    "    slope_data = load_data(slope_data, img_shape)\n",
    "    aspect_data = load_data(aspect_data, img_shape)\n",
    "    path_length_data = load_data(path_length_data, img_shape)\n",
    "    solar_zenith_data = load_data(solar_zenith_data, img_shape)\n",
    "    solar_azimuth_data = load_data(solar_azimuth_data, img_shape)\n",
    "\n",
    "    # Set sensor angles to zero matrices matching image shape\n",
    "    sensor_zenith_data = np.zeros(img_shape)\n",
    "    sensor_azimuth_data = np.zeros(img_shape)\n",
    "\n",
    "    # Extract spatial metadata\n",
    "    proj = reflectance_ds.GetProjection()\n",
    "    geo_transform = reflectance_ds.GetGeoTransform()\n",
    "    \n",
    "    spatial_ref = osr.SpatialReference()\n",
    "    spatial_ref.ImportFromWkt(proj)\n",
    "    epsg_code = spatial_ref.GetAuthorityCode(None)\n",
    "\n",
    "    utm_zone = int(remove_non_numeric(proj.split(\"UTM zone \")[1].split(\",\")[0])) if \"UTM zone\" in proj else None\n",
    "    map_info_string = f\"UTM, 1.000, 1.000, {geo_transform[0]:.3f}, {geo_transform[3]:.3f}, {geo_transform[1]:.3f}, {geo_transform[5]:.3f}, {utm_zone}, North, WGS-84, units=Meters, 0\"\n",
    "\n",
    "    # HDF5 Structure\n",
    "    h5_data = {\n",
    "        'NIWO': {\n",
    "            'Reflectance': {\n",
    "                'Metadata': {\n",
    "                    'GeoTransform': geo_transform,  # Add GeoTransform here\n",
    "                    'Coordinate_System': {\n",
    "                        'Coordinate_System_String': proj,\n",
    "                        'EPSG Code': epsg_code,\n",
    "                        'Map_Info': map_info_string,\n",
    "                        'Proj4': spatial_ref.ExportToProj4()\n",
    "                    },\n",
    "                    \"Ancillary_Imagery\": {\n",
    "                        \"Path_Length\": path_length_data,  # Now uses raster-based path length\n",
    "                        \"Slope\": slope_data,\n",
    "                        \"Aspect\": aspect_data\n",
    "                    },\n",
    "                    \"Logs\": {\n",
    "                        \"Solar_Azimuth_Angle\": solar_azimuth_data,\n",
    "                        \"Solar_Zenith_Angle\": solar_zenith_data\n",
    "                    },\n",
    "                    \"to-sensor_Azimuth_Angle\": sensor_azimuth_data,\n",
    "                    \"to-sensor_Zenith_Angle\": sensor_zenith_data,\n",
    "                    \"Spectral_Data\": {\n",
    "                        'FWHM': fwhm,\n",
    "                        'Wavelength': wavelengths\n",
    "                    }\n",
    "                },\n",
    "                'Reflectance_Data': np.transpose(reflectance_data, axes=(1, 2, 0))  # Convert to (rows, cols, bands)\n",
    "                # 'Reflectance_Data': np.flipud(np.transpose(reflectance_data, axes=(1, 2, 0)))  # Flip along Y-axis if needed\n",
    "\n",
    "            }\n",
    "        }\n",
    "    }\n",
    "\n",
    "    # Save to HDF5\n",
    "    h5_filename = os.path.join(data_dir,'NEON_D13_NIWO_' + os.path.basename(reflectance_tiff_path).replace('.tif', '.h5'))\n",
    "    with h5py.File(h5_filename, \"w\") as hdf_file:\n",
    "        create_h5_file_from_dict(h5_data, hdf_file)\n",
    "\n",
    "    print(f\"HDF5 file created: {h5_filename}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "1915d977-1666-4b1b-a533-9176a6879b65",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Band1 has shape (835, 1179)\n",
      "eastings shape (835, 1179)\n"
     ]
    }
   ],
   "source": [
    "import ephem\n",
    "import rasterio\n",
    "import numpy as np\n",
    "from rasterio import Affine as A\n",
    "from rasterio.warp import reproject, Resampling\n",
    "import h5py\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "# from osgeo import gdal\n",
    "import seaborn as sns\n",
    "import rasterio\n",
    "from pyproj import Proj, transform\n",
    "from rasterio.transform import from_origin\n",
    "import ephem\n",
    "import datetime\n",
    "import math\n",
    "import warnings\n",
    "\n",
    "def pixel_to_coord(file_path):\n",
    "     with rasterio.open(file_path) as src:\n",
    "            band1 = src.read(1)\n",
    "            print('Band1 has shape', band1.shape)\n",
    "            height = band1.shape[0]\n",
    "            width = band1.shape[1]\n",
    "            cols, rows = np.meshgrid(np.arange(width), np.arange(height))\n",
    "            xs, ys = rasterio.transform.xy(src.transform, rows, cols)\n",
    "            eastings= np.array(xs)\n",
    "            northings = np.array(ys)\n",
    "            print('eastings shape', eastings.shape)\n",
    "            p1 = Proj(src.crs)\n",
    "            p2 = Proj(proj='latlong', datum='WGS84')\n",
    "            lons, lats = transform(p1, p2, eastings, northings)\n",
    "     return lons, lats, cols, rows,\n",
    "\n",
    "\n",
    "file_path = file_path #r'/home/jovyan/data-store/cross-sensor-cal/data_set2/aligned_orthomosaic.tif'\n",
    "longitudes, latitudes, cols, rows = pixel_to_coord(file_path)\n",
    "\n",
    "# # # Function to convert lat/lon to row/col\n",
    "# # def latlon_to_rowcol(transform, lat, lon):\n",
    "# #     col, row = ~transform * (lon, lat)\n",
    "# #     return int(row), int(col)\n",
    "\n",
    "# Convert the date and time to UTC. The time given is 2:34 PM, which is 14:34 in 24-hour format\n",
    "date_time_str = date_time_str #'2023-08-01 21:34:00'\n",
    "\n",
    "\n",
    "with rasterio.open(file_path) as src:\n",
    "    # Get the affine transform for the raster\n",
    "    transform = src.transform\n",
    "    \n",
    "    # Create arrays to hold the azimuth and zenith values\n",
    "    azimuth = np.zeros((src.height, src.width), dtype=np.float32)\n",
    "    zenith = np.zeros((src.height, src.width), dtype=np.float32)\n",
    "    \n",
    "    # Assume a date and time for the Sun position calculation\n",
    "    observer = ephem.Observer()\n",
    "    observer.date = ephem.date(date_time_str)\n",
    "    \n",
    "    # Iterate over each pixel in the raster\n",
    "    for row in range(latitudes.shape[0]):\n",
    "        for col in range(latitudes.shape[1]):\n",
    "            #lon, lat = pixel_to_coord(transform, row, col)\n",
    "            observer.lat, observer.lon = latitudes[row,col], longitudes[row,col]\n",
    "            \n",
    "            sun = ephem.Sun(observer)\n",
    "            \n",
    "            # Convert azimuth and altitude (zenith angle is 90 - altitude) to degrees\n",
    "            az = np.degrees(sun.az)\n",
    "            #az = sun.az\n",
    "            alt = np.degrees(sun.alt)\n",
    "            zen = 90 - sun.alt\n",
    "            \n",
    "            azimuth[row, col] = az\n",
    "            zenith[row, col] = zen\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "7d7777cf-1a43-4a94-b484-43946a947c6c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "File exists: True\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "\n",
    "file_path = \"/home/jovyan/data-store/cross-sensor-cal/data_set2/aligned_orthomosaic.tif\"\n",
    "print(\"File exists:\", os.path.exists(file_path))"
   ]
  },
  {
   "cell_type": "raw",
   "id": "0425bf19-dbea-4148-b511-0928c6b6febf",
   "metadata": {},
   "source": [
    "import ephem\n",
    "import rasterio\n",
    "import numpy as np\n",
    "from rasterio.warp import reproject, Resampling\n",
    "from pyproj import Transformer\n",
    "import datetime\n",
    "\n",
    "def pixel_to_coord(file_path):\n",
    "    with rasterio.open(file_path) as src:\n",
    "        band1 = src.read(1)\n",
    "        print('Band1 has shape:', band1.shape)\n",
    "\n",
    "        height, width = band1.shape\n",
    "        cols, rows = np.meshgrid(np.arange(width), np.arange(height))\n",
    "        xs, ys = rasterio.transform.xy(src.transform, rows, cols)\n",
    "\n",
    "        eastings = np.array(xs)\n",
    "        northings = np.array(ys)\n",
    "        \n",
    "        print('Eastings shape:', eastings.shape)\n",
    "\n",
    "        # Define the transformation\n",
    "        transformer = Transformer.from_crs(src.crs, \"EPSG:4326\", always_xy=True)\n",
    "        lons, lats = transformer.transform(eastings, northings)\n",
    "\n",
    "        # Flip latitude if needed (common issue with GeoTIFFs)\n",
    "        if lats[0, 0] < lats[-1, 0]:  # Check if latitudes are inverted\n",
    "            lats = np.flipud(lats)\n",
    "\n",
    "    return np.array(lons), np.array(lats), cols, rows\n",
    "\n",
    "\n",
    "file_path = file_path  # Provide the correct file path\n",
    "date_time_str = date_time_str  # Correct date format\n",
    "\n",
    "longitudes, latitudes, cols, rows = pixel_to_coord(file_path)\n",
    "\n",
    "# Open the raster again to get the transform\n",
    "with rasterio.open(file_path) as src:\n",
    "    transform = src.transform\n",
    "\n",
    "    # Create arrays to hold the azimuth and zenith values\n",
    "    azimuth = np.zeros((src.height, src.width), dtype=np.float32)\n",
    "    zenith = np.zeros((src.height, src.width), dtype=np.float32)\n",
    "\n",
    "    # Set up the observer for ephem\n",
    "    observer = ephem.Observer()\n",
    "    observer.date = ephem.date(date_time_str)\n",
    "\n",
    "    # Iterate over each pixel in the raster\n",
    "    for row in range(latitudes.shape[0]):\n",
    "        for col in range(latitudes.shape[1]):\n",
    "            observer.lat, observer.lon = str(latitudes[row, col]), str(longitudes[row, col])\n",
    "\n",
    "            sun = ephem.Sun(observer)\n",
    "\n",
    "            # Convert azimuth and altitude (zenith angle is 90 - altitude) to degrees\n",
    "            az = np.degrees(sun.az)\n",
    "            alt = np.degrees(sun.alt)\n",
    "            zen = 90 - alt\n",
    "\n",
    "            azimuth[row, col] = az\n",
    "            zenith[row, col] = zen\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "449f3ca2-6560-4641-a8c8-69d6d1c0a0a7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HDF5 file created: /home/jovyan/data-store/cross-sensor-cal/data_set2/NEON_D13_NIWO_aligned_orthomosaic.h5\n",
      "6\n"
     ]
    }
   ],
   "source": [
    "tiff_to_h5(file_path, slope_path, aspect_path, chm_path, zenith, azimuth)\n",
    "print(\"6\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "4ee07684-b48a-4dd1-bd5e-e01e59357617",
   "metadata": {},
   "outputs": [],
   "source": [
    "NEON_image_path = '/home/jovyan/data-store/cross-sensor-cal/data_set2/NEON_D13_NIWO_aligned_orthomosaic.h5'\n",
    "output_dir = '/home/jovyan/data-store/cross-sensor-cal/output'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "1eb9d0b5-28bc-4ad5-8197-8d979edbe1f4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "print(os.path.exists(NEON_image_path))  # Should return True\n",
    "print(os.path.isdir(output_dir))  # Should return True\n"
   ]
  },
  {
   "cell_type": "raw",
   "id": "4d7b433b-a780-4421-a571-58f6ffec9f06",
   "metadata": {},
   "source": [
    "row_image_paths = '/home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_testaligned_orthomosaic/NEON_D13_NIWO_testaligned_orthomosaic'\n",
    "\n",
    "el_spectral.show_rgb(file_paths)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "86e92e58-88be-473c-8dcd-6d617def52bb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using script: /home/jovyan/data-store/cross-sensor-cal/neon2envi2_Drone.py\n"
     ]
    }
   ],
   "source": [
    "import subprocess\n",
    "import os\n",
    "\n",
    "neon_script = os.path.abspath(\"neon2envi2_Drone.py\")  # Get absolute path\n",
    "print(f\"Using script: {neon_script}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "ef886ef2-ab98-4412-b852-b176ebae4ee3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import subprocess\n",
    "import os\n",
    "\n",
    "def process_hdf5_with_neon2envi(image_path, output_dir):\n",
    "    \"\"\"\n",
    "    Runs the modified neon2envi2_generic.py script.\n",
    "    \"\"\"\n",
    "    # Ensure we're using the correct Python script\n",
    "    neon_script = os.path.abspath(\"neon2envi2_generic.py\")  # Get absolute path\n",
    "    print(f\"Using script: {neon_script}\")\n",
    "\n",
    "    # Define the command\n",
    "    command = [\n",
    "        \"/opt/conda/envs/macrosystems/bin/python\", neon_script,\n",
    "        \"-anc\",\n",
    "        \"--images\", image_path,\n",
    "        \"--output_dir\", output_dir\n",
    "    ]\n",
    "\n",
    "    print(\"Executing command:\", \" \".join(command))  # Debugging line\n",
    "\n",
    "    try:\n",
    "        with subprocess.Popen(command, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, text=True) as proc:\n",
    "            for line in proc.stdout:\n",
    "                print(line, end='')  # Print output as it's received\n",
    "    except subprocess.CalledProcessError as e:\n",
    "        print(f\"Error executing command: {e}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "357dcc39-5926-46c4-99c0-9048a5f4610b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using script: /home/jovyan/data-store/cross-sensor-cal/neon2envi2_generic.py\n",
      "Executing command: /opt/conda/envs/macrosystems/bin/python /home/jovyan/data-store/cross-sensor-cal/neon2envi2_generic.py -anc --images /home/jovyan/data-store/cross-sensor-cal/data_set2/NEON_D13_NIWO_aligned_orthomosaic.h5 --output_dir /home/jovyan/data-store/cross-sensor-cal/output\n",
      "2025-01-31 02:48:27,542\tWARNING services.py:1996 -- WARNING: The object store is using /tmp instead of /dev/shm because /dev/shm has only 67108864 bytes available. This will harm performance! You may be able to free up space by deleting files in /dev/shm. If you are inside a Docker container, you can increase /dev/shm size by passing '--shm-size=10.24gb' to 'docker run' (or add it to the run_options list in a Ray cluster config). Make sure to set this to more than 30% of available RAM.\n",
      "2025-01-31 02:48:28,676\tINFO worker.py:1673 -- Started a local Ray instance.\n",
      "Here we GO!\n",
      "\n",
      "Exporting ancillary data\n",
      "\u001b[36m(HyTools pid=39860)\u001b[0m Exported NEON_D13_NIWO_aligned_orthomosaic to ENVI format.\n",
      "\u001b[36m(HyTools pid=39860)\u001b[0m Ancillary_Imagery/Path_Length /NIWO/Reflectance/Metadata/Ancillary_Imagery/Path_Length\n",
      "\u001b[36m(HyTools pid=39860)\u001b[0m to-sensor_Azimuth_Angle /NIWO/Reflectance/Metadata/to-sensor_Azimuth_Angle\n",
      "\u001b[36m(HyTools pid=39860)\u001b[0m to-sensor_Zenith_Angle /NIWO/Reflectance/Metadata/to-sensor_Zenith_Angle\n",
      "\u001b[36m(HyTools pid=39860)\u001b[0m Logs/Solar_Azimuth_Angle /NIWO/Reflectance/Metadata/Logs/Solar_Azimuth_Angle\n",
      "\u001b[36m(HyTools pid=39860)\u001b[0m Logs/Solar_Zenith_Angle /NIWO/Reflectance/Metadata/Logs/Solar_Zenith_Angle\n",
      "\u001b[36m(HyTools pid=39860)\u001b[0m Ancillary_Imagery/Slope /NIWO/Reflectance/Metadata/Ancillary_Imagery/Slope\n",
      "\u001b[36m(HyTools pid=39860)\u001b[0m Ancillary_Imagery/Aspect /NIWO/Reflectance/Metadata/Ancillary_Imagery/Aspect\n",
      "\u001b[36m(HyTools pid=39860)\u001b[0m [984465, 984465, 984465, 984465, 984465, 984465, 984465]\n",
      "Processing complete.\n",
      "\u001b[36m(HyTools pid=39860)\u001b[0m Ancillary data exported to /home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/NEON_D13_NIWO_aligned_orthomosaic_ancillary\n",
      "done\n"
     ]
    }
   ],
   "source": [
    "process_hdf5_with_neon2envi(NEON_image_path, output_dir)\n",
    "print('done')"
   ]
  },
  {
   "cell_type": "raw",
   "id": "368868c1-0dca-4156-8f20-1e39bdb904de",
   "metadata": {},
   "source": [
    "import rasterio\n",
    "import numpy as np\n",
    "from rasterio.transform import Affine\n",
    "\n",
    "def flip_envi_and_update_transform(envi_path, output_path, tiff_transform):\n",
    "    \"\"\"\n",
    "    Flip an ENVI image vertically and update its geotransformation to match the TIFF file.\n",
    "    \n",
    "    Parameters:\n",
    "        envi_path (str): Path to the ENVI file.\n",
    "        output_path (str): Path to save the corrected ENVI file.\n",
    "        tiff_transform (Affine): The correct transformation from the TIFF file.\n",
    "    \"\"\"\n",
    "    with rasterio.open(envi_path) as src_envi:\n",
    "        # Read ENVI data and metadata\n",
    "        envi_data = src_envi.read()\n",
    "        envi_meta = src_envi.meta.copy()\n",
    "\n",
    "    # Check if ENVI needs flipping (y-resolution positive)\n",
    "    if envi_meta['transform'][5] > 0:\n",
    "        print(\"Flipping ENVI image to match TIFF orientation...\")\n",
    "        # Flip the data vertically\n",
    "        flipped_data = np.flip(envi_data, axis=1)\n",
    "\n",
    "        # Update the transformation to match the TIFF file\n",
    "        new_transform = Affine(\n",
    "            tiff_transform.a, tiff_transform.b, tiff_transform.c,\n",
    "            tiff_transform.d, tiff_transform.e, tiff_transform.f\n",
    "        )\n",
    "        envi_meta.update({\"transform\": new_transform})\n",
    "    else:\n",
    "        flipped_data = envi_data  # No flipping needed\n",
    "        print(\"ENVI image does not need flipping.\")\n",
    "\n",
    "    # Save the corrected ENVI file\n",
    "    with rasterio.open(output_path, \"w\", **envi_meta) as dst:\n",
    "        dst.write(flipped_data)\n",
    "\n",
    "    print(f\"Corrected ENVI file saved to {output_path}\")\n",
    "\n",
    "# Example usage\n",
    "envi_path = raster_path\n",
    "output_path = raster_path\n",
    "tiff_transform = Affine(0.25, 0.0, 457602.75, 0.0, -0.25, 4432862.25)  # From TIFF metadata\n",
    "\n",
    "flip_envi_and_update_transform(envi_path, output_path, tiff_transform)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "8cc723b3-9ba5-48cb-b2cb-00763136bb8c",
   "metadata": {},
   "source": [
    "import rasterio\n",
    "import numpy as np\n",
    "\n",
    "def fix_envi_orientation(envi_path, output_path):\n",
    "    \"\"\"\n",
    "    Reflip the ENVI image vertically while preserving its corrected spatial position.\n",
    "    \n",
    "    Parameters:\n",
    "        envi_path (str): Path to the ENVI file.\n",
    "        output_path (str): Path to save the final corrected ENVI file.\n",
    "    \"\"\"\n",
    "    with rasterio.open(envi_path) as src:\n",
    "        # Read ENVI data and metadata\n",
    "        envi_data = src.read()\n",
    "        envi_meta = src.meta.copy()\n",
    "\n",
    "        # Reflip the ENVI image vertically\n",
    "        print(\"Reflipping ENVI image to restore correct orientation...\")\n",
    "        reflipped_data = np.flip(envi_data, axis=1)  # Flip along the y-axis (rows)\n",
    "\n",
    "    # Save the final corrected ENVI file\n",
    "    with rasterio.open(output_path, \"w\", **envi_meta) as dst:\n",
    "        dst.write(reflipped_data)\n",
    "\n",
    "    print(f\"Final corrected ENVI file saved to {output_path}\")\n",
    "\n",
    "# Example usage\n",
    "envi_path = raster_path\n",
    "output_path = raster_path\n",
    "\n",
    "fix_envi_orientation(envi_path, output_path)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "68c099a1-435d-4afb-ac56-da556240c559",
   "metadata": {},
   "source": [
    "with rasterio.open(file_path) as src:\n",
    "    print(\"TIFF Metadata:\", src.meta)\n",
    "\n",
    "with rasterio.open('/home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/NEON_D13_NIWO_aligned_orthomosaic') as src:\n",
    "    print(\"ENVI Metadata:\", src.meta)\n"
   ]
  },
  {
   "cell_type": "raw",
   "id": "34002191-50c4-469f-97cf-63aec41cc2b4",
   "metadata": {},
   "source": [
    "import rasterio\n",
    "import numpy as np\n",
    "import hytools as ht\n",
    "import re\n",
    "import os\n",
    "\n",
    "def modify_envi_to_match_tiff(envi_path, output_path, tiff_path):\n",
    "    \"\"\"\n",
    "    Modify an ENVI file to match the geotransformation of a TIFF file.\n",
    "    \n",
    "    Parameters:\n",
    "        envi_path (str): Path to the ENVI data file (without .hdr extension).\n",
    "        output_path (str): Path to save the modified ENVI file.\n",
    "        tiff_path (str): Path to the reference TIFF file.\n",
    "    \"\"\"\n",
    "    # Load the ENVI file using HyTools\n",
    "    hy_obj = ht.HyTools()\n",
    "    hy_obj.read_file(envi_path)  # Corrected function call\n",
    "    hy_obj.load_data()\n",
    "\n",
    "    # Read the TIFF file's geotransformation\n",
    "    with rasterio.open(tiff_path) as tiff_src:\n",
    "        tiff_transform = tiff_src.transform\n",
    "        tiff_crs = tiff_src.crs\n",
    "        tiff_meta = tiff_src.meta\n",
    "\n",
    "    # Extract map info from TIFF (needed for ENVI header)\n",
    "    tiff_map_info = f\"UTM, 1.000, 1.000, {tiff_transform.c}, {tiff_transform.f}, {tiff_transform.a}, {tiff_transform.e}, {tiff_crs.to_epsg()}, North, WGS-84, units=Meters\"\n",
    "\n",
    "        # Extract map info correctly\n",
    "    envi_map_info_list = hy_obj.get_header().get('map info', [])\n",
    "    \n",
    "    if envi_map_info_list:  # Ensure it's not empty\n",
    "        envi_map_info = envi_map_info_list[0]  # Extract the first string element\n",
    "        envi_map_info_parts = envi_map_info.split(',')  # Split by commas\n",
    "    \n",
    "        if len(envi_map_info_parts) > 5:  # Ensure the list has at least 6 elements\n",
    "            match = re.search(r\"([-+]?\\d*\\.\\d+|\\d+)\", envi_map_info_parts[5])  # Extract y-resolution\n",
    "            if match:\n",
    "                y_res_envi = float(match.group(0))\n",
    "                if y_res_envi > 0:  # If y-resolution is positive, flip the ENVI data\n",
    "                    print(\"Flipping ENVI data to match TIFF orientation...\")\n",
    "                    hy_obj.data = np.flip(hy_obj.data, axis=0)  # Flip vertically\n",
    "        else:\n",
    "            print(\"Warning: 'map info' does not contain enough values to extract y-resolution.\")\n",
    "\n",
    "\n",
    "    # Save modified ENVI file\n",
    "    output_hdr_path = output_path + \".hdr\"\n",
    "    output_data_path = output_path\n",
    "\n",
    "    print(f\"Saving modified ENVI file to {output_data_path}\")\n",
    "\n",
    "    # Write modified ENVI file\n",
    "    hy_obj.write_envi(output_data_path)\n",
    "\n",
    "    # Modify the ENVI header manually to match the TIFF geotransformation\n",
    "    envi_hdr_path = envi_path + \".hdr\"\n",
    "    if os.path.exists(envi_hdr_path):\n",
    "        with open(envi_hdr_path, \"r\") as hdr_file:\n",
    "            header_lines = hdr_file.readlines()\n",
    "\n",
    "        with open(output_hdr_path, \"w\") as hdr_file:\n",
    "            for line in header_lines:\n",
    "                if \"map info\" in line:\n",
    "                    hdr_file.write(f\"map info = {tiff_map_info}\\n\")\n",
    "                else:\n",
    "                    hdr_file.write(line)\n",
    "\n",
    "    print(f\"Modified ENVI file saved at: {output_data_path}\")\n",
    "    print(f\"Updated ENVI header saved at: {output_hdr_path}\")\n",
    "\n",
    "# Example usage\n",
    "envi_path = \"/home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/NEON_D13_NIWO_aligned_orthomosaic\"  # ENVI file (without .hdr extension)\n",
    "output_path = \"/home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/NEON_D13_NIWO_aligned_orthomosaic_modified\"  # Path to save the modified ENVI file\n",
    "tiff_path = file_path  # Reference TIFF file\n",
    "\n",
    "# Run the function\n",
    "modify_envi_to_match_tiff(envi_path, output_path, tiff_path)\n"
   ]
  },
  {
   "cell_type": "raw",
   "id": "434f24e8-a6ac-44c7-b16e-234da39dd1a4",
   "metadata": {},
   "source": [
    "# Input and output paths\n",
    "envi_path = \"/home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/NEON_D13_NIWO_aligned_orthomosaic\"  # Path to the original ENVI data file\n",
    "envi_hdr_path = \"/home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/NEON_D13_NIWO_aligned_orthomosaic.hdr\"  # Path to the ENVI header file\n",
    "output_path = \"/home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/NEON_D13_NIWO_aligned_orthomosaic_modified\"  # Path to save the modified ENVI file\n",
    "\n",
    "\n",
    "# Run the function\n",
    "modify_envi_to_match_tiff(envi_path, envi_hdr_path, output_path, file_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "e758772a-a53b-4348-b253-5385952beef5",
   "metadata": {},
   "outputs": [],
   "source": [
    "directory = '/home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "b1cd251d-f6b0-4cf1-a2c3-c000edb71d22",
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "import glob\n",
    "import numpy as np\n",
    "import os\n",
    "\n",
    "# Define your configuration settings\n",
    "bad_bands = []\n",
    "file_type = 'envi'\n",
    "directory = directory\n",
    "main_image_name = os.path.basename(directory)\n",
    "main_image_file = os.path.join(directory, main_image_name )  # Assuming the main image file has .h5 extension\n",
    "\n",
    "# Glob pattern to find ancillary files within the same directory\n",
    "anc_files_pattern = os.path.join(directory, \"*_ancillary*\")\n",
    "anc_files = glob.glob(anc_files_pattern)\n",
    "anc_files.sort()\n",
    "\n",
    "aviris_anc_names = ['path_length', 'sensor_az', 'sensor_zn', 'solar_az', 'solar_zn',  'slope', 'aspect', 'phase', 'cosine_i']\n",
    "\n",
    "suffix_labels = [\"reflectance\", \"anc\"]  # Define suffixes for different types of files\n",
    "\n",
    "# Loop through each ancillary file and create a separate config file\n",
    "for i, anc_file in enumerate(anc_files):\n",
    "    if i == 0:\n",
    "        suffix_label = f\"_reflectance\"\n",
    "    elif i == 1:\n",
    "        suffix_label = f\"_anc\"\n",
    "    else:\n",
    "        suffix_label = f\"_{i}\"  # Fallback for unexpected 'i' values\n",
    "\n",
    "    \n",
    "    \n",
    "    config_dict = {}\n",
    "\n",
    "    config_dict['bad_bands'] = bad_bands\n",
    "    config_dict['file_type'] = file_type\n",
    "    config_dict[\"input_files\"] = [main_image_file]\n",
    "    \n",
    "    config_dict[\"anc_files\"] = {\n",
    "        main_image_file: dict(zip(aviris_anc_names, [[anc_file, a] for a in range(len(aviris_anc_names))]))\n",
    "    }\n",
    "\n",
    "    # Export settings\n",
    "    config_dict['export'] = {}\n",
    "    config_dict['export']['coeffs'] = True\n",
    "    config_dict['export']['image'] = True\n",
    "    config_dict['export']['masks'] = True\n",
    "    config_dict['export']['subset_waves'] = []\n",
    "    config_dict['export']['output_dir'] = os.path.join(directory)\n",
    "    config_dict['export'][\"suffix\"] = suffix_label\n",
    "\n",
    "\n",
    "\n",
    "    # Detailed settings for export options, TOPO and BRDF corrections\n",
    "    # These settings include parameters like the type of correction, calculation and application of masks, \n",
    "    # coefficients to be used, output directory, suffixes for output files, and various specific parameters \n",
    "    # for TOPO and BRDF correction methods.\n",
    "    # Input format: Nested dictionaries with specific keys and values as per the correction algorithm requirements\n",
    "    # Example settings include:\n",
    "    # - 'export': Dictionary of export settings like coefficients, image, masks, output directory, etc.\n",
    "    # - 'topo': Dictionary of topographic correction settings including types, masks, coefficients, etc.\n",
    "    # - 'brdf': Dictionary of BRDF correction settings including solar zenith type, geometric model, volume model, etc.\n",
    "\n",
    "    # Additional settings can be added as needed for specific correction algorithms and export requirements.\n",
    "\n",
    "    # TOPO Correction options\n",
    "    config_dict[\"corrections\"] = ['topo','brdf']\n",
    "    config_dict[\"topo\"] =  {}\n",
    "    config_dict[\"topo\"]['type'] = 'scs+c'\n",
    "    config_dict[\"topo\"]['calc_mask'] = [[\"ndi\", {'band_1': 862,'band_2': 668,\n",
    "                                                'min': 0.1,'max': 1.0}],\n",
    "                                        ['ancillary',{'name':'slope',\n",
    "                                                    'min': np.radians(5),'max':'+inf' }],\n",
    "                                        ['ancillary',{'name':'cosine_i',\n",
    "                                                    'min': 0.12,'max':'+inf' }]]\n",
    "\n",
    "    config_dict[\"topo\"]['apply_mask'] = [[\"ndi\", {'band_1': 862,'band_2': 668,\n",
    "                                                'min': 0.1,'max': 1.0}],\n",
    "                                        ['ancillary',{'name':'slope',\n",
    "                                                    'min': np.radians(5),'max':'+inf' }],\n",
    "                                        ['ancillary',{'name':'cosine_i',\n",
    "                                                    'min': 0.12,'max':'+inf' }]]\n",
    "    config_dict[\"topo\"]['c_fit_type'] = 'nnls'\n",
    "\n",
    "    # config_dict[\"topo\"]['type'] =  'precomputed'\n",
    "    # config_dict[\"brdf\"]['coeff_files'] =  {}\n",
    "\n",
    "    # BRDF Correction options\n",
    "    config_dict[\"brdf\"] = {}\n",
    "    config_dict[\"brdf\"]['solar_zn_type'] ='scene'\n",
    "    config_dict[\"brdf\"]['type'] = 'flex'\n",
    "    config_dict[\"brdf\"]['grouped'] = True\n",
    "    config_dict[\"brdf\"]['sample_perc'] = 0.1\n",
    "    config_dict[\"brdf\"]['geometric'] = 'li_dense_r'\n",
    "    config_dict[\"brdf\"]['volume'] = 'ross_thick'\n",
    "    config_dict[\"brdf\"][\"b/r\"] = 10  #these may need updating. These constants pulled from literature. \n",
    "    config_dict[\"brdf\"][\"h/b\"] = 2  # These may need updating. These contanstants pulled from literature.\n",
    "    config_dict[\"brdf\"]['interp_kind'] = 'linear'\n",
    "    config_dict[\"brdf\"]['calc_mask'] = [[\"ndi\", {'band_1': 862, 'band_2': 668, 'min': 0.1, 'max': 1.0}]]\n",
    "    config_dict[\"brdf\"]['apply_mask'] = [[\"ndi\", {'band_1': 862, 'band_2': 668, 'min': 0.1, 'max': 1.0}]]\n",
    "    config_dict[\"brdf\"]['diagnostic_plots'] = True\n",
    "    config_dict[\"brdf\"]['diagnostic_waves'] = [475, 862, 705, 668]\n",
    "\n",
    "    # ## Flex dynamic NDVI params\n",
    "    config_dict[\"brdf\"]['bin_type'] = 'dynamic'\n",
    "    config_dict[\"brdf\"]['num_bins'] = 25\n",
    "    config_dict[\"brdf\"]['ndvi_bin_min'] = 0.05\n",
    "    config_dict[\"brdf\"]['ndvi_bin_max'] = 1.0\n",
    "    config_dict[\"brdf\"]['ndvi_perc_min'] = 10\n",
    "    config_dict[\"brdf\"]['ndvi_perc_max'] = 95\n",
    "\n",
    "    # Define the number of CPUs to be used (considering the number of image-ancillary pairs)\n",
    "    config_dict['num_cpus'] = 1\n",
    "    \n",
    "    # Output path for configuration file\n",
    "    # Assuming you want to include the suffix in the filename:\n",
    "    suffix = config_dict['export'][\"suffix\"]\n",
    "    output_dir = config_dict['export']['output_dir']\n",
    "    config_file_name = f\"{suffix}.json\"\n",
    "    config_file_path = os.path.join(output_dir, config_file_name)\n",
    "\n",
    "    config_dict[\"resample\"]  = False\n",
    "    config_dict[\"resampler\"]  = {}\n",
    "    config_dict[\"resampler\"]['type'] =  'cubic'\n",
    "    config_dict[\"resampler\"]['out_waves'] = []\n",
    "    config_dict[\"resampler\"]['out_fwhm'] = []\n",
    "\n",
    "    # # Remove bad bands from output waves\n",
    "    # for wavelength in range(450,660,100):\n",
    "    #     bad=False\n",
    "    #     for start,end in config_dict['bad_bands']:\n",
    "    #         bad = ((wavelength >= start) & (wavelength <=end)) or bad\n",
    "    #     if not bad:\n",
    "    #         config_dict[\"resampler\"]['out_waves'].append(wavelength)\n",
    "\n",
    "\n",
    "    # Construct the filename for the configuration JSON\n",
    "    config_filename = f\"{main_image_name}_config_{suffix}.json\"\n",
    "    config_file_path = os.path.join(directory, config_filename)\n",
    "\n",
    "    # Save the configuration to a JSON file\n",
    "    with open(config_file_path, 'w') as outfile:\n",
    "        json.dump(config_dict, outfile, indent=3)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "515fde35-3086-4214-b149-aebb595b4f78",
   "metadata": {},
   "outputs": [],
   "source": [
    "import subprocess\n",
    "import os\n",
    "\n",
    "def brdf_topo_correct_image(jason_path):\n",
    "    \"\"\"\n",
    "    Runs the modified neon2envi2_generic.py script.\n",
    "    \"\"\"\n",
    "    # Ensure we're using the correct Python script\n",
    "    neon_script = os.path.abspath(\"image_correct.py\")  # Get absolute path\n",
    "    print(f\"Using script: {neon_script}\")\n",
    "\n",
    "    # Define the command\n",
    "    command = [\n",
    "        \"/opt/conda/envs/macrosystems/bin/python\", neon_script,\n",
    "        jason_path\n",
    "    ]\n",
    "\n",
    "    print(\"Executing command:\", \" \".join(command))  # Debugging line\n",
    "\n",
    "    try:\n",
    "        with subprocess.Popen(command, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, text=True) as proc:\n",
    "            for line in proc.stdout:\n",
    "                print(line, end='')  # Print output as it's received\n",
    "    except subprocess.CalledProcessError as e:\n",
    "        print(f\"Error executing command: {e}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "94d371e9-ffe3-49ce-bcb6-ea698622c436",
   "metadata": {},
   "outputs": [],
   "source": [
    "jason_path = '/home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/NEON_D13_NIWO_aligned_orthomosaic_config__reflectance.json'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "23ea3369-3c11-4955-b468-5fd6ef0a2d12",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using script: /home/jovyan/data-store/cross-sensor-cal/image_correct.py\n",
      "Executing command: /opt/conda/envs/macrosystems/bin/python /home/jovyan/data-store/cross-sensor-cal/image_correct.py /home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/NEON_D13_NIWO_aligned_orthomosaic_config__reflectance.json\n",
      "2025-01-31 02:48:35,284\tWARNING services.py:1996 -- WARNING: The object store is using /tmp instead of /dev/shm because /dev/shm has only 67108864 bytes available. This will harm performance! You may be able to free up space by deleting files in /dev/shm. If you are inside a Docker container, you can increase /dev/shm size by passing '--shm-size=10.24gb' to 'docker run' (or add it to the run_options list in a Ray cluster config). Make sure to set this to more than 30% of available RAM.\n",
      "2025-01-31 02:48:36,421\tINFO worker.py:1673 -- Started a local Ray instance.\n",
      "/home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/NEON_D13_NIWO_aligned_orthomosaic_config__reflectance.json\n",
      "Using 1 CPUs.\n",
      "Calculating topographic coefficients.\n",
      "Scene average solar zenith angle : 89.023 degrees\n",
      "Calculating BRDF coefficients\n",
      "\n",
      " ##########------------------------------------------------------------------------------------------ [ 10.00%]\n",
      " ####################-------------------------------------------------------------------------------- [ 20.00%]\n",
      " ##############################---------------------------------------------------------------------- [ 30.00%]\n",
      " ########################################------------------------------------------------------------ [ 40.00%]\n",
      " ##################################################-------------------------------------------------- [ 50.00%]\n",
      " ############################################################---------------------------------------- [ 60.00%]\n",
      " ######################################################################------------------------------ [ 70.00%]\n",
      " ################################################################################-------------------- [ 80.00%]\n",
      " ##########################################################################################---------- [ 90.00%]\n",
      " #################################################################################################### [100.00%]\n",
      "\n",
      "Exporting correction coefficients.\n",
      "Exporting corrected images.\n",
      "done\n"
     ]
    }
   ],
   "source": [
    "brdf_topo_correct_image(jason_path)\n",
    "print('done')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "3d31749b-4f3a-4ef2-b630-b3451d1ffa91",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/NEON_D13_NIWO_aligned_orthomosaic__reflectance_corrected'"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "# chnage here based on AOP underflights or not. If underflights, the polygons are at the AOP geojson file. If not within the AOP, the polygons are at uas geojson file\n",
    "# geojson_path = '/home/jovyan/data-store/cross-sensor-cal/uas_polygons_1_24_2025_analysis_ready.geojson'\n",
    "# geojson_path_half = '/home/jovyan/data-store/cross-sensor-cal/uas_polygons_1_24_2025_analysis_ready_half_diam.geojson'\n",
    "\n",
    "geojson_path = '/home/jovyan/data-store/cross-sensor-cal/aop_macrosystems_data_1_7_25.geojson'\n",
    "geojson_path_half = '/home/jovyan/data-store/cross-sensor-cal/aop_polygons_1_24_2025_analysis_ready_half_diam.geojson'\n",
    "\n",
    "\n",
    "main_image_name = os.path.basename(directory)\n",
    "raster_path = os.path.join(directory, main_image_name + '__reflectance_corrected' )\n",
    "\n",
    "output_csv = \"/home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/NEON_D13_NIWO_reflectance_per_pixel.csv\"\n",
    "output_csv_half = \"/home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/NEON_D13_NIWO_reflectance_per_pixel_half.csv\"\n",
    "\n",
    "raster_path_2 = os.path.join(directory, main_image_name + '__reflectance_corrected_flipped' )\n",
    "\n",
    "\n",
    "raster_path"
   ]
  },
  {
   "cell_type": "raw",
   "id": "5d8f6118-103b-4f05-96a1-0c705f748a98",
   "metadata": {},
   "source": [
    "import rasterio\n",
    "import numpy as np\n",
    "from rasterio.transform import Affine\n",
    "\n",
    "def flip_envi_and_update_transform(envi_path, output_path, tiff_transform):\n",
    "    \"\"\"\n",
    "    Flip an ENVI image vertically and update its geotransformation to match the TIFF file.\n",
    "    \n",
    "    Parameters:\n",
    "        envi_path (str): Path to the ENVI file.\n",
    "        output_path (str): Path to save the corrected ENVI file.\n",
    "        tiff_transform (Affine): The correct transformation from the TIFF file.\n",
    "    \"\"\"\n",
    "    with rasterio.open(envi_path) as src_envi:\n",
    "        # Read ENVI data and metadata\n",
    "        envi_data = src_envi.read()\n",
    "        envi_meta = src_envi.meta.copy()\n",
    "\n",
    "    # Check if ENVI needs flipping (y-resolution positive)\n",
    "    if envi_meta['transform'][5] > 0:\n",
    "        print(\"Flipping ENVI image to match TIFF orientation...\")\n",
    "        # Flip the data vertically\n",
    "        flipped_data = np.flip(envi_data, axis=1)\n",
    "\n",
    "        # Update the transformation to match the TIFF file\n",
    "        new_transform = Affine(\n",
    "            tiff_transform.a, tiff_transform.b, tiff_transform.c,\n",
    "            tiff_transform.d, tiff_transform.e, tiff_transform.f\n",
    "        )\n",
    "        envi_meta.update({\"transform\": new_transform})\n",
    "    else:\n",
    "        flipped_data = envi_data  # No flipping needed\n",
    "        print(\"ENVI image does not need flipping.\")\n",
    "\n",
    "    # Save the corrected ENVI file\n",
    "    with rasterio.open(output_path, \"w\", **envi_meta) as dst:\n",
    "        dst.write(flipped_data)\n",
    "\n",
    "    print(f\"Corrected ENVI file saved to {output_path}\")\n",
    "\n",
    "# Example usage\n",
    "envi_path = raster_path\n",
    "output_path = raster_path_2\n",
    "tiff_transform = Affine(0.25, 0.0, 457602.75, 0.0, -0.25, 4432862.25)  # From TIFF metadata\n",
    "\n",
    "flip_envi_and_update_transform(envi_path, output_path, tiff_transform)\n"
   ]
  },
  {
   "cell_type": "raw",
   "id": "a5363cea-8533-4dc2-af6f-f49625a4c178",
   "metadata": {},
   "source": [
    "import rasterio\n",
    "import numpy as np\n",
    "\n",
    "def fix_envi_orientation(envi_path, output_path):\n",
    "    \"\"\"\n",
    "    Reflip the ENVI image vertically while preserving its corrected spatial position.\n",
    "    \n",
    "    Parameters:\n",
    "        envi_path (str): Path to the ENVI file.\n",
    "        output_path (str): Path to save the final corrected ENVI file.\n",
    "    \"\"\"\n",
    "    with rasterio.open(envi_path) as src:\n",
    "        # Read ENVI data and metadata\n",
    "        envi_data = src.read()\n",
    "        envi_meta = src.meta.copy()\n",
    "\n",
    "        # Reflip the ENVI image vertically\n",
    "        print(\"Reflipping ENVI image to restore correct orientation...\")\n",
    "        reflipped_data = np.flip(envi_data, axis=1)  # Flip along the y-axis (rows)\n",
    "\n",
    "    # Save the final corrected ENVI file\n",
    "    with rasterio.open(output_path, \"w\", **envi_meta) as dst:\n",
    "        dst.write(reflipped_data)\n",
    "\n",
    "    print(f\"Final corrected ENVI file saved to {output_path}\")\n",
    "\n",
    "# Example usage\n",
    "envi_path = raster_path_2\n",
    "output_path = raster_path_2\n",
    "\n",
    "fix_envi_orientation(envi_path, output_path)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "ed73a901-8e52-4759-a7a7-af177ed825b6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TIFF Transform: | 0.25, 0.00, 457268.25|\n",
      "| 0.00,-0.25, 4429034.75|\n",
      "| 0.00, 0.00, 1.00|\n"
     ]
    }
   ],
   "source": [
    "import rasterio\n",
    "\n",
    "def get_tiff_transform(tiff_path):\n",
    "    \"\"\"\n",
    "    Extract the transform (Affine) from a TIFF file.\n",
    "\n",
    "    Parameters:\n",
    "        tiff_path (str): Path to the TIFF file.\n",
    "\n",
    "    Returns:\n",
    "        Affine: The geotransformation of the TIFF file.\n",
    "    \"\"\"\n",
    "    with rasterio.open(tiff_path) as src:\n",
    "        return src.transform\n",
    "\n",
    "# Example usage\n",
    "tiff_path = file_path  # Path to your reference TIFF file\n",
    "tiff_transform = get_tiff_transform(tiff_path)\n",
    "\n",
    "print(\"TIFF Transform:\", tiff_transform)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "2c417ec3-b6a1-464e-8f01-d13dea19246f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio\n",
    "import numpy as np\n",
    "from rasterio.transform import Affine\n",
    "\n",
    "def flip_envi_and_preserve_metadata(envi_path, output_path, tiff_transform):\n",
    "    \"\"\"\n",
    "    Flip an ENVI image vertically, update its geotransformation, and preserve all metadata.\n",
    "    \n",
    "    Parameters:\n",
    "        envi_path (str): Path to the ENVI file.\n",
    "        output_path (str): Path to save the corrected ENVI file.\n",
    "        tiff_transform (Affine): The correct transformation from the TIFF file.\n",
    "    \"\"\"\n",
    "    with rasterio.open(envi_path) as src_envi:\n",
    "        # Read ENVI data and metadata\n",
    "        envi_data = src_envi.read()\n",
    "        envi_meta = src_envi.meta.copy()\n",
    "        \n",
    "        # Extract additional metadata like wavelengths\n",
    "        wavelengths = src_envi.descriptions  # Band descriptions, e.g., wavelengths\n",
    "        wavelengths_units = src_envi.tags().get(\"wavelength units\", \"nanometers\")  # Default to nanometers\n",
    "\n",
    "    # Check if ENVI needs flipping (y-resolution positive)\n",
    "    if envi_meta['transform'][5] > 0:\n",
    "        print(\"Flipping ENVI image to match TIFF orientation...\")\n",
    "        # Flip the data vertically\n",
    "        flipped_data = np.flip(envi_data, axis=1)\n",
    "\n",
    "        # Update the transformation to match the TIFF file\n",
    "        new_transform = Affine(\n",
    "            tiff_transform.a, tiff_transform.b, tiff_transform.c,\n",
    "            tiff_transform.d, tiff_transform.e, tiff_transform.f\n",
    "        )\n",
    "        envi_meta.update({\"transform\": new_transform})\n",
    "    else:\n",
    "        flipped_data = envi_data  # No flipping needed\n",
    "        print(\"ENVI image does not need flipping.\")\n",
    "\n",
    "    # Save the corrected ENVI file\n",
    "    with rasterio.open(output_path, \"w\", **envi_meta) as dst:\n",
    "        dst.write(flipped_data)\n",
    "\n",
    "        # Add wavelengths back to metadata\n",
    "        if wavelengths:\n",
    "            dst.descriptions = wavelengths\n",
    "\n",
    "        # Add wavelength units if available\n",
    "        if wavelengths_units:\n",
    "            dst.update_tags(**{\"wavelength units\": wavelengths_units})\n",
    "\n",
    "    print(f\"Corrected ENVI file saved to {output_path} with preserved metadata.\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "99f700d1-e68b-45f4-b1db-9b25997e6402",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Flipping ENVI image to match TIFF orientation...\n",
      "Corrected ENVI file saved to /home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/NEON_D13_NIWO_aligned_orthomosaic__reflectance_corrected_flipped with preserved metadata.\n"
     ]
    }
   ],
   "source": [
    "envi_path = raster_path\n",
    "output_path = raster_path_2\n",
    "\n",
    "flip_envi_and_preserve_metadata(envi_path, output_path, tiff_transform)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "4371b0b8-a1fa-434c-9402-0d6dd1466adf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reflipping ENVI image to restore correct orientation...\n",
      "Corrected ENVI file saved to /home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/NEON_D13_NIWO_aligned_orthomosaic__reflectance_corrected_flipped with preserved metadata.\n"
     ]
    }
   ],
   "source": [
    "import rasterio\n",
    "import numpy as np\n",
    "\n",
    "def fix_envi_orientation(envi_path, output_path):\n",
    "    \"\"\"\n",
    "    Reflip the ENVI image vertically while preserving its corrected spatial position.\n",
    "    \n",
    "    Parameters:\n",
    "        envi_path (str): Path to the ENVI file.\n",
    "        output_path (str): Path to save the final corrected ENVI file.\n",
    "    \"\"\"\n",
    "    with rasterio.open(envi_path) as src_envi:\n",
    "        # Read ENVI data and metadata\n",
    "        envi_data = src_envi.read()\n",
    "        envi_meta = src_envi.meta.copy()\n",
    "        \n",
    "        # Extract additional metadata like wavelengths\n",
    "        wavelengths = src_envi.descriptions  # Band descriptions, e.g., wavelengths\n",
    "        wavelengths_units = src_envi.tags().get(\"wavelength units\", \"nanometers\")  # Default to nanometers\n",
    "\n",
    "        \n",
    "        # envi_data = src.read()\n",
    "        # envi_meta = src.meta.copy()\n",
    "\n",
    "        # Reflip the ENVI image vertically\n",
    "        print(\"Reflipping ENVI image to restore correct orientation...\")\n",
    "        reflipped_data = np.flip(envi_data, axis=1)  # Flip along the y-axis (rows)\n",
    "\n",
    "    # Save the corrected ENVI file\n",
    "    with rasterio.open(output_path, \"w\", **envi_meta) as dst:\n",
    "        dst.write(reflipped_data)\n",
    "\n",
    "        # Add wavelengths back to metadata\n",
    "        if wavelengths:\n",
    "            dst.descriptions = wavelengths\n",
    "\n",
    "        # Add wavelength units if available\n",
    "        if wavelengths_units:\n",
    "            dst.update_tags(**{\"wavelength units\": wavelengths_units})\n",
    "\n",
    "    print(f\"Corrected ENVI file saved to {output_path} with preserved metadata.\")\n",
    "\n",
    "# Example usage\n",
    "envi_path = raster_path_2\n",
    "output_path = raster_path_2\n",
    "\n",
    "fix_envi_orientation(envi_path, output_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "65b4f44a-cfdf-455d-8837-972c648d1ff5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "import rasterio\n",
    "import rasterio.mask\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from shapely.geometry import box\n",
    "from shapely.geometry import mapping\n",
    "\n",
    "def check_and_reproject(geojson_path, raster_path):\n",
    "    \"\"\"\n",
    "    Checks if the GeoJSON polygons and the ENVI raster have the same CRS.\n",
    "    If not, reprojects the polygons to match the raster's CRS.\n",
    "    \"\"\"\n",
    "    # Load GeoJSON\n",
    "    polygons = gpd.read_file(geojson_path)\n",
    "\n",
    "    # Load raster CRS\n",
    "    with rasterio.open(raster_path) as src:\n",
    "        raster_crs = src.crs\n",
    "\n",
    "    # Check if CRS matches\n",
    "    if polygons.crs != raster_crs:\n",
    "        print(f\"🔄 Reprojecting GeoJSON from {polygons.crs} to {raster_crs}\")\n",
    "        polygons = polygons.to_crs(raster_crs)\n",
    "    else:\n",
    "        print(\"✅ GeoJSON and Raster have the same CRS\")\n",
    "\n",
    "    return polygons\n",
    "\n",
    "def extract_pixel_reflectance(geojson_path, raster_path, output_csv):\n",
    "    \"\"\"\n",
    "    Extracts reflectance values for each pixel inside each polygon and saves to CSV.\n",
    "    \"\"\"\n",
    "    # Reproject GeoJSON if needed\n",
    "    polygons = check_and_reproject(geojson_path, raster_path)\n",
    "\n",
    "    # Open raster and get metadata\n",
    "    with rasterio.open(raster_path) as src:\n",
    "        wavelengths = src.descriptions  # Get actual wavelengths (e.g., [\"444nm\", \"475nm\", ...])\n",
    "        # # Get band descriptions or assign default names if descriptions are missing\n",
    "        # if src.descriptions and all(src.descriptions):\n",
    "        #     wavelengths = src.descriptions  # Use actual descriptions if available\n",
    "        # else:\n",
    "        #     wavelengths = [f\"Band_{i+1}\" for i in range(src.count)]  # Assign default names\n",
    "        \n",
    "        raster_crs = src.crs\n",
    "        raster_bounds = src.bounds\n",
    "        raster_transform = src.transform  # Affine transform for pixel-to-geo mapping\n",
    "\n",
    "    # Convert wavelengths to valid column names\n",
    "    wavelengths = [w.replace(\" \", \"_\") for w in wavelengths]  # Remove spaces if any\n",
    "\n",
    "    # Convert raster bounds to a polygon\n",
    "    raster_extent = box(*raster_bounds)\n",
    "\n",
    "    # Filter polygons that are within the raster extent\n",
    "    polygons = polygons[polygons.geometry.intersects(raster_extent)]\n",
    "\n",
    "    # If no polygons remain, exit\n",
    "    if polygons.empty:\n",
    "        print(\"❌ No polygons found within the raster extent!\")\n",
    "        return\n",
    "\n",
    "    print(f\"✅ {len(polygons)} polygons found within raster extent.\")\n",
    "\n",
    "    # Prepare a list to store extracted data\n",
    "    extracted_data = []\n",
    "\n",
    "    # Loop through each polygon\n",
    "    for poly_idx, polygon in polygons.iterrows():\n",
    "        polygon_id = polygon.get(\"OBJECTID\", poly_idx)  # Use an ID field if available\n",
    "\n",
    "        # Mask the raster to extract pixel values within the polygon\n",
    "        with rasterio.open(raster_path) as src:\n",
    "            out_image, out_transform = rasterio.mask.mask(src, [mapping(polygon.geometry)], crop=True)\n",
    "            out_image = out_image.astype(np.float32)  # Convert to float for precision\n",
    "\n",
    "            # Get pixel coordinates\n",
    "            num_rows, num_cols = out_image.shape[1], out_image.shape[2]\n",
    "            x_coords = np.arange(num_cols) * out_transform[0] + out_transform[2]\n",
    "            y_coords = np.arange(num_rows) * out_transform[4] + out_transform[5]\n",
    "\n",
    "            # Loop over each pixel\n",
    "            pixel_id = 0\n",
    "            for i in range(num_rows):\n",
    "                for j in range(num_cols):\n",
    "                    # Extract reflectance for all bands\n",
    "                    reflectance_values = out_image[:, i, j]\n",
    "\n",
    "                    # Check if pixel is valid (not nodata)\n",
    "                    if np.any(reflectance_values > 0):  # Assuming negative values are NoData\n",
    "                        row = {\n",
    "                            \"Polygon_ID\": polygon_id,\n",
    "                            \"Pixel_ID\": f\"{polygon_id}_{pixel_id}\",  # Unique ID for each pixel\n",
    "                            \"X_Coordinate\": x_coords[j],\n",
    "                            \"Y_Coordinate\": y_coords[i],\n",
    "                            **{wavelengths[b]: reflectance_values[b] for b in range(len(wavelengths))}\n",
    "                        }\n",
    "                        extracted_data.append(row)\n",
    "                        pixel_id += 1  # Increment pixel counter\n",
    "\n",
    "    # Convert to DataFrame\n",
    "    df = pd.DataFrame(extracted_data)\n",
    "\n",
    "    # Save to CSV\n",
    "    df.to_csv(output_csv, index=False)\n",
    "    print(f\"📂 Reflectance data saved to: {output_csv}\")\n",
    "\n",
    "\n",
    "    # # 🔥 PLOTTING: Raster + Overlapping Polygons\n",
    "    # print(\"📊 Plotting raster and polygons...\")\n",
    "\n",
    "    # fig, ax = plt.subplots(figsize=(10, 8))\n",
    "\n",
    "    # # Plot raster (use first band)\n",
    "    # with rasterio.open(raster_path) as src:\n",
    "    #     rasterio.plot.show(src.read(1), transform=src.transform, ax=ax, cmap=\"gray\")\n",
    "\n",
    "    # # Plot polygons\n",
    "    # polygons.boundary.plot(ax=ax, edgecolor=\"red\", linewidth=1.5)\n",
    "\n",
    "    # # Label polygons with IDs\n",
    "    # for idx, row in polygons.iterrows():\n",
    "    #     x, y = row.geometry.centroid.x, row.geometry.centroid.y\n",
    "    #     ax.text(x, y, str(row[\"OBJECTID\"]), fontsize=12, color=\"blue\", weight=\"bold\", ha=\"center\")\n",
    "\n",
    "    # ax.set_title(\"Raster with Overlapping Polygons\")\n",
    "    # plt.show()\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "fd7ff12e-7511-42a3-8a8e-8cfc39eca514",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "🔄 Reprojecting GeoJSON from EPSG:4326 to EPSG:26913\n"
     ]
    }
   ],
   "source": [
    "polygons = check_and_reproject(geojson_path, raster_path_2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "3d18c0d7-a0ef-4fa2-b600-cd72ede95dcf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "🔄 Reprojecting GeoJSON from EPSG:4326 to EPSG:26913\n",
      "✅ 58 polygons found within raster extent.\n",
      "📂 Reflectance data saved to: /home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/NEON_D13_NIWO_reflectance_per_pixel.csv\n"
     ]
    }
   ],
   "source": [
    "extract_pixel_reflectance(geojson_path, raster_path_2, output_csv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "5d89bd75-1bad-452b-826d-b27954ab6c3f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "🔄 Reprojecting GeoJSON from EPSG:5070 to EPSG:26913\n",
      "✅ 58 polygons found within raster extent.\n",
      "📂 Reflectance data saved to: /home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/NEON_D13_NIWO_reflectance_per_pixel_half.csv\n"
     ]
    }
   ],
   "source": [
    "extract_pixel_reflectance(geojson_path_half, raster_path_2, output_csv_half)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "a1a5a5c9-6910-469f-86f7-ff1405e76b53",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Copied: aligned_chm.tif\n",
      "Copied: aligned_slope.tif\n",
      "Copied: aligned_aspect.tif\n",
      "Copied: aligned_orthomosaic.tif\n",
      "Copied: NEON_D13_NIWO_aligned_orthomosaic.h5\n"
     ]
    },
    {
     "ename": "Error",
     "evalue": "[('/home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/.ipynb_checkpoints/NEON_D13_NIWO_reflectance_per_pixel-checkpoint.csv', '/data-store/iplant/home/shared/earthlab/macrosystems/Drone_topo_corr/corrected_new/NEON_D13_NIWO_aligned_orthomosaic/.ipynb_checkpoints/NEON_D13_NIWO_reflectance_per_pixel-checkpoint.csv', \"[Errno 121] Remote I/O error: '/data-store/iplant/home/shared/earthlab/macrosystems/Drone_topo_corr/corrected_new/NEON_D13_NIWO_aligned_orthomosaic/.ipynb_checkpoints/NEON_D13_NIWO_reflectance_per_pixel-checkpoint.csv'\")]",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mError\u001b[0m                                     Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[32], line 56\u001b[0m\n\u001b[1;32m     53\u001b[0m src_directory \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m/home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m     54\u001b[0m dest_directory \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m/data-store/iplant/home/shared/earthlab/macrosystems/Drone_topo_corr/corrected_new\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m---> 56\u001b[0m \u001b[43mcopy_directory\u001b[49m\u001b[43m(\u001b[49m\u001b[43msrc_directory\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdest_directory\u001b[49m\u001b[43m)\u001b[49m\n",
      "Cell \u001b[0;32mIn[32], line 48\u001b[0m, in \u001b[0;36mcopy_directory\u001b[0;34m(src_dir, dest_dir)\u001b[0m\n\u001b[1;32m     45\u001b[0m dest_dir_with_src_name \u001b[38;5;241m=\u001b[39m os\u001b[38;5;241m.\u001b[39mpath\u001b[38;5;241m.\u001b[39mjoin(dest_dir, os\u001b[38;5;241m.\u001b[39mpath\u001b[38;5;241m.\u001b[39mbasename(src_dir))\n\u001b[1;32m     47\u001b[0m \u001b[38;5;66;03m# Copy the source directory into the destination\u001b[39;00m\n\u001b[0;32m---> 48\u001b[0m \u001b[43mshutil\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcopytree\u001b[49m\u001b[43m(\u001b[49m\u001b[43msrc_dir\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdest_dir_with_src_name\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdirs_exist_ok\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m)\u001b[49m\n\u001b[1;32m     49\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mCopied directory: \u001b[39m\u001b[38;5;132;01m{\u001b[39;00msrc_dir\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m to \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mdest_dir_with_src_name\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m)\n",
      "File \u001b[0;32m/opt/conda/envs/macrosystems/lib/python3.10/shutil.py:559\u001b[0m, in \u001b[0;36mcopytree\u001b[0;34m(src, dst, symlinks, ignore, copy_function, ignore_dangling_symlinks, dirs_exist_ok)\u001b[0m\n\u001b[1;32m    557\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m os\u001b[38;5;241m.\u001b[39mscandir(src) \u001b[38;5;28;01mas\u001b[39;00m itr:\n\u001b[1;32m    558\u001b[0m     entries \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mlist\u001b[39m(itr)\n\u001b[0;32m--> 559\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43m_copytree\u001b[49m\u001b[43m(\u001b[49m\u001b[43mentries\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mentries\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43msrc\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43msrc\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdst\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdst\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43msymlinks\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43msymlinks\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    560\u001b[0m \u001b[43m                 \u001b[49m\u001b[43mignore\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mignore\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcopy_function\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mcopy_function\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    561\u001b[0m \u001b[43m                 \u001b[49m\u001b[43mignore_dangling_symlinks\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mignore_dangling_symlinks\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    562\u001b[0m \u001b[43m                 \u001b[49m\u001b[43mdirs_exist_ok\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdirs_exist_ok\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m/opt/conda/envs/macrosystems/lib/python3.10/shutil.py:513\u001b[0m, in \u001b[0;36m_copytree\u001b[0;34m(entries, src, dst, symlinks, ignore, copy_function, ignore_dangling_symlinks, dirs_exist_ok)\u001b[0m\n\u001b[1;32m    511\u001b[0m         errors\u001b[38;5;241m.\u001b[39mappend((src, dst, \u001b[38;5;28mstr\u001b[39m(why)))\n\u001b[1;32m    512\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m errors:\n\u001b[0;32m--> 513\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m Error(errors)\n\u001b[1;32m    514\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m dst\n",
      "\u001b[0;31mError\u001b[0m: [('/home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/.ipynb_checkpoints/NEON_D13_NIWO_reflectance_per_pixel-checkpoint.csv', '/data-store/iplant/home/shared/earthlab/macrosystems/Drone_topo_corr/corrected_new/NEON_D13_NIWO_aligned_orthomosaic/.ipynb_checkpoints/NEON_D13_NIWO_reflectance_per_pixel-checkpoint.csv', \"[Errno 121] Remote I/O error: '/data-store/iplant/home/shared/earthlab/macrosystems/Drone_topo_corr/corrected_new/NEON_D13_NIWO_aligned_orthomosaic/.ipynb_checkpoints/NEON_D13_NIWO_reflectance_per_pixel-checkpoint.csv'\")]"
     ]
    }
   ],
   "source": [
    "import shutil\n",
    "import os\n",
    "\n",
    "def copy_files(src_dir, dest_dir):\n",
    "    \"\"\"\n",
    "    Copies all files from src_dir to dest_dir.\n",
    "    \n",
    "    Parameters:\n",
    "        src_dir (str): Path to the source directory.\n",
    "        dest_dir (str): Path to the destination directory.\n",
    "    \"\"\"\n",
    "    # Ensure destination directory exists\n",
    "    if not os.path.exists(dest_dir):\n",
    "        os.makedirs(dest_dir)\n",
    "\n",
    "    # Copy each file from source to destination\n",
    "    for filename in os.listdir(src_dir):\n",
    "        src_path = os.path.join(src_dir, filename)\n",
    "        dest_path = os.path.join(dest_dir, filename)\n",
    "        \n",
    "        # Copy only if it's a file\n",
    "        if os.path.isfile(src_path):\n",
    "            shutil.copy(src_path, dest_path)\n",
    "            print(f\"Copied: {filename}\")\n",
    "\n",
    "# Example usage\n",
    "src_directory = \"/home/jovyan/data-store/cross-sensor-cal/data_set2\"\n",
    "dest_directory = \"/home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic\"\n",
    "\n",
    "copy_files(src_directory, dest_directory)\n",
    "\n",
    "def copy_directory(src_dir, dest_dir):\n",
    "    \"\"\"\n",
    "    Copies an entire directory (and its contents) into another directory.\n",
    "\n",
    "    Parameters:\n",
    "        src_dir (str): Path to the source directory.\n",
    "        dest_dir (str): Path to the destination directory.\n",
    "    \"\"\"\n",
    "    # Ensure the destination directory exists\n",
    "    if not os.path.exists(dest_dir):\n",
    "        os.makedirs(dest_dir)\n",
    "\n",
    "    # Add the source directory name to the destination path\n",
    "    dest_dir_with_src_name = os.path.join(dest_dir, os.path.basename(src_dir))\n",
    "\n",
    "    # Copy the source directory into the destination\n",
    "    shutil.copytree(src_dir, dest_dir_with_src_name, dirs_exist_ok=True)\n",
    "    print(f\"Copied directory: {src_dir} to {dest_dir_with_src_name}\")\n",
    "\n",
    "\n",
    "# Example usage\n",
    "src_directory = \"/home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic\"\n",
    "dest_directory = \"/data-store/iplant/home/shared/earthlab/macrosystems/Drone_topo_corr/corrected_new\"\n",
    "\n",
    "copy_directory(src_directory, dest_directory)\n"
   ]
  },
  {
   "cell_type": "raw",
   "id": "628daed2-d48c-4a91-92ae-60c22b0fbd7b",
   "metadata": {},
   "source": [
    "import hytools as ht\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os"
   ]
  },
  {
   "cell_type": "raw",
   "id": "f5227c80-df12-446d-966a-267cf310bc4a",
   "metadata": {},
   "source": [
    "os.getcwd()"
   ]
  },
  {
   "cell_type": "raw",
   "id": "8172fae3-67a8-4d2a-8385-7e2a46cccb63",
   "metadata": {},
   "source": [
    "envi_image= '/home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/NEON_D13_NIWO_aligned_orthomosaic__reflectance_corrected'\n",
    "envi = ht.HyTools()\n",
    "envi.read_file(envi_image,'envi')\n"
   ]
  },
  {
   "cell_type": "raw",
   "id": "7b64df5a-86ac-45b6-8d21-13d13eb33dfe",
   "metadata": {},
   "source": [
    "original_image= '/home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/NEON_D13_NIWO_aligned_orthomosaic'\n",
    "\n",
    "og_envi = ht.HyTools()\n",
    "og_envi.read_file(original_image,'envi')"
   ]
  },
  {
   "cell_type": "raw",
   "id": "49d9099b-a36c-4f1b-aa11-4040a5c6bf1c",
   "metadata": {},
   "source": [
    "line = envi.get_line(100, corrections = [])\n",
    "band = envi.get_band(5)\n",
    "wave = envi.get_wave(668)\n",
    "pixel = envi.get_pixels([100],[200])\n",
    "#anc = envi.get_anc('sensor_zn',radians =False)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "4db76b42-38e9-430f-b8ef-e0a351417c44",
   "metadata": {},
   "source": [
    "topo_coeffs = '/home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/NEON_D13_NIWO_aligned_orthomosaic_topo_coeffs__reflectance.json'\n",
    "envi.load_coeffs(topo_coeffs,'topo')\n",
    "\n",
    "brdf_coeffs = '/home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/NEON_D13_NIWO_aligned_orthomosaic_brdf_coeffs__reflectance.json'\n",
    "envi.load_coeffs(brdf_coeffs,'brdf')"
   ]
  },
  {
   "cell_type": "raw",
   "id": "b1a2632a-7f4f-4c38-a1ce-3985c8d1413a",
   "metadata": {},
   "source": [
    "def show_rgb(hy_obj,r=668,g=560,b=444, correct= []):\n",
    "\n",
    "    rgb=  np.stack([hy_obj.get_wave(r,corrections= correct),\n",
    "                    hy_obj.get_wave(g,corrections= correct),\n",
    "                    hy_obj.get_wave(b,corrections= correct)])\n",
    "    rgb = np.moveaxis(rgb,0,-1).astype(float)\n",
    "    rgb[rgb ==hy_obj.no_data] = np.nan\n",
    "\n",
    "    bottom = np.nanpercentile(rgb,5,axis = (0,1))\n",
    "    top = np.nanpercentile(rgb,95,axis = (0,1))\n",
    "    rgb = np.clip(rgb,bottom,top)\n",
    "\n",
    "    rgb = (rgb-np.nanmin(rgb,axis=(0,1)))/(np.nanmax(rgb,axis= (0,1))-np.nanmin(rgb,axis= (0,1)))\n",
    "\n",
    "    height = int(hy_obj.lines/hy_obj.columns)\n",
    "\n",
    "    fig  = plt.figure(figsize = (7,7) )\n",
    "    plt.imshow(rgb)\n",
    "    plt.show()\n",
    "    plt.close()\n",
    "    \n"
   ]
  },
  {
   "cell_type": "raw",
   "id": "431a175b-3540-4047-9e71-87bc06f9c3ef",
   "metadata": {},
   "source": [
    "show_rgb(envi, correct= [])"
   ]
  },
  {
   "cell_type": "raw",
   "id": "a9584149-1887-4621-a4f4-a1c0c7e0ba8f",
   "metadata": {},
   "source": [
    "show_rgb(og_envi)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "5116a817-30d4-45a0-be02-5e0f2d3660c8",
   "metadata": {},
   "source": [
    "show_rgb(og_envi, correct = ['brdf'])"
   ]
  },
  {
   "cell_type": "raw",
   "id": "493c145c-b558-4229-a5c8-97c9baed1bde",
   "metadata": {},
   "source": [
    "og_envi.get_wave(862, corrections= ['topo', 'brdf'])"
   ]
  },
  {
   "cell_type": "raw",
   "id": "614cc82f-c2fb-4c56-84ab-7d037d7d2de5",
   "metadata": {},
   "source": [
    "difference = (envi.get_wave(862)) - og_envi.get_wave(862)\n",
    "plt.matshow(difference, cmap='viridis')  # 'viridis' is the default colormap in Matplotlib\n",
    "\n",
    "# Add a colorbar to act as a legend for the values\n",
    "cbar = plt.colorbar()\n",
    "cbar.set_label('Difference Values', fontsize=12)  # Label for the colorbar\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "raw",
   "id": "f25c7f73-e399-45f4-8cfb-bbb96e9a8393",
   "metadata": {},
   "source": [
    "envi.wavelengths"
   ]
  },
  {
   "cell_type": "raw",
   "id": "a5a320a4-e7c7-41e1-a2db-05e2939a86f9",
   "metadata": {},
   "source": [
    "corrected = envi.get_pixels([200],[300])\n",
    "corrected\n"
   ]
  },
  {
   "cell_type": "raw",
   "id": "5948f2dc-dbd3-450c-b35a-4820d7e73a19",
   "metadata": {},
   "source": [
    "reflectance = corrected.flatten()\n",
    "reflectance"
   ]
  },
  {
   "cell_type": "raw",
   "id": "3706570f-1f33-4000-ba9b-7fa138c0cffb",
   "metadata": {},
   "source": [
    "og_envi2 = og_envi.get_pixels([200],[300])\n",
    "og_envi2"
   ]
  },
  {
   "cell_type": "raw",
   "id": "f995474b-b6cc-40d5-af32-a0f4de37246b",
   "metadata": {},
   "source": [
    "old_ref = og_envi2.flatten()\n",
    "old_ref"
   ]
  },
  {
   "cell_type": "raw",
   "id": "66483ef3-57b1-4b6c-87e0-206cf649df47",
   "metadata": {},
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Create the figure and axes\n",
    "plt.figure(figsize=(8, 5))  # Set figure size\n",
    "\n",
    "# Plot both reflectance curves on the same graph\n",
    "plt.plot(envi.wavelengths, reflectance, marker='o', linestyle='-', color='b', label='Corrected Reflectance')\n",
    "plt.plot(envi.wavelengths, old_ref, marker='s', linestyle='--', color='r', label='Original Reflectance')\n",
    "\n",
    "# Labels and title\n",
    "plt.xlabel(\"Wavelength (nm)\", fontsize=12)\n",
    "plt.ylabel(\"Reflectance\", fontsize=12)\n",
    "plt.title(\"Reflectance vs. Wavelength\", fontsize=14)\n",
    "\n",
    "# Legend and grid\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "\n",
    "# Show the plot\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "raw",
   "id": "cddc84e2-35ad-4d6f-a2e4-ffc8e2a44e32",
   "metadata": {},
   "source": [
    "from spectral import open_image\n",
    "\n",
    "def read_envi_metadata(envi_path):\n",
    "    \"\"\"\n",
    "    Reads metadata from an ENVI raster file and prints key attributes.\n",
    "    \"\"\"\n",
    "    # Open the ENVI file\n",
    "    envi_image = open_image(envi_path)\n",
    "\n",
    "    # Print all metadata\n",
    "    print(\"ENVI Metadata:\")\n",
    "    for key, value in envi_image.metadata.items():\n",
    "        print(f\"{key}: {value}\")\n",
    "\n",
    "    # Extract wavelengths if available\n",
    "    if 'wavelength' in envi_image.metadata:\n",
    "        wavelengths = envi_image.metadata['wavelength']\n",
    "        print(\"\\nWavelengths:\")\n",
    "        print(wavelengths)\n",
    "\n",
    "# Example usage\n",
    "\n",
    "#hdr_path_2 = os.path.join(directory, main_image_name + '__reflectance_corrected_flipped.hdr' )\n",
    "hdr_path_2 = '/home/jovyan/data-store/cross-sensor-cal/output/NEON_D13_NIWO_aligned_orthomosaic/NEON_D13_NIWO_aligned_orthomosaic__reflectance_corrected_flipped.hdr'\n",
    "read_envi_metadata(hdr_path_2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6738ffa-228e-47a6-8680-f9146b7c63e7",
   "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
}