README: Thermal Correction for Hyperspectral Imaging Data

Overview

This Python script performs thermal correction on hyperspectral imaging (HSI) data, retrieving emissivity, surface temperature, and disk function values for each pixel. It utilizes Planck's Law and a forward modeling approach to separate reflected and emitted radiance, optimizing these parameters using observed radiance values.

---

Features

1. Planck's Law Implementation: Calculates spectral radiance based on wavelength and temperature.
2. Forward Model: Models radiance as a combination of emitted and reflected components.
3. Optimization: Uses the scipy.optimize.least_squares function to fit the model to observed radiance.
4. Solar Flux Interpolation: Adapts solar flux data to match the wavelengths of the HSI image.
5. Window-Based Processing: Processes a sliding 3x3 pixel window for improved stability in optimization.
6. Outputs:
   • Emissivity map (H x W x Bands)
   • Surface temperature map (H x W)
   • Disk function map (H x W)
   • Residual cost map (H x W)

---

Requirements

Libraries

• Python 3.8+
• NumPy
• SciPy
• Rasterio
• Matplotlib (optional for visualization)

Input Files

1. Hyperspectral Image (GeoTIFF): Contains radiance data with dimensions (Height, Width, Bands).
2. Solar Flux File (Text): A text file with two columns:
   • Wavelength (in micrometers or nanometers)
   • Flux values

---

Workflow

1. Load Input Data

• The script reads the HSI image and extracts the radiance data. The dimensions are transposed to match (Height, Width, Bands).
• Reads the solar flux file and interpolates it to match the HSI image wavelengths.

2. Initialize Parameters

• Emissivity: Defaulted to 0.8 for all bands.
• Surface Temperature: Defaulted to 300 K.
• Disk Function: Defaulted to 1.0.
• Parameter bounds:
  • Emissivity: [0, 1]
  • Temperature: [250 K, 400 K]
  • Disk Function: [0.5, 2.0]

3. Perform Thermal Correction

• The script processes the image using a 3x3 sliding window approach:
  • Observed radiance for the window is reshaped for optimization.
  • The forward model calculates modeled radiance using the current state vector.
  • Residuals are minimized to fit the observed radiance.

4. Outputs

• Emissivity Map: (H, W, Bands) array of retrieved emissivity values.
• Temperature Map: (H, W) array of surface temperatures.
• Disk Function Map: (H, W) array of disk function values.
• Residual Costs: (H, W) array showing optimization costs for each pixel.

---

How to Run

1. Update file paths:

   • Replace clipped_img.tif with the path to your HSI image.
   • Replace ch2_iirs_solar_flux.txt with the path to your solar flux file.

2. Run the script:

   python thermal_correction.py

3. Outputs:

   • The script will print progress during processing.
   • Outputs can be saved to disk or visualized.

---

Visualization (Optional)

You can visualize the results using libraries like Matplotlib:

Residual Plot for a Pixel

import matplotlib.pyplot as plt

# Example for pixel (50, 50)
pixel_index = (50, 50)
observed = image[:, pixel_index[0], pixel_index[1]]
modeled = forward_model(
    [output_emissivity[pixel_index[0], pixel_index[1], :], 
     output_temperature[pixel_index[0], pixel_index[1]], 
     output_disk_function[pixel_index[0], pixel_index[1]]], 
    image_wavelengths, 
    J_lambda
)
residuals = observed - modeled

plt.plot(image_wavelengths, residuals, marker='o')
plt.axhline(0, color='r', linestyle='--')
plt.xlabel("Wavelength (nm)")
plt.ylabel("Residual (Observed - Modeled)")
plt.title("Residual Plot for Pixel (50, 50)")
plt.show()

---

Key Functions

1. planck_law(wavelength, T):
   Calculates spectral radiance for a given wavelength and temperature using Planck's Law.

2. forward_model(state, wavelengths, J_lambda):
   Models the total radiance as the sum of reflected and emitted radiance.

3. cost_function(state, wavelengths, J_lambda, observed_radiance):
   Computes residuals for optimization.

4. process_window(window_radiance, wavelengths, J_lambda, init_state, bounds):
   Performs optimization for a single 3x3 pixel window.

5. process_image(image, wavelengths, J_lambda, init_state, bounds):
   Processes the entire image using a sliding window approach.

---

Notes

• Ensure that the solar flux file covers the wavelength range of your HSI image.
• Use appropriate initial guesses and bounds for the optimization process.
• The 3x3 window approach helps stabilize optimization but may slightly increase computation time.

For any questions or troubleshooting, feel free to ask!