diff --git a/noa-coherence/README.md b/noa-coherence/README.md new file mode 100644 index 0000000..7cba61b --- /dev/null +++ b/noa-coherence/README.md @@ -0,0 +1,243 @@ +# Sentinel-1 InSAR Coherence Processing Pipeline + +## Table of Contents + +- [Introduction](#introduction) +- [What is Coherence?](#what-is-coherence) +- [Prerequisites](#prerequisites) +- [Installation](#installation) +- [Usage](#usage) + - [Command-Line Arguments](#command-line-arguments) + - [Example](#example) +- [Functionality](#functionality) +- [Contributing](#contributing) +- [License](#license) + +## Introduction + +The **Sentinel-1 InSAR Coherence Processing Pipeline** is a Python-based tool designed to process Sentinel-1 Single Look Complex (SLC) products using the [SNAP](http://step.esa.int/main/download/) (Sentinel Application Platform) framework via the SNAPPY python library. This pipeline facilitates the generation of coherence products, which are essential for various interferometric analyses, including surface deformation studies and terrain mapping. + +## What is Coherence? + +In the context of Interferometric Synthetic Aperture Radar (InSAR), **coherence** is a measure of the similarity between two SAR images taken from slightly different positions or at different times. It quantifies the degree of correlation between the phase and amplitude of the backscattered radar signals from two Synthetic Aperture Radar (SAR) images. High coherence values indicate stable scattering surfaces with minimal changes between acquisitions, while low coherence suggests significant changes or noise. + +Coherence is crucial for: + +- **Interferogram Quality:** High coherence ensures reliable phase information for interferometric measurements. +- **Surface Deformation Analysis:** Identifying areas with consistent scattering for accurate displacement measurements. +- **Terrain Mapping:** Enhancing the precision of digital elevation models (DEMs). + +## Prerequisites + +Before running the InSAR Coherence Processing Pipeline, ensure that the following software and dependencies are installed: + +- **Python 3.9 or higher** +- **SNAP (Sentinel Application Platform) 8.x or higher** +- **Java Runtime Environment (JRE) 8 or higher** +- **Required Python Libraries:** + - `jpy` + - `snappy` + - `argparse` + - `sys` + +### Installing SNAP + +Download and install SNAP from the [official ESA website](http://step.esa.int/main/download/). Follow the installation instructions specific to your operating system. + +### Setting Up Python Environment + +It's recommended to use a virtual environment to manage dependencies: + +```bash +# Create a virtual environment +python3 -m venv insar_env + +# Activate the virtual environment +# On Windows: +insar_env\Scripts\activate +# On macOS/Linux: +source insar_env/bin/activate + +# Upgrade pip +pip install --upgrade pip + +# Install required Python libraries +pip install jpy snappy argparse +``` + +> **Note:** The `snappy` library requires SNAP to be correctly installed and configured. Ensure that the `SNAP_HOME` environment variable points to your SNAP installation directory. + +## Installation + +1. **Clone the Repository:** + + ```bash + git clone https://github.com/yourusername/sentinel1-insar-pipeline.git + cd sentinel1-insar-pipeline + ``` + +2. **Install Dependencies:** + + Ensure that all prerequisites are met as outlined above. + +3. **Configure SNAP:** + + Make sure that SNAP is properly installed and that the `snappy` module is correctly configured. You may need to run the following command to set up `snappy`: + + ```bash + # Navigate to the snappy folder + cd $SNAP_HOME/snappy + + # Install snappy + python setup.py install + ``` + +## Usage + +The InSAR Coherence Processing Pipeline requires two Sentinel-1 SLC products (master and slave) to generate coherence products. The script accepts several command-line arguments to customize the processing parameters. + +### Command-Line Arguments + +```bash +usage: insar_pipeline.py [-h] master slave outpath minX minY maxX maxY +``` + +**Positional Arguments:** + +1. **master** (`str`): + Path to the master Sentinel-1 SLC product file. + +2. **slave** (`str`): + Path to the slave Sentinel-1 SLC product file. + +3. **outpath** (`str`): + Directory path where the output products will be saved. + +4. **minX** (`float`): + Minimum longitude of the bounding box (in decimal degrees). + +5. **minY** (`float`): + Minimum latitude of the bounding box (in decimal degrees). + +6. **maxX** (`float`): + Maximum longitude of the bounding box (in decimal degrees). + +7. **maxY** (`float`): + Maximum latitude of the bounding box (in decimal degrees). + +### Example + +```bash +python insar_pipeline.py \ + /path/to/master_SLC.dim \ + /path/to/slave_SLC.dim \ + /path/to/output_directory \ + -123.45 37.77 -122.45 38.77 \ +``` + +**Explanation:** + +- **master_SLC.dim**: Path to the master Sentinel-1 SLC product. +- **slave_SLC.dim**: Path to the slave Sentinel-1 SLC product. +- **EVT12345**: Event identifier. +- **/path/to/output_directory**: Directory where outputs will be saved. +- **-123.45 37.77 -122.45 38.77**: Bounding box coordinates (minX, minY, maxX, maxY). +- **additional_parameter_value**: Placeholder for any additional parameter required by the pipeline. + +## Functionality + +The pipeline consists of several functions that handle different processing steps. Below is an overview of the key functions: + +### 1. `read(filename: str) -> Product` + +Reads a Sentinel-1 product from the specified filename. + +### 2. `write(product, filename: str)` + +Writes a product to a file in GeoTIFF format. + +### 3. `topsar_split(product, swath: str, pol: str) -> Product` + +Splits a TOPSAR product into subswaths and polarizations. + +### 4. `topsar_merge(product, pol: str) -> Product` + +Merges split TOPSAR products based on polarization. + +### 5. `apply_orbit_file(product, mode: int) -> Product` + +Applies orbit file corrections to the product. +**Parameters:** +- `mode`: `1` for Restituted orbit, else Precise orbit. + +### 6. `back_geocoding(product) -> Product` + +Performs back-geocoding on the product using a Digital Elevation Model (DEM). + +### 7. `topsar_deburst(product, pol: str) -> Product` + +Removes burst overlaps from TOPSAR products based on polarization. + +### 8. `topophase_removal(product) -> Product` + +Removes topographic phase from the interferogram. + +### 9. `coherence_generation(product) -> Product` + +Generates coherence from the product. + +### 10. `geometric_correction(product, polarisation: str, to_print: bool = True) -> Product` + +Performs geometric correction on the product. + +### 11. `subset(product, bbox: list, to_print: bool = True) -> Optional[Product]` + +Subsets the product using a bounding box. + +### 12. `insar_pipeline(master: str, slave: str, outpath: str, event_id: str, bbox: list, other_param)` + +Main InSAR processing pipeline that orchestrates the entire workflow using the above functions. + +## Contributing + +Contributions are welcome! Please follow these steps to contribute: + +1. **Fork the Repository** + +2. **Create a New Branch** + + ```bash + git checkout -b feature/YourFeatureName + ``` + +3. **Commit Your Changes** + + Follow the [Conventional Commits](https://www.conventionalcommits.org/en/v1.0.0/) specification for commit messages. + +4. **Push to the Branch** + + ```bash + git push origin feature/YourFeatureName + ``` + +5. **Open a Pull Request** + + Provide a clear description of your changes and reference any related issues. + +## License + +This project is licensed under the [MIT License](LICENSE). + +--- + +**Note:** Ensure that all paths and parameters are correctly specified when running the script. For any issues or feature requests, please open an [issue](https://github.com/yourusername/sentinel1-insar-pipeline/issues) on GitHub. + +--- + +### Additional Tips: + +- **Verify Markdown Rendering:** After pasting the content into your `README.md`, preview it on GitHub or your preferred Markdown viewer to ensure that all sections, links, and code blocks render correctly. +- **Update Repository Links:** Replace `https://github.com/yourusername/sentinel1-insar-pipeline.git` and other placeholder links with your actual repository URLs. +- **Customize as Needed:** Feel free to modify the content to better fit your project's specifics, such as adding more detailed usage examples, screenshots, or additional sections relevant to your users. + +By following these steps, your `README.md` should display properly and provide comprehensive information to users and contributors of your project. diff --git a/noa-coherence/coherence.py b/noa-coherence/coherence.py new file mode 100644 index 0000000..5761ecf --- /dev/null +++ b/noa-coherence/coherence.py @@ -0,0 +1,329 @@ +import os,sys +import datetime as dt +import argparse + +# Change the path below +sys.path.append('/home//.snap/snap-python') + +from snappy import GPF +from snappy import ProductIO +from snappy import HashMap +from snappy import jpy +from snappy import WKTReader +from snappy import File +from snappy import ProgressMonitor +from time import sleep +import gc + + +def read(filename): + reader = ProductIO.getProductReader('SENTINEL-1') + + return reader.readProductNodes(filename,None) + +def write(product, filename): + ProductIO.writeProduct(product, filename, "GeoTIFF") + + +def topsar_split(product, swath: str, pol: str): + """ + Splits a TOPSAR product into subswaths and polarizations. + + Parameters: + product: The input TOPSAR product. + swath (str): The subswath to process (e.g., 'IW1', 'IW2', 'IW3'). + pol (str): The polarization to select (e.g., 'VV', 'VH'). + + Returns: + Product: The split product. + """ + parameters = HashMap() + parameters.put('subswath', swath) + parameters.put('selectedPolarisations', pol) + + return GPF.createProduct("TOPSAR-Split", parameters, product) + +def topsar_merge(product, pol: str): + """ + Merges split TOPSAR products. + + Parameters: + product: The input product. + pol (str): The polarization to select. + + Returns: + Product: The merged product. + """ + parameters = HashMap() + parameters.put('selectedPolarisations', pol) + + return GPF.createProduct("TOPSAR-Merge", parameters, product) + + +def apply_orbit_file(product, mode: int): + """ + Applies orbit file corrections to the product. + + Parameters: + product: The input product. + mode (int): Orbit file mode (1 for Restituted, else Precise). + + Returns: + Product: The orbit-corrected product. + """ + parameters = HashMap() + + if mode == 1: + parameters.put("Orbit State Vectors", "Sentinel Restituted (Auto Download)") + else: + parameters.put("Orbit State Vectors", "Sentinel Precise (Auto Download)") + + parameters.put("Polynomial Degree", 3) + + return GPF.createProduct("Apply-Orbit-File", parameters, product) + +def back_geocoding(product): + """ + Performs back-geocoding on the product. + + Parameters: + product: The input product. + + Returns: + Product: The back-geocoded product. + """ + parameters = HashMap() + parameters.put("Digital Elevation Model", "SRTM 1Sec HGT (Auto Download)") + parameters.put("DEM Resampling Method", "NEAREST_NEIGHBOUR") + parameters.put("Resampling Type", "BICUBIC_INTERPOLATION") + parameters.put("Mask out areas with no elevation", True) + parameters.put("Output Deramp and Demod Phase", False) + + return GPF.createProduct("Back-Geocoding", parameters, product) + +def topsar_deburst(product, pol: str): + """ + Removes burst overlaps from TOPSAR products. + + Parameters: + product: The input product. + pol (str): The polarization to process. + + Returns: + Product: The debursted product. + """ + parameters = HashMap() + parameters.put("Polarisations", pol) + return GPF.createProduct("TOPSAR-Deburst", parameters, product) + +def coherence_generation(product): + """ + Generates coherence from the product. + + Parameters: + product: The input product. + + Returns: + Product: The coherence product. + """ + parameters = HashMap() + parameters.put("Subtract flat-earth phase", True) + parameters.put("Degree of \"Flat Earth\" polynomial", 5) + parameters.put("Number of \"Flat Earth\" estimation points", 501) + parameters.put("Orbit interpolation degree", 3) + parameters.put("Include coherence estimation", True) + parameters.put("Square Pixel", True) + parameters.put("Independent Window Sizes", False) + parameters.put("Coherence Azimuth Window Size", 3) + parameters.put("Coherence Range Window Size", 10) + return GPF.createProduct("Coherence", parameters, product) + +def write_snaphu(product, filename): + ProductIO.writeProduct(product, filename, "Snaphu") + + + +def geometric_correction(product, polarisation: str, to_print: bool = True): + """ + Performs geometric correction on the product. + + Parameters: + product: The input product. + polarisation (str): The polarization to process. + to_print (bool): Whether to print band names. + + Returns: + Product: The geometrically corrected product. + """ + parameters = HashMap() + parameters.put('demName', 'SRTM 1Sec HGT')# 'SRTM 3Sec') # SRTM 1Sec HGT + parameters.put('externalDEMNoDataValue', 0.0) + parameters.put('demResamplingMethod', "BILINEAR_INTERPOLATION") + parameters.put('imgResamplingMethod', "BILINEAR_INTERPOLATION") + parameters.put("Mask out areas with no elevation", True) + parameters.put('sourceBands', list(speckle.getBandNames())[0]) + parameters.put('nodataValueAtSea', False) + parameters.put('pixelSpacingInMeter', 10.0) + parameters.put('mapProjection',"WGS84(DD)") + terrain = GPF.createProduct('Terrain-Correction', parameters, speckle) + if to_print: + print("Bands: %s" % (list(terrain.getBandNames()))) + return terrain + +def subset(product, bbox: list, to_print: bool = True): + """ + Subsets the product using a bounding box. + + Parameters: + product: The input product. + bbox (list): Bounding box coordinates [min_lon, min_lat, max_lon, max_lat]. + to_print (bool): Whether to print band names. + + Returns: + Product or None: The subsetted product or None if no bands. + """ + parameters = HashMap() + min_lon, min_lat, max_lon, max_lat = bbox[0], bbox[1], bbox[2], bbox[3] + wkt_polygon = f"POLYGON(({min_lon} {min_lat}, {max_lon} {min_lat}, {max_lon} {max_lat}, {min_lon} {max_lat}, {min_lon} {min_lat}))" + geom = WKTReader().read(wkt_polygon) + parameters.put('copyMetadata', True) + parameters.put('geoRegion', geom) + subset = GPF.createProduct('Subset', parameters, product) + + if to_print: + print("Bands: %s" % (list(subset.getBandNames()))) + if len(list(subset.getBandNames())) == 0: + return None + else: + return subset + + +def subset_with_bbox(source_product, min_lon: float, min_lat: float, max_lon: float, max_lat: float): + """ + Subsets the product using specified bounding box coordinates. + + Parameters: + source_product: The input product. + min_lon (float): Minimum longitude. + min_lat (float): Minimum latitude. + max_lon (float): Maximum longitude. + max_lat (float): Maximum latitude. + + Returns: + Product: The subsetted product. + """ + parameters = HashMap() + wkt = f"POLYGON(({min_lon} {min_lat}, {max_lon} {min_lat}, {max_lon} {max_lat}, {min_lon} {max_lat}, {min_lon} {min_lat}))" + geom = WKTReader().read(wkt) + parameters.put('copyMetadata', True) + parameters.put('geoRegion', geom) + subset_product = GPF.createProduct('Subset', parameters, source_product) + return subset_product + + +def insar_pipeline(filename_1: str, filename_2: str, outpath: str, bbox: list, selectedPol: list): + """ + Main InSAR processing pipeline for Sentinel-1 data. + + Parameters: + master (str): Path to the master image. + slave (str): Path to the slave image. + outpath (str): Output directory path. + event_id (str): Event identifier. + bbox (list): Bounding box coordinates [min_lon, min_lat, max_lon, max_lat]. + other_param: Additional parameter for processing. + """ + for pol in selectedPol: + product_1 = read(filename_1) + product_2 = read(filename_2) + + swaths = ['IW1', 'IW2', 'IW3'] + ind = 0 + final_products = [] + failedPreprocessing = False + for swath in swaths: + try: + print('TOPSAR-Split') + product_TOPSAR_1 = topsar_split(product_1, swath, pol) + product_TOPSAR_2 = topsar_split(product_2, swath, pol) + except Exception as e: + print(e) + failedPreprocessing = True + break + product_orbitFile_1 = apply_orbit_file(product_TOPSAR_1,2) + product_orbitFile_2 = apply_orbit_file(product_TOPSAR_2,2) + + backGeocoding = back_geocoding([product_orbitFile_1, product_orbitFile_2]) + coherence = coherence_generation(backGeocoding) + TOPSAR_deburst = topsar_deburst(coherence, pol) + final_products.append(TOPSAR_deburst) + + if failedPreprocessing: + return + + merged = topsar_merge(final_products, pol) + terrain = geometric_correction(merged, pol) + band_names = list(terrain.getBandNames()) + + extension = f"{band_names[0]}_{str(pol)}_buffer0_window10_clipped.tif" + + output = os.path.join(outpath, extension) + + + subset_data = subset(terrain,bbox) + tmp = os.path.join(output) + parameters = HashMap() + parameters.put('file', tmp) + parameters.put('formatName', 'GeoTIFF') + + if not os.path.exists(output): + WriteOp = jpy.get_type('org.esa.snap.core.gpf.common.WriteOp') + w = WriteOp(subset_data, File(tmp), 'GeoTIFF') + w.writeProduct(ProgressMonitor.NULL) + + subset_data.dispose() + terrain.dispose() + product_1.dispose() + product_2.dispose() + product_TOPSAR_1.dispose() + product_TOPSAR_2.dispose() + product_orbitFile_2.dispose() + product_orbitFile_1.dispose() + backGeocoding.dispose() + coherence.dispose() + TOPSAR_deburst.dispose() + merged.dispose() + gc.collect() + + +def main(): + """ + Main function to process Sentinel-1 data using the SNAP toolbox via snappy. + """ + + parser = argparse.ArgumentParser(description='Process Sentinel-1 data using SNAP.') + parser.add_argument('master', help='Path to the master image') + parser.add_argument('slave', help='Path to the slave image') + parser.add_argument('outpath', help='Output path') + parser.add_argument('bbox', nargs=4, type=float, + help='Bounding box coordinates: [minLon, minLat, maxLon, maxLat]') + + # TODO add the polarizations as parameter + pols = ['VH','VV'] + + args = parser.parse_args() + + # Extract arguments + master = args.master + slave = args.slave + outpath = args.outpath + bbox = args.bbox + + insar_pipeline(master, slave, outpath,bbox,pols) + + +if __name__ == "__main__": + main() + + +