|
| 1 | +# icesat-2_kdph_py |
| 2 | + |
| 3 | + |
| 4 | +This repository contains Python code developed in collaboration to replicate the functionality of the original MATLAB scripts (https://github.com/emilyeidam/icesat-2_kdph |
| 5 | +) by Dr. Emily Eidam (emily.eidam@oregonstate.edu). The code facilitates processing ICESat-2 data to calculate the diffuse attenuation coefficient (Kd) based on space-based lidar photon profiles, following methods described in the original MATLAB code. Please cite appropriately both Python and the MATLAB version under the GNU GPLv3 license if you use this code in your research. |
| 6 | + |
| 7 | +## Repository Structure |
| 8 | + |
| 9 | +Current structure: |
| 10 | + |
| 11 | +```text |
| 12 | +icesat-2_kdph_py/ |
| 13 | ++-- main.py |
| 14 | ++-- config.py |
| 15 | ++-- requirements.txt |
| 16 | ++-- README.md |
| 17 | ++-- kd_utils/ |
| 18 | + +-- __init__.py |
| 19 | + +-- data_processing.py |
| 20 | + +-- sea_photons_analysis.py |
| 21 | + +-- bathy_processing.py |
| 22 | + +-- Kd_analysis.py |
| 23 | + +-- interpolation.py |
| 24 | + +-- visualization.py |
| 25 | + +-- sliderule_adapter.py |
| 26 | + +-- SeaSurfaceFlattening.py |
| 27 | + +-- SolarBckgrd/ |
| 28 | +``` |
| 29 | + |
| 30 | +Key modules: |
| 31 | + |
| 32 | +- `main.py`: End-to-end pipeline entrypoint for local ATL03 workflows. |
| 33 | +- `config.py`: All CLI switches, paths, and thresholds. |
| 34 | +- `kd_utils/data_processing.py`: ATL03 ingestion, land/sea masking, optional pre-filters. |
| 35 | +- `kd_utils/sea_photons_analysis.py`: Binning and sea-surface detection. |
| 36 | +- `kd_utils/bathy_processing.py`: Subsurface filtering and optional gap steps (GEBCO/ATL24/refraction/etc.). |
| 37 | +- `kd_utils/Kd_analysis.py`: Kd fitting. |
| 38 | +- `kd_utils/sliderule_adapter.py`: SlideRule ATL03 fetch + conversion to framework schema for notebook/script imports. |
| 39 | + |
| 40 | +### Getting Started |
| 41 | + |
| 42 | +### Prerequisites |
| 43 | + |
| 44 | +- Python 3.x |
| 45 | +- Libraries: `numpy`, `pandas`, `scipy`, `h5py`, `matplotlib` |
| 46 | +- Ensure you have ICESat-2 ATL03 `.h5` files available in the appropriate directory, as specified in `config.py`, and also refers to the original MATLAB code. |
| 47 | + |
| 48 | +### Installation |
| 49 | + |
| 50 | +1. Clone the repository: |
| 51 | + |
| 52 | + ```bash |
| 53 | + git clone https://github.com/ChaoEcohydroRS/icesat-2_kdph_py.git |
| 54 | + cd IS2_Kd_Proj |
| 55 | + ``` |
| 56 | + |
| 57 | + |
| 58 | +2. Install the required packages: |
| 59 | + ``` |
| 60 | + pip install -r requirements.txt |
| 61 | + ``` |
| 62 | + |
| 63 | +3. Ensure .h5 files are placed in the appropriate default folder structure. |
| 64 | + |
| 65 | +### Usage |
| 66 | +1. Edit the config.py file to set up file paths and any parameters. |
| 67 | +2. Run the entire workflow through main.py: |
| 68 | + |
| 69 | + ```bash |
| 70 | + python main.py |
| 71 | + ``` |
| 72 | + |
| 73 | +### SlideRule Ingestion (Notebook Import Path) |
| 74 | + |
| 75 | +You can reuse the same framework without local ATL03 `.h5` files by importing the SlideRule adapter: |
| 76 | + |
| 77 | +```python |
| 78 | +from kd_utils.sliderule_adapter import fetch_and_prepare_sliderule_dataset |
| 79 | +from kd_utils.sea_photons_analysis import process_sea_photon_binning |
| 80 | +from kd_utils.bathy_processing import process_subsurface_photon_filtering |
| 81 | +from kd_utils.Kd_analysis import process_kd_calculation |
| 82 | + |
| 83 | +region = [ |
| 84 | + {"lon": -84.029, "lat": 29.732}, |
| 85 | + {"lon": -84.029, "lat": 29.954}, |
| 86 | + {"lon": -83.969, "lat": 29.954}, |
| 87 | + {"lon": -83.969, "lat": 29.732}, |
| 88 | + {"lon": -84.029, "lat": 29.732}, |
| 89 | +] |
| 90 | + |
| 91 | +sea_photon_dataset = fetch_and_prepare_sliderule_dataset( |
| 92 | + region=region, |
| 93 | + t0="2025-05-22T00:00:00Z", |
| 94 | + t1="2025-05-24T00:00:00Z", |
| 95 | + rgt=1033, |
| 96 | + strong_beams_only=True, |
| 97 | +) |
| 98 | + |
| 99 | +binned = process_sea_photon_binning(sea_photon_dataset, horizontal_res=500, vertical_res=0.25) |
| 100 | +sea_h, sea_lbl, subsurface = process_subsurface_photon_filtering( |
| 101 | + binned, GEBCO_paths=[], subsurface_thresh=1.0, Ignore_Subsurface_Height_Thres=-6 |
| 102 | +) |
| 103 | +kd_df = process_kd_calculation(subsurface) |
| 104 | +``` |
| 105 | + |
| 106 | +This path reuses your existing processing framework; only ingestion switches from local HDF5 to SlideRule API. |
| 107 | + |
| 108 | +### Optional Gap Steps (Sensitivity On/Off) |
| 109 | + |
| 110 | +All extra processing steps are optional and disabled by default. |
| 111 | + |
| 112 | +- IR/AP proxy photon filtering: |
| 113 | + ```bash |
| 114 | + python main.py --enable_ir_ap_filter |
| 115 | + ``` |
| 116 | +- Solar background filtering: |
| 117 | + ```bash |
| 118 | + python main.py --enable_solar_background_filter |
| 119 | + ``` |
| 120 | +- GEBCO seafloor filtering: |
| 121 | + ```bash |
| 122 | + python main.py --enable_gebco_filter |
| 123 | + ``` |
| 124 | +- ATL24-first with GEBCO fallback: |
| 125 | + ```bash |
| 126 | + python main.py --enable_atl24_filter --atl24_file path/to/atl24_points.csv --enable_gebco_filter |
| 127 | + ``` |
| 128 | +- Sea surface flattening: |
| 129 | + ```bash |
| 130 | + python main.py --enable_sea_surface_flattening |
| 131 | + ``` |
| 132 | +- Convex hull reasonableness filter: |
| 133 | + ```bash |
| 134 | + python main.py --enable_convex_hull_filter |
| 135 | + ``` |
| 136 | +- Paired-beam combine before Kd fit: |
| 137 | + ```bash |
| 138 | + python main.py --enable_paired_beam_combine |
| 139 | + ``` |
| 140 | +- Histogram quality filter (<5% signal at 6-7 m depth band): |
| 141 | + ```bash |
| 142 | + python main.py --enable_histogram_quality_filter |
| 143 | + ``` |
| 144 | +- Surface Gaussian sigma filter: |
| 145 | + ```bash |
| 146 | + python main.py --enable_surface_sigma_filter |
| 147 | + ``` |
| 148 | +- Refraction correction: |
| 149 | + ```bash |
| 150 | + python main.py --enable_refraction_correction |
| 151 | + ``` |
| 152 | +- Post-refraction rebuild + re-fit: |
| 153 | + ```bash |
| 154 | + python main.py --enable_refraction_correction --enable_post_refraction_refit |
| 155 | + ``` |
| 156 | + |
| 157 | +### Gap Implementation Progress |
| 158 | + |
| 159 | +The flowchart gaps are being implemented one-by-one with optional switches. |
| 160 | + |
| 161 | +- [x] Gap 1: Remove photons flagged as IR/AP (proxy implementation via `quality_ph` and `photon_conf`) |
| 162 | + - Switch: `--enable_ir_ap_filter` |
| 163 | + - Tuning: `--ir_ap_quality_max`, `--ir_ap_min_signal_conf` |
| 164 | +- [x] Gap 2: Histogram quality review (<5% signal at 6-7 m) |
| 165 | + - Switch: `--enable_histogram_quality_filter` |
| 166 | + - Tuning: |
| 167 | + - `--histogram_quality_min_ratio` (default `0.05`) |
| 168 | + - `--histogram_quality_depth_min` (default `6.0`) |
| 169 | + - `--histogram_quality_depth_max` (default `7.0`) |
| 170 | +- [x] Gap 3: Discard bins with high surface Gaussian std dev |
| 171 | + - Switch: `--enable_surface_sigma_filter` |
| 172 | + - Tuning: `--surface_sigma_max` (default `0.5`) |
| 173 | +- [x] Gap 4: Refraction correction integration |
| 174 | + - Switch: `--enable_refraction_correction` |
| 175 | + - Tuning: |
| 176 | + - `--refraction_water_temp_c` (default `20.0`) |
| 177 | + - `--refraction_wavelength_nm` (default `532.0`) |
| 178 | +- [x] Gap 5: Rebuild histogram and re-fit after refraction |
| 179 | + - Switch: `--enable_post_refraction_refit` (requires `--enable_refraction_correction`) |
| 180 | +- [x] Gap 6: ATL24 query + fallback orchestration |
| 181 | + - Switch: `--enable_atl24_filter` |
| 182 | + - Note: current implementation uses local ATL24 point files provided via `--atl24_file`. |
| 183 | + - Inputs: |
| 184 | + - `--atl24_file` (CSV/Parquet/GPKG/SHP with lon/lat + bathymetry field) |
| 185 | + - `--atl24_max_match_distance_deg` (default `0.01`) |
| 186 | + - Fallback: if ATL24 has no usable matches and `--enable_gebco_filter` is on, GEBCO filtering is used. |
| 187 | +- [x] Gap 7: Combine paired beams |
| 188 | + - Switch: `--enable_paired_beam_combine` |
| 189 | + - Method: maps left/right beams to `gt1_pair`, `gt2_pair`, `gt3_pair` and re-bins by `relative_AT_dist` using `horizontal_res`. |
| 190 | + |
| 191 | +### Citation |
| 192 | +If you use this code, please cite: |
| 193 | + |
| 194 | +Eidam, E.F., K. Bisson, C. Wang, C. Walker, and A. Gibbons (2024). ICESat-2 and ocean particulates: A roadmap for calculating Kd from space-based lidar photon profiles. Remote Sensing of Environment. Vol. 311. https://doi.org/10.1016/j.rse.2024.114222 |
| 195 | + |
| 196 | +### License |
| 197 | +This project is licensed under the GNU GPLv3 license. |
0 commit comments