diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 8af99a25..a0184d63 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -4,6 +4,8 @@ on: pull_request: branches: - main + - 'v?*b' + - 'v?*beta' push: branches: - main @@ -15,40 +17,13 @@ permissions: contents: write env: - SKIP_DEPLOY: false - SKIP_INSTALL: true + SKIP_DEPLOY: true + SKIP_INSTALL: false SKIP_TEST: false jobs: - before_script: - runs-on: self-hosted - - strategy: - matrix: - node-version: [ 20.x ] - - steps: - - name: Checkout code - uses: actions/checkout@v4 - - - name: Set up Python - uses: actions/setup-python@v5 - with: - python-version: '3.x' - - - name: Install dependencies - run: | - source /opt/conda/etc/profile.d/conda.sh - conda init bash - source ~/.bashrc - conda env remove --name ci_temp_env --yes || echo "Environment ci_temp_env does not exist" - conda create --prefix /home/runneruser/.conda/envs/ci_temp_env --clone ci_env - conda activate ci_temp_env - shell: bash - test_sarvey: runs-on: self-hosted - needs: before_script steps: - name: Checkout code uses: actions/checkout@v4 @@ -64,11 +39,12 @@ jobs: source /opt/conda/etc/profile.d/conda.sh conda init bash source ~/.bashrc - conda activate ci_temp_env + conda activate ci_env rm -rf tests/testdata # wget -nv -c -O testdata.zip https://seafile.projekt.uni-hannover.de/f/4b3be399dffa488e98db/?dl=1 wget -nv -c -O testdata.zip https://seafile.projekt.uni-hannover.de/f/104b499f6f7e4360877d/?dl=1 unzip testdata.zip + mkdir -p test mv testdata tests/ mamba list make pytest @@ -79,7 +55,7 @@ jobs: run: | conda init bash source ~/.bashrc - source activate ci_temp_env + source activate ci_env # replace documentation address for tags befro make docs IFS='/' read -r OWNER REPO <<< "${GITHUB_REPOSITORY}" @@ -165,7 +141,7 @@ jobs: source /opt/conda/etc/profile.d/conda.sh conda init bash source ~/.bashrc - conda activate ci_temp_env + conda activate ci_env make lint shell: bash @@ -208,7 +184,7 @@ jobs: source /opt/conda/etc/profile.d/conda.sh conda init bash source ~/.bashrc - source activate ci_temp_env + source activate ci_env make urlcheck shell: bash @@ -228,92 +204,112 @@ jobs: run: | source /opt/conda/etc/profile.d/conda.sh conda activate base - mamba env remove --name sarvey_testinstall --yes || echo "Environment sarvey_testinstall does not exist" - mamba clean --index-cache --tarballs --packages -y - pip install conda-merge - export PATH=$HOME/.local/bin:$PATH # workaround conda-merge not found! - wget https://raw.githubusercontent.com/insarlab/MiaplPy/main/conda-env.yml - conda-merge conda-env.yml tests/CI_docker/context/environment_sarvey.yml > env.yml - mamba env create --name sarvey_testinstall -f env.yml + + conda create -n sarvey_testinstall python=3.10 pip -y conda activate sarvey_testinstall + conda install conda-forge::pysolid -y + conda install conda-forge::gdal pip install git+https://github.com/insarlab/MiaplPy.git - pip install . + + # get current branch for installation + IFS='/' read -r OWNER REPO <<< "${GITHUB_REPOSITORY}" + + URL="https://github.com/${OWNER}/${REPO}" + URL_IO="https://${OWNER}.github.io/${REPO}" + + echo "Repository URL: ${URL}" + echo "Repository Documentation URL: ${URL}" + + if [[ "$GITHUB_REF" == refs/tags/* ]]; then + current_ref="${GITHUB_REF##*/}" + elif [[ "$GITHUB_REF" == refs/heads/* ]]; then + current_ref="${GITHUB_REF#refs/heads/}" + elif [[ "$GITHUB_REF" == refs/pull/* ]]; then + current_ref="${GITHUB_HEAD_REF}" + else + current_ref="${GITHUB_REF}" + fi + + echo "Current branch/tag: ${URL}.git@$current_ref" + pip install git+${URL}.git@$current_ref + pip install sarvey[dev] + OUTPUT=$(pip check) || { echo "$OUTPUT"; true; } - cd .. + conda list python -c "import sarvey; print(sarvey)" shell: bash - deploy_pages: - runs-on: self-hosted - - needs: - - test_sarvey - - test_urls - - test_styles - if: github.ref == 'refs/heads/main' || startsWith(github.ref, 'refs/tags/') - steps: - - name: Checkout code - uses: actions/checkout@v4 - - - name: Set up Python - uses: actions/setup-python@v5 - with: - python-version: '3.x' - - - name: Download docs - uses: actions/download-artifact@v4 - with: - name: docs - path: docs/_build/html/ - - - name: Download coverage report - uses: actions/download-artifact@v4 - with: - name: coverage-report - path: htmlcov/ - - - name: Download report.html - uses: actions/download-artifact@v4 - with: - name: test-report - - - name: Deploy to GitHub Pages - # trigger if merged into the main branch || published new tag - if: env.SKIP_DEPLOY == 'false' && github.event_name != 'pull_request' && (github.ref == 'refs/heads/main' || startsWith(github.ref, 'refs/tags/')) - run: | - rm -rf public + # deploy_pages: + # runs-on: self-hosted + + # needs: + # - test_sarvey + # - test_urls + # - test_styles + # if: github.ref == 'refs/heads/main' || startsWith(github.ref, 'refs/tags/') + # steps: + # - name: Checkout code + # uses: actions/checkout@v4 + + # - name: Set up Python + # uses: actions/setup-python@v5 + # with: + # python-version: '3.x' + + # - name: Download docs + # uses: actions/download-artifact@v4 + # with: + # name: docs + # path: docs/_build/html/ + + # - name: Download coverage report + # uses: actions/download-artifact@v4 + # with: + # name: coverage-report + # path: htmlcov/ + + # - name: Download report.html + # uses: actions/download-artifact@v4 + # with: + # name: test-report + + # - name: Deploy to GitHub Pages + # # trigger if merged into the main branch || published new tag + # if: env.SKIP_DEPLOY == 'false' && github.event_name != 'pull_request' && (github.ref == 'refs/heads/main' || startsWith(github.ref, 'refs/tags/')) + # run: | + # rm -rf public - git clone --branch gh-pages https://github.com/${{ github.repository }} public - - if [[ "${GITHUB_REF}" == refs/tags/* ]]; then - TAG_NAME=${GITHUB_REF##*/} - echo "Deploying to GitHub Pages for version tag: $TAG_NAME" - DOCS_PATH=public/tags/$TAG_NAME - else - echo "Deploying to GitHub Pages for main branch" - DOCS_PATH=public/main - fi - - rm -rf $DOCS_PATH - mkdir -p $DOCS_PATH/docs - mkdir -p $DOCS_PATH/images - mkdir -p $DOCS_PATH/coverage - mkdir -p $DOCS_PATH/test_reports + # git clone --branch gh-pages https://github.com/${{ github.repository }} public + + # if [[ "${GITHUB_REF}" == refs/tags/* ]]; then + # TAG_NAME=${GITHUB_REF##*/} + # echo "Deploying to GitHub Pages for version tag: $TAG_NAME" + # DOCS_PATH=public/tags/$TAG_NAME + # else + # echo "Deploying to GitHub Pages for main branch" + # DOCS_PATH=public/main + # fi + + # rm -rf $DOCS_PATH + # mkdir -p $DOCS_PATH/docs + # mkdir -p $DOCS_PATH/images + # mkdir -p $DOCS_PATH/coverage + # mkdir -p $DOCS_PATH/test_reports - cp -r docs/_build/html/* $DOCS_PATH/docs - cp -r htmlcov/* $DOCS_PATH/coverage/ - cp report.html $DOCS_PATH/test_reports/ - - ls -al $DOCS_PATH - ls -al $DOCS_PATH/docs - ls -al $DOCS_PATH/coverage - ls -al $DOCS_PATH/test_reports - - shell: bash - - - name: Upload to GitHub Pages - if: env.SKIP_DEPLOY == 'false' - uses: peaceiris/actions-gh-pages@v4 - with: - github_token: ${{ secrets.GITHUB_TOKEN }} - publish_dir: ./public + # cp -r docs/_build/html/* $DOCS_PATH/docs + # cp -r htmlcov/* $DOCS_PATH/coverage/ + # cp report.html $DOCS_PATH/test_reports/ + + # ls -al $DOCS_PATH + # ls -al $DOCS_PATH/docs + # ls -al $DOCS_PATH/coverage + # ls -al $DOCS_PATH/test_reports + + # shell: bash + + # - name: Upload to GitHub Pages + # if: env.SKIP_DEPLOY == 'false' + # uses: peaceiris/actions-gh-pages@v4 + # with: + # github_token: ${{ secrets.GITHUB_TOKEN }} + # publish_dir: ./public diff --git a/HISTORY.rst b/HISTORY.rst index adebb52d..c5a64fa5 100644 --- a/HISTORY.rst +++ b/HISTORY.rst @@ -2,10 +2,21 @@ History ======= +Future major version +-------------------- + +* Export data to GIS format in WGS84 coordinates. +* Change file name coordinates_UTM.h5 to coordinates_map.h5. +* Use Transverse Mercator instead of UTM as map coordinates. + Future minor version (release soon) ----------------------------------- +* Update CI docker builder +* Update runner to test installation +* Update documentation with new instruction for installation including pip +* Fix numerical problems when computing grid size. 1.2.0 (2025-02-19) ------------------ diff --git a/README.rst b/README.rst index 6d14f21a..f0c2243f 100644 --- a/README.rst +++ b/README.rst @@ -10,7 +10,7 @@ Open-source InSAR time series analysis software developed within the project SAR Documentation ------------- The documentation with installation instructions, processing steps, and examples with a demo dataset can be found at: -https://luhipi.github.io/sarvey/main/docs/ +https://luhipi.github.io/sarvey/main Discussion ---------- @@ -26,16 +26,12 @@ Status :target: https://github.com/luhipi/sarvey/actions :alt: Pipelines .. image:: https://img.shields.io/static/v1?label=Documentation&message=GitHub%20Pages&color=blue - :target: https://luhipi.github.io/sarvey/main/docs/ + :target: https://luhipi.github.io/sarvey/main :alt: Documentation .. image:: https://zenodo.org/badge/DOI/10.5281/zenodo.12544130.svg :target: https://doi.org/10.5281/zenodo.12544130 :alt: DOI - -See also the latest coverage_ report and the pytest_ HTML report. - - License ------- @@ -191,15 +187,13 @@ This package was created with Cookiecutter_ and the `fernlab/cookiecutter-python .. _Cookiecutter: https://github.com/audreyr/cookiecutter .. _`fernlab/cookiecutter-python-package`: https://git.gfz-potsdam.de/fernlab/products/cookiecutters/cookiecutter-python-package -.. _coverage: https://luhipi.github.io/sarvey/main/coverage/ -.. _pytest: https://luhipi.github.io/sarvey/main/test_reports/report.html -.. _processing: https://luhipi.github.io/sarvey/main/docs/processing.html -.. _`processing steps`: https://luhipi.github.io/sarvey/main/docs/processing.html#processing-steps-for-two-step-unwrapping-workflow -.. _preparation: https://luhipi.github.io/sarvey/main/docs/preparation.html -.. _`Two-step unwrapping`: https://luhipi.github.io/sarvey/main/docs/processing.html#processing-steps-for-two-step-unwrapping-workflow -.. _`One-step unwrapping`: https://luhipi.github.io/sarvey/main/docs/processing.html#processing-steps-for-one-step-unwrapping-workflow -.. _`installation instruction`: https://luhipi.github.io/sarvey/main/docs/installation.html -.. _`history`: https://luhipi.github.io/sarvey/main/docs/history.html +.. _processing: https://luhipi.github.io/sarvey/main/processing.html +.. _`processing steps`: https://luhipi.github.io/sarvey/main/processing.html#processing-steps-for-two-step-unwrapping-workflow +.. _preparation: https://luhipi.github.io/sarvey/main/preparation.html +.. _`Two-step unwrapping`: https://luhipi.github.io/sarvey/main/processing.html#processing-steps-for-two-step-unwrapping-workflow +.. _`One-step unwrapping`: https://luhipi.github.io/sarvey/main/processing.html#processing-steps-for-one-step-unwrapping-workflow +.. _`installation instruction`: https://luhipi.github.io/sarvey/main/installation.html +.. _`history`: https://luhipi.github.io/sarvey/main/history.html .. _MiaplPy: https://github.com/insarlab/MiaplPy .. _MintPy: https://github.com/insarlab/MintPy .. _ISCE: https://github.com/isce-framework/isce2 diff --git a/docs/demo/demo_masjed_dam_detailed_guide.rst b/docs/demo/demo_masjed_dam_detailed_guide.rst index bb4df9df..d88cda51 100644 --- a/docs/demo/demo_masjed_dam_detailed_guide.rst +++ b/docs/demo/demo_masjed_dam_detailed_guide.rst @@ -174,7 +174,7 @@ After running this step, a `sbas` directory is created. Inside this directory, y ├── temporal_coherence.h5 ├── ifg_stack.h5 ├── ifg_network.h5 - ├── coordinates_utm.h5 + ├── coordinates_map.h5 ├── config.json ├── background_map.h5 └── pic/ @@ -261,15 +261,12 @@ Outputs of this step are: .. code-block:: none outputs/ - ├── p2_coh80_ifg_wr.h5 - ├── p2_coh80_aps.h5 ├── p1_aps.h5 ├── p1_ts_filt.h5 + ├── aps_parameters.h5 └── pic/ ├── step_3_temporal_autocorrelation.png - ├── step_3_stable_points.png - ├── selected_pixels_temp_coh_0.8.png - └── step_3_mask_p2_coh80.png + └── step_3_stable_points.png Step 2.4: Run Step 4 of SARvey @@ -279,7 +276,18 @@ Step 2.4: Run Step 4 of SARvey sarvey -f config.json 4 4 -.. outputs directory structure to be added +Outputs of this step are: + +.. code-block:: none + + outputs/ + ├── coh80_ifg_wr.h5 + ├── coh80_aps.h5 + ├── coh80_ifg_unw.h5 + ├── coh80_ts.h5 + └── pic/ + ├── selected_pixels_temp_coh_0.8.png + └── step_4_mask_coh80.png The results of step 4 of SARvey, including the time series, are stored in the `p2_coh80_ts.h5` file. The file is named based on the `coherence_p2` parameter in the config.json file. @@ -310,11 +318,11 @@ Step 4: Modify Config File and Rerun SARvey Modify the config.json file and change **coherence_p2** from 0.8 to 0.7. -Run steps 3 and 4 using the following command: +Rerun steps 4 using the following command: .. code-block:: bash - sarvey -f config.json 3 4 + sarvey -f config.json 4 4 A new file `p2_coh70_ts.h5` is created. You can now visualize this file that has a higher point density. diff --git a/docs/demo/demo_masjed_dam_fast_track.rst b/docs/demo/demo_masjed_dam_fast_track.rst index bededb6c..360f11aa 100644 --- a/docs/demo/demo_masjed_dam_fast_track.rst +++ b/docs/demo/demo_masjed_dam_fast_track.rst @@ -51,7 +51,7 @@ You can run each step individually or a range of steps by specifying the first a Check Outputs """"""""""""" -First, check the output snapshots in the `outputs/pics` directory. You can also use **`sarvey_plot`** to plot various products to assess the quality of the results and decide how to adjust parameters. Modify the parameters in the config file and rerun the corresponding steps of `sarvey` to improve the results. For instance, changing **`coherence_p2`** from 0.8 to 0.7 and rerunning steps 3 and 4 can increase the density of the final set of points. However, be cautious that reducing the value too much may include noisy points of low quality in the analysis, potentially leading to poor final results. You can check the details of all parameters using the -p flag in `sarvey` and decide how to tune them. For more explanations, please refer to :ref:`processing` +First, check the output snapshots in the `outputs/pics` directory. You can also use **`sarvey_plot`** to plot various products to assess the quality of the results and decide how to adjust parameters. Modify the parameters in the config file and rerun the corresponding steps of `sarvey` to improve the results. For instance, changing **`coherence_p2`** from 0.8 to 0.7 and rerunning step 4 can increase the density of the final set of points. However, be cautious that reducing the value too much may include noisy points of low quality in the analysis, potentially leading to poor final results. You can check the details of all parameters using the -p flag in `sarvey` and decide how to tune them. For more explanations, please refer to :ref:`processing` diff --git a/docs/installation.rst b/docs/installation.rst index fb2d49e3..9e50d8f2 100644 --- a/docs/installation.rst +++ b/docs/installation.rst @@ -13,7 +13,7 @@ SARvey is a cross-platform python-based software and can be installed on Linux ----- -On Linux, SARvey can be installed `Using Mamba (recommended)`_ or `Using Anaconda or Miniconda`_. +On Linux, SARvey can be installed `Using Mamba (recommended)`_ or `Using Anaconda or Miniconda`_ or `Using Pip`_. Using Mamba (recommended) ^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -21,7 +21,7 @@ Using Mamba (recommended) Using mamba_ (latest version recommended), **SARvey** is installed as follows: -1. Clone the SARvey source code and install SARvey and all dependencies from the environment_sarvey.yml file: +1. Clone the SARvey source code and install SARvey and all dependencies from the environment.yml file: .. code-block:: bash @@ -35,7 +35,7 @@ Using mamba_ (latest version recommended), **SARvey** is installed as follows: pip install conda-merge wget https://raw.githubusercontent.com/insarlab/MiaplPy/main/conda-env.yml - conda-merge conda-env.yml tests/CI_docker/context/environment_sarvey.yml > env.yml + conda-merge conda-env.yml environment.yml > env.yml mamba env create -n sarvey -f env.yml rm env.yml conda-env.yml mamba activate sarvey @@ -43,17 +43,13 @@ Using mamba_ (latest version recommended), **SARvey** is installed as follows: pip install . -This is the preferred method to install **SARvey**, as it always installs the most recent stable release and -automatically resolves all the dependencies. - - Using Anaconda or Miniconda ^^^^^^^^^^^^^^^^^^^^^^^^^^^ Using conda_ (latest version recommended), **SARvey** is installed as follows: -1. Then clone the **SARvey** source code and install **SARvey** and all dependencies from the environment_sarvey.yml file: +1. Then clone the **SARvey** source code and install **SARvey** and all dependencies from the environment.yml file: .. code-block:: bash @@ -67,7 +63,7 @@ Using conda_ (latest version recommended), **SARvey** is installed as follows: pip install conda-merge wget https://raw.githubusercontent.com/insarlab/MiaplPy/main/conda-env.yml - conda-merge conda-env.yml tests/CI_docker/context/environment_sarvey.yml > env.yml + conda-merge conda-env.yml environment.yml > env.yml conda env create -n sarvey -f env.yml rm env.yml conda-env.yml conda activate sarvey @@ -75,6 +71,45 @@ Using conda_ (latest version recommended), **SARvey** is installed as follows: pip install . +Using pip +^^^^^^^^^ + +Using pip_ (latest version recommended), **SARvey** is installed as follows: + +1. Create a new environment for **SARvey** (optional but recommended): + + .. code-block:: bash + + conda create -n sarvey pip -y + conda activate sarvey + +2. Install dependencies + + .. code-block:: bash + + conda install -c conda-forge pysolid gdal + + +3. Install **Miaplpy** + + .. code-block:: bash + + pip install git+https://github.com/insarlab/MiaplPy.git + +4. Install **SARvey** + + .. code-block:: bash + + pip install git+https://github.com/luhipi/sarvey.git + + +If your are a developer, install the development requirements using the following command. + + .. code-block:: bash + + pip install sarvey[dev] + + MacOS ARM (Apple Silicon M2) ---------------------------- @@ -120,12 +155,12 @@ Using conda_ (latest version recommended), SARvey is installed as follows: git clone https://github.com/luhipi/sarvey.git cd sarvey - 2.2 Open `tests/CI_docker/context/environment_sarvey.yml` in an editor of your choice and comment out the lines `isce2` and `gcc_linux-64`. Alternatively, you can run the following commands. + 2.2 Open `environment.yml` in an editor of your choice and comment out the lines `isce2` and `gcc_linux-64`. Alternatively, you can run the following commands. .. code-block:: bash - sed -i '' '/isce2/s/^/# /' tests/CI_docker/context/environment_sarvey.yml - sed -i '' '/gcc_linux-64/s/^/# /' tests/CI_docker/context/environment_sarvey.yml + sed -i '' '/isce2/s/^/# /' environment.yml + sed -i '' '/gcc_linux-64/s/^/# /' environment.yml Note: As of the time of creation of this document, `isce2` for MacOS ARM64 is not available in Conda repositories. Therefore, it is skipped, but it should not cause any problems for running SARvey. Also, `gcc_linux-64` is not required on ARM64. @@ -133,7 +168,7 @@ Using conda_ (latest version recommended), SARvey is installed as follows: .. code-block:: bash - conda env update --name sarvey -f tests/CI_docker/context/environment_sarvey.yml + conda env update --name sarvey -f environment.yml conda activate sarvey pip install . @@ -183,4 +218,5 @@ On Windows, SARvey is tested on Windows Subsystem for Linux (WSL_) version 2. Pl .. _mamba: https://github.com/mamba-org/mamba .. _Conda for Mac: https://docs.conda.io/projects/conda/en/latest/user-guide/install/macos.html .. _WSL: https://learn.microsoft.com/en-us/windows/wsl/ +.. _MiaplPy: https://github.com/insarlab/MiaplPy diff --git a/docs/processing.rst b/docs/processing.rst index f14e1b3f..59a0ac82 100644 --- a/docs/processing.rst +++ b/docs/processing.rst @@ -109,12 +109,12 @@ Step 0: Preparation - Estimating the temporal coherence: The phase noise of each pixel is approximated by the estimation of the temporal phase coherence (Zhao and Mallorqui 2019). Thereby, a low-pass filter with a certain window size is used (**preparation:filter_window_size**). - The temporal coherence is used to select the first- and second-order points in the later steps (**consistency_check:coherence_p1** and **filtering:coherence_p2**). + The temporal coherence is used to select the first- and second-order points in the later steps (**consistency_check:coherence_p1** and **densification:coherence_p2**). - Output of this step - background_map.h5 - ifg_stack.h5 - - coordinates_utm.h5 + - coordinates_map.h5 - ifg_network.h5 - temporal_coherence.h5 @@ -223,17 +223,6 @@ However, the step 3 has to be executed as the second-order points are selected d Points with a non-linear displacement behaviour are removed by a threshold on the temporal autocorrelation of the displacement time series (**filtering:max_temporal_autocorrelation**) (Crosetto et al. 2018). A regular grid (**filtering:grid_size** in [m]) is applied to select the first-order points with the lowest temporal autocorrelation to reduce the computational complexity during filtering. -- Selecting second-order points: - Second-order points are selected based on a temporal coherence threshold (**filtering:coherence_p2**) on the temporal phase coherence computed during step 0. - A mask file can be specified (**filtering:mask_p2_file**) to limit the second-order points to the given area of interest. - Second-order points can also be selected based on the results of phase-linking (set **phase_linking:use_phase_linking_results** to True) implemented in MiaplPy (Mirzaee et al. 2023). - More information on Miaplpy and phase-linking can be found `here `_. - The number of siblings (**phase_linking:num_siblings**) used during phase-linking within MiaplPy processing needs to be specified to identify the distributed scatterers (DS) among the pixels selected by MiaplPy. - A mask file can be specified (**phase_linking:mask_phase_linking_file**) to limit the phase-linking to the given area of interest. - MiaplPy also provides a selection of persistent scatterers (PS) which can be included as second-order points (set **phase_linking:use_ps** to True) and also specify the path to the maskPS.h5 (**phase_linking:mask_ps_file**) which is also an output of MiaplPy. - In case the second-order points are selected among the results from MiaplPy, the filtered interferometric phase (MiaplPy result) is used for the respective points. - The DS pixels from MiaplPy and the pixels selected with the temporal phase coherence from step 0 are both selected with the same coherence threshold (**filtering:coherence_p2**). - - Estimating the atmospheric phase screen (APS): The estimation of the APS takes place in time-domain and not interferogram-domain to reduce the computational time. The phase contributions are removed from the first-order points which were selected for atmospheric filtering. @@ -244,27 +233,41 @@ However, the step 3 has to be executed as the second-order points are selected d - Output of this step - p1_ts_filt.h5 - p1_aps.h5 - - p2_cohXX_aps.h5 - - p2_cohXX_ifg_wr.h5 + - aps_parameters.h5 -The placeholder XX depends on the threshold for the temporal coherence used for selecting the second-order points. -For example, a threshold of 0.8 would result in p2_coh80_aps.h5 and p2_coh80_ifg_wr.h5. Step 4: Densification ^^^^^^^^^^^^^^^^^^^^^ The densification step is the last step of the two-step unwrapping workflow. So far, the displacement was only estimated at the sparse locations of the first-order points. -The second-order points selected during step 3 (filtering) are added to the first-order points to densify the final set of points. -During the densification, first, the estimated APS is removed from both first- and second-order points and second, the displacement time series are retrieved by unwrapping phases of the points jointly. +Now, second order points are selected and added to the first-order points to densify the final set of points. +Afterwards, the APS, estimated in step 3, is interpolated to the location of the second-order points and removed from first- and second-order points. +Lastly, the displacement time series are retrieved by unwrapping phases of the points jointly. Two unwrapping options (**general:apply_temporal_unwrapping**, also applies to step 2) are implemented and should be chosen based on the characteristics of the displacement (spatial extend, magnitude, temporal behaviour). +- Selecting second-order points: + Second-order points are selected based on a temporal coherence threshold (**densification:coherence_p2**) on the temporal phase coherence computed during step 0. + A mask file can be specified (**densification:mask_p2_file**) to limit the second-order points to the given area of interest. + Second-order points can also be selected based on the results of phase-linking (set **phase_linking:use_phase_linking_results** to True) implemented in MiaplPy (Mirzaee et al. 2023). + More information on Miaplpy and phase-linking can be found `here `_. + The number of siblings (**phase_linking:num_siblings**) used during phase-linking within MiaplPy processing needs to be specified to identify the distributed scatterers (DS) among the pixels selected by MiaplPy. + A mask file can be specified (**phase_linking:mask_phase_linking_file**) to limit the phase-linking to the given area of interest. + MiaplPy also provides a selection of persistent scatterers (PS) which can be included as second-order points (set **phase_linking:use_ps** to True) and also specify the path to the maskPS.h5 (**phase_linking:mask_ps_file**) which is also an output of MiaplPy. + In case the second-order points are selected among the results from MiaplPy, the filtered interferometric phase (MiaplPy result) is used for the respective points. + The DS pixels from MiaplPy and the pixels selected with the temporal phase coherence from step 0 are both selected with the same coherence threshold (**filtering:coherence_p2**). + +- Apply estimated atmospheric phase screen (APS) to second-order points: + The APS, estimated in step 3, is interpolated to the location of the second-order points. + - Output of this step + - p2_cohXX_ifg_wr.h5 + - p2_cohXX_aps.h5 - p2_cohXX_ifg_unw.h5 - p2_cohXX_ts.h5 -The placeholder XX depends on the threshold for the temporal coherence used for selecting the second-order points during filtering in step 3. +The placeholder XX depends on the threshold for the temporal coherence used for selecting the second-order points. For example, a threshold of 0.8 would result in p2_coh80_ifg_unw.h5 and p2_coh80_ts.h5. Option 1: Unwrapping in time and space diff --git a/tests/CI_docker/context/environment_sarvey.yml b/environment.yml similarity index 100% rename from tests/CI_docker/context/environment_sarvey.yml rename to environment.yml diff --git a/sarvey/config.py b/sarvey/config.py index e746db84..9405989d 100644 --- a/sarvey/config.py +++ b/sarvey/config.py @@ -467,12 +467,6 @@ class Unwrapping(BaseModel, extra=Extra.forbid): class Filtering(BaseModel, extra=Extra.forbid): """Template for filtering settings in config file.""" - coherence_p2: float = Field( - title="Temporal coherence threshold", - description="Set the temporal coherence threshold for the filtering step.", - default=0.8 - ) - apply_aps_filtering: bool = Field( title="Apply atmosphere filtering", description="Set whether to filter atmosphere or to skip it.", @@ -491,12 +485,6 @@ class Filtering(BaseModel, extra=Extra.forbid): default=1000 ) - mask_p2_file: Optional[str] = Field( - title="Path to spatial mask file for second-order points.", - description="Path to the mask file, e.g. created by sarvey_mask.", - default="" - ) - use_moving_points: bool = Field( title="Use moving points", description="Set whether to use moving points in the filtering step.", @@ -509,15 +497,6 @@ class Filtering(BaseModel, extra=Extra.forbid): default=0.3 ) - @validator('coherence_p2') - def checkTempCohThrsh2(cls, v): - """Check if the temporal coherence threshold is valid.""" - if v < 0: - raise ValueError("Temporal coherence threshold cannot be negative.") - if v > 1: - raise ValueError("Temporal coherence threshold cannot be greater than 1.") - return v - @validator('interpolation_method') def checkInterpolationMethod(cls, v): """Check if the interpolation method is valid.""" @@ -534,16 +513,6 @@ def checkGridSize(cls, v): else: return v - @validator('mask_p2_file') - def checkSpatialMaskPath(cls, v): - """Check if the path is correct.""" - if v == "" or v is None: - return None - else: - if not os.path.exists(os.path.abspath(v)): - raise ValueError(f"mask_p2_file path is invalid: {v}") - return v - @validator('max_temporal_autocorrelation') def checkMaxAutoCorr(cls, v): """Check if the value is correct.""" @@ -555,6 +524,24 @@ def checkMaxAutoCorr(cls, v): class Densification(BaseModel, extra=Extra.forbid): """Template for densification settings in config file.""" + coherence_p2: float = Field( + title="Temporal coherence threshold", + description="Set the temporal coherence threshold for the densification step.", + default=0.8 + ) + + coherence_threshold: float = Field( + title="Coherence threshold for densification", + description="Set coherence threshold for densification.", + default=0.5 + ) + + mask_p2_file: Optional[str] = Field( + title="Path to spatial mask file for second-order points.", + description="Path to the mask file, e.g. created by sarvey_mask.", + default="" + ) + num_connections_to_p1: int = Field( title="Number of connections in temporal unwrapping.", description="Set number of connections between second-order point and closest first-order points for temporal " @@ -593,11 +580,30 @@ class Densification(BaseModel, extra=Extra.forbid): default=0.5 ) + @validator('coherence_p2') + def checkTempCohThrsh2(cls, v): + """Check if the temporal coherence threshold is valid.""" + if v < 0: + raise ValueError("Temporal coherence threshold cannot be negative.") + if v > 1: + raise ValueError("Temporal coherence threshold cannot be greater than 1.") + return v + + @validator('mask_p2_file') + def checkSpatialMaskPath(cls, v): + """Check if the path is correct.""" + if v == "" or v is None: + return None + else: + if not os.path.exists(os.path.abspath(v)): + raise ValueError(f"mask_p2_file path is invalid: {v}") + return v + @validator('num_connections_to_p1') def checkNumConn1(cls, v): - """Check if num_connections_p1 are valid.""" + """Check if num_connections_to_p1 are valid.""" if v <= 0: - raise ValueError(f"num_connections_p1 must be greater than 0: {v}") + raise ValueError(f"num_connections_to_p1 must be greater than 0: {v}") return v @validator('max_distance_to_p1') diff --git a/sarvey/densification.py b/sarvey/densification.py index d22f05e8..cc64ecf3 100644 --- a/sarvey/densification.py +++ b/sarvey/densification.py @@ -28,17 +28,25 @@ # with this program. If not, see . """Densification module for SARvey.""" +from os.path import join import time import multiprocessing import numpy as np from scipy.spatial import KDTree from logging import Logger +import matplotlib.pyplot as plt +import cmcrameri as cmc + +from mintpy.utils import ptime, readfile +from miaplpy.objects.slcStack import slcStack -from mintpy.utils import ptime from sarvey.unwrapping import oneDimSearchTemporalCoherence -from sarvey.objects import Points +from sarvey.objects import Points, AmplitudeImage, BaseStack, ApsParameters import sarvey.utils as ut +from sarvey.preparation import selectPixels, createTimeMaskFromDates +from sarvey.filtering import applySimpleInterpolationToP2, applySpatialFilteringToP2 +from sarvey.config import Config def densificationInitializer(tree_p1: KDTree, point2_obj: Points, demod_phase1: np.ndarray): @@ -120,8 +128,8 @@ def launchDensifyNetworkConsistencyCheck(args: tuple): for idx in range(num_points): p2 = idx_range[idx] # nearest points in p1 - dist, nearest_p1 = global_tree_p1.query([global_point2_obj.coord_utm[p2, 0], - global_point2_obj.coord_utm[p2, 1]], k=num_conn_p1) + dist, nearest_p1 = global_tree_p1.query([global_point2_obj.coord_map[p2, 0], + global_point2_obj.coord_map[p2, 1]], k=num_conn_p1) mask = (dist < max_dist_p1) & (dist != 0) mask[:3] = True # ensure that always at least the three closest points are used nearest_p1 = nearest_p1[mask] @@ -197,7 +205,7 @@ def densifyNetwork(*, point1_obj: Points, vel_p1: np.ndarray, demerr_p1: np.ndar start_time = time.time() # find the closest points from first-order network - tree_p1 = KDTree(data=point1_obj.coord_utm) + tree_p1 = KDTree(data=point1_obj.coord_map) # remove parameters from wrapped phase pred_phase_demerr, pred_phase_vel = ut.predictPhase( @@ -272,3 +280,237 @@ def densifyNetwork(*, point1_obj: Points, vel_p1: np.ndarray, demerr_p1: np.ndar vel_p2 = vel_p2[sort_idx] gamma_p2 = gamma_p2[sort_idx] return demerr_p2, vel_p2, gamma_p2 + + +def selectP2(*, output_path: str, config: Config, logger: Logger): + """Select second order points and interpolate APS to them. + + Parameters + ---------- + output_path : str + path to the directory with processing files. + config : dict + dictionay containing parameters. + logger : Logger + Logger instance for logging messages. + + Returns + ------- + point2_obj : Points + Points object with second-order points + aps2_phase: np.ndarray + APS estimation for second-order points + """ + # this function is mainly adapted from sarvey.processing.runFiltering + # TODO: Directly pass input parameters instead of config dictionary + coherence_p2 = config.densification.coherence_p2 + logger.debug("Start Selecting 2nd order points.") + + coh_value = int(coherence_p2 * 100) + + p1_file = join(output_path, "p1_ts_filt.h5") + logger.debug(f"Create a mask based on points in file {p1_file}.") + point1_obj = Points(file_path=p1_file, logger=logger) + point1_obj.open(input_path=config.general.input_path) + p1_mask = point1_obj.createMask() # used later for selecting psCand2 when a spatial mask AOI is given. + logger.debug(f"Mask created with {p1_mask.sum()} selected points.") + + bmap_file = join(output_path, "background_map.h5") + logger.debug(f"Use file {bmap_file} to create a mask invalid points.") + bmap_obj = AmplitudeImage(file_path=bmap_file) + + point_id_img = np.arange(0, point1_obj.length * point1_obj.width).reshape( + (point1_obj.length, point1_obj.width)) + + logger.debug(f"Select second-order points based on temporal coherence {coherence_p2:.2f}.") + cand_mask2 = selectPixels( + path=output_path, selection_method="temp_coh", + thrsh=coherence_p2, + grid_size=None, bool_plot=True, + logger=logger + ) # first-order points are included in second-order points + logger.debug(f"{cand_mask2.sum()} second-order points selected based on temporal coherence {coherence_p2:.2f}.") + + if config.phase_linking.use_phase_linking_results: + # read PL results + pl_file = join(config.phase_linking.inverted_path, "phase_series.h5") + logger.debug(f"Read phase linking results from file {pl_file}.") + pl_coh = readfile.read(pl_file, + datasetName='temporalCoherence')[0] + pl_coh = pl_coh[1, :, :] + siblings = readfile.read(pl_file, + datasetName='shp')[0] + + if config.phase_linking.use_ps: + logger.debug(f"Read phase linking PS mask from file {config.phase_linking.mask_ps_file}.") + mask_ps = readfile.read(config.phase_linking.mask_ps_file, + datasetName='mask')[0].astype(np.bool_) + cand_mask_pl = (pl_coh > coherence_p2) | mask_ps + else: + cand_mask_pl = (pl_coh > coherence_p2) + # remove ps candidates, because the ps detection strategy in miaplpy seems to be biased. + cand_mask_pl[siblings <= config.phase_linking.num_siblings] = False + + if config.phase_linking.mask_phase_linking_file is not None: + path_mask_pl_aoi = join(config.phase_linking.mask_phase_linking_file) + logger.debug(f"Load mask for area of interest from file {path_mask_pl_aoi}.") + mask_pl_aoi = readfile.read(path_mask_pl_aoi, datasetName='mask')[0].astype(np.bool_) + + fig = plt.figure(figsize=(15, 5)) + ax = fig.add_subplot() + ax.imshow(mask_pl_aoi, cmap=cmc.cm.cmaps["grayC"], alpha=0.5, zorder=10, vmin=0, vmax=1) + bmap_obj.plot(ax=ax, logger=logger) + coord_xy = np.array(np.where(cand_mask_pl)).transpose() + val = np.ones_like(cand_mask_pl) + sc = ax.scatter(coord_xy[:, 1], coord_xy[:, 0], c=val[cand_mask_pl], s=0.5, + cmap=cmc.cm.cmaps["lajolla_r"], + vmin=1, vmax=2) # set min, max to ensure that points are yellow + cbar = plt.colorbar(sc, pad=0.03, shrink=0.5) + cbar.ax.set_visible(False) # make size of axis consistent with all others + plt.tight_layout() + plt.title("Mask for phase linking points") + fig.savefig(join(output_path, "pic", "step_4_mask_p2_coh{}_phase_linking.png".format(coh_value)), dpi=300) + plt.close(fig) + + # mask points after plotting, so that available coherent points are visible in figure + cand_mask_pl[~mask_pl_aoi] = False + + # combine phase linking coherence with TPC cand_mask2 + initial_num_cand_mask2 = cand_mask2.sum() + cand_mask2 = cand_mask2 | cand_mask_pl + logger.debug(f"{cand_mask2.sum()-initial_num_cand_mask2} additional pixels selected by phase linking.") + + mask_valid_area = ut.detectValidAreas(bmap_obj=bmap_obj, logger=logger) + + if config.densification.mask_p2_file is not None: + path_mask_aoi = join(config.densification.mask_p2_file) + logger.info(f"Load mask for area of interest from file {path_mask_aoi}.") + mask_aoi = readfile.read(path_mask_aoi, datasetName='mask')[0].astype(np.bool_) + mask_valid_area &= mask_aoi + # todo: add unstable points from p1 for densification + else: + logger.info(msg="No mask for area of interest given.") + + cand_mask2[p1_mask] = True # add all selected 1.order points to avoid spatial gaps in 2D unwrapping + # cand_mask2[cand_mask_sparse] = True # add only stable points from 1.order points + + cand_mask2 &= mask_valid_area + logger.info(f"Final number of selected second-order points: {cand_mask2.sum()}.") + + fig = plt.figure(figsize=(15, 5)) + ax = fig.add_subplot() + ax.imshow(mask_valid_area, cmap=cmc.cm.cmaps["grayC"], alpha=0.5, zorder=10, vmin=0, vmax=1) + bmap_obj.plot(ax=ax, logger=logger) + coord_xy = np.array(np.where(cand_mask2)).transpose() + val = np.ones_like(cand_mask2) + sc = ax.scatter(coord_xy[:, 1], coord_xy[:, 0], c=val[cand_mask2], s=0.5, cmap=cmc.cm.cmaps["lajolla_r"], + vmin=0, vmax=10) # set min, max to ensure that points are yellow + cbar = plt.colorbar(sc, pad=0.03, shrink=0.5) + cbar.ax.set_visible(False) # make size of axis consistent with all others + plt.tight_layout() + plt.title("Mask for dense point set") + fig.savefig(join(output_path, "pic", "step_4_mask_p2_coh{}.png".format(coh_value)), dpi=300) + plt.close(fig) + + p2_file = join(output_path, f"p2_coh{coh_value}_ifg_wr.h5") + logger.debug(f"Prepare second-order points object to save to file {p2_file}.") + point2_obj = Points(file_path=p2_file, logger=logger) + coord_xy = np.array(np.where(cand_mask2)).transpose() + point_id2 = point_id_img[cand_mask2] + point2_obj.prepare( + point_id=point_id2, + coord_xy=coord_xy, + input_path=config.general.input_path + ) + + ifg_stack_file = join(output_path, "ifg_stack.h5") + logger.debug(f"Read interferogram phase for selected second-order points from file {ifg_stack_file}.") + ifg_stack_obj = BaseStack(file=ifg_stack_file, logger=logger) + point2_obj.phase = ut.readPhasePatchwise(stack_obj=ifg_stack_obj, dataset_name="ifgs", + num_patches=config.general.num_patches, cand_mask=cand_mask2, + point_id_img=point_id_img, logger=logger) + + if config.phase_linking.use_phase_linking_results: + pl_inverted_file = join(config.phase_linking.inverted_path, "phase_series.h5") + logger.debug(f"Read interferogram phase from MiaplPy results from file {pl_inverted_file}.") + phase_linking_obj = BaseStack(file=pl_inverted_file, + logger=logger) + + pl_phase = ut.readPhasePatchwise( + stack_obj=phase_linking_obj, dataset_name="phase", + num_patches=config.general.num_patches, + cand_mask=cand_mask2, + point_id_img=point_id_img, logger=logger + ) + + # subset to time span + slc_stack_obj = slcStack(join(config.general.input_path, "slcStack.h5")) + slc_stack_obj.open() + time_mask = createTimeMaskFromDates( + start_date=config.preparation.start_date, + stop_date=config.preparation.end_date, + date_list=slc_stack_obj.dateList, + logger=logger + )[0] + pl_phase = pl_phase[:, time_mask] + + pl_ifgs = np.zeros((point2_obj.num_points, point2_obj.ifg_net_obj.num_ifgs), dtype=np.float32) + + c = 0 + for i, j in np.asarray(point1_obj.ifg_net_obj.ifg_list): + pl_ifgs[:, c] = np.angle(np.exp(1j * pl_phase[:, i]) * np.conjugate(np.exp(1j * pl_phase[:, j]))) + c += 1 + + # change only phase for good phase linking pixels and keep original phase for good tpc pixels + mask_pl = cand_mask_pl[cand_mask2] + point2_obj.phase[mask_pl] = pl_ifgs[mask_pl] + + logger.debug(f"Write second-order points to file {point2_obj.file_path}.") + point2_obj.writeToFile() + del ifg_stack_obj + + p1_aps_file = join(output_path, "p1_aps.h5") + logger.debug(f"Read APS for first order point from file {p1_aps_file}.") + aps1_obj = Points(file_path=p1_aps_file, logger=logger) + aps1_obj.open(input_path=config.general.input_path) + + aps_params_file = join(output_path, "aps_parameters.h5") + logger.debug(f"Read APS parameters from file {aps_params_file}.") + aps_params_obj = ApsParameters(file_path=aps_params_file, logger=logger) + aps_params_obj.open() + + filtering_method = aps_params_obj.method + + if filtering_method == "None": + logger.info("Skip APS filtering. To enable APS filtering, update the config file and rerun step 3.") + aps2_phase = np.zeros((point2_obj.num_points, point2_obj.ifg_net_obj.num_images)) + else: + logger.info("Apply APS filtering to second-order points.") + + if filtering_method == "kriging": + aps2_phase = applySpatialFilteringToP2(coord_map1=aps1_obj.coord_map, + residuals=aps_params_obj.phase, + coord_map2=point2_obj.coord_map, + model_name=aps_params_obj.model_name, + model_params=aps_params_obj.model_params, + logger=logger) + elif filtering_method == "simple": + aps2_phase = applySimpleInterpolationToP2(residuals=aps_params_obj.phase, + coord_map1=aps1_obj.coord_map, + coord_map2=point2_obj.coord_map, + logger=logger, + interp_method=aps_params_obj.params_obj.model_name) + + # TODO: delete this part after testing the aps estimation + # p2_aps_file = join(output_path, f"p2_coh{coh_value}_aps.h5") + # logger.info(f"Write APS for second-order points to file {p2_aps_file}.") + # aps2_obj = Points(file_path=p2_aps_file, logger=logger) + # aps2_obj.open( + # other_file_path=join(output_path, f"p2_coh{coh_value}_ifg_wr.h5"), + # input_path=config.general.input_path + # ) + # aps2_obj.phase = aps2_phase + # aps2_obj.writeToFile() + + logger.debug("Finished selecting second-order points.") + return point2_obj, aps2_phase diff --git a/sarvey/filtering.py b/sarvey/filtering.py index b4e3d280..411ea754 100644 --- a/sarvey/filtering.py +++ b/sarvey/filtering.py @@ -41,6 +41,122 @@ import sarvey.utils as ut +def applySpatialFilteringToP2(*, coord_map1, residuals, coord_map2, model_name, model_params, logger): + """ + Apply spatial filtering to second-order points using kriging. + + Parameters + ---------- + coord_map1 : np.ndarray + Map coordinates of the first-order points in UTM. + residuals : np.ndarray + Residuals of the first-order points that were used to estimate APS parameters. + coord_map2 : np.ndarray + Map coordinates of the second-order points in UTM. + model_name : str + Name of the variogram model used in APS estimation. + model_params : np.ndarray + Parameters of variogram model used in APS estimation. + logger : Logger + Logger instance for logging messages. + + Returns + ------- + aps2 : np.ndarray + Interpolated phase values for the second-order points. + """ + num_time = model_params.shape[1] + num_points2 = coord_map2.shape[0] + logger.debug(f"Start applying {num_time} time APS to {num_points2} second-order points using Kriging" + f" interpolation.") + + if model_name == 'stable': + logger.debug(f"Recreate {model_name} gs model using {model_params.shape[0]} parameters.") + model = gs.Stable(dim=2) + if model_params.shape[0] != 4: + error_msg = (f"Number of model parameters does not match the number of parameters required by " + f"{model_name} model") + logger.error(msg=error_msg) + raise ValueError(error_msg) + else: + error_msg = f"Model {model_name} not implemented yet." + logger.warning(msg=error_msg) + raise ValueError(error_msg) + + coord_point1 = [coord_map1[:, 0], coord_map1[:, 1]] + + aps2 = np.zeros((num_points2, num_time), dtype=np.float32) + + logger.debug(msg=f"Start looping for kriging parameters estimation and APS interpolation for {num_time} snapshots.") + for i in range(num_time): + try: + model = applyModelParams(model, model_params[:, i], logger) + except Exception as e: + warning_msg = f"Model parameters for time step {i} are not provided: {e}" + logger.warning(msg=warning_msg) + raise ValueError(warning_msg) + + field = residuals[:, i].astype(np.float32) + + # estimate kriging parameters based on the variogram model and the residuals of first order points + sk = gs.krige.Simple(model=model, cond_pos=coord_point1, cond_val=field) + + # evaluate the kriging model at new locations + aps2[:, i], var_sk_new = sk((coord_map2[:, 1], coord_map2[:, 0]), return_var=True) + + logger.debug(msg="Finished applying spatial filtering to second-order points using Kriging interpolation.") + return aps2 + + +def applySimpleInterpolationToP2(*, residuals: np.ndarray, coord_map1: np.ndarray, coord_map2: np.ndarray, + logger: Logger, interp_method: str = "linear",): + """ + Apply simple interpolation to second-order points. + + Parameters + ---------- + logger + residuals : np.ndarray + Residual phase values from the first-order points (size: num_points_p1 x num_ifgs) + coord_map1 : np.ndarray + Map coordinates of the points for which the residuals are given (size: num_points_p1 x 2) + coord_map2 : np.ndarray + Map coordinates of the second-order points. (size: num_points_p2 x 2) + logger: Logger + Logger instance for logging messages. + + interp_method : str + Interpolation method (default: "linear"; options: "linear", "cubic") + + Returns + ------- + interpolated_phase : np.ndarray + Atmospheric phase screen for the second order points (size: num_points_p2 x num_images) + """ + num_points2 = coord_map2.shape[0] + num_images = residuals.shape[1] + logger.debug(f"Start applying {num_images} APS to {num_points2} second-order points using" + f" {interp_method} interpolation.") + + aps2 = np.zeros((num_points2, num_images), dtype=np.float32) + + logger.debug(msg=f"Start looping for APS interpolation for {num_images} snapshots.") + for i in range(num_images): + aps2[:, i] = griddata(coord_map1, residuals[:, i], coord_map2, method=interp_method) + # interpolation with 'linear' or 'cubic' yields nan values for pixel that need to be extrapolated. + # interpolation with 'knn' solves this problem. + mask_extrapolate = np.isnan(aps2[:, i]) + aps2[mask_extrapolate, i] = griddata( + coord_map1, + residuals[:, i], + coord_map2[mask_extrapolate, :], + method='nearest' + ) + + logger.debug(msg=f"Finished applying spatial filtering to second-order points using {interp_method} interpolation.") + return aps2 + + def launchSpatialFiltering(parameters: tuple): """Launch_spatial_filtering. @@ -57,10 +173,8 @@ def launchSpatialFiltering(parameters: tuple): number of time steps residuals: np.ndarray residual phase (size: num_points x num_ifgs) - coord_utm1: np.ndarray - coordinates in UTM of the first-order points for which the residuals are given (size: num_points_p1 x 2) - coord_utm2: np.ndarray - coordinates in UTM of the new points which shall be interpolated (size: num_points_p2 x 2) + coord_map1: np.ndarray + map coordinates of the first-order points for which the residuals are given (size: num_points_p1 x 2) bins: np.ndarray bin edges for the variogram bool_plot: bool @@ -74,22 +188,28 @@ def launchSpatialFiltering(parameters: tuple): range of indices for the time series aps1: np.ndarray atmospheric phase screen for the known points (size: num_points_p1 x num_ifgs) - aps2: np.ndarray - atmospheric phase screen for the new points (size: num_points_p2 x num_ifgs) """ # Unpack the parameters - (idx_range, num_time, residuals, coord_utm1, coord_utm2, bins, bool_plot, logger) = parameters + (idx_range, num_time, residuals, coord_map1, bins, bool_plot, logger) = parameters - x = coord_utm1[:, 1] - y = coord_utm1[:, 0] - x_new = coord_utm2[:, 1] - y_new = coord_utm2[:, 0] + x = coord_map1[:, 1] + y = coord_map1[:, 0] - aps1 = np.zeros((coord_utm1.shape[0], num_time), dtype=np.float32) - aps2 = np.zeros((coord_utm2.shape[0], num_time), dtype=np.float32) + aps1 = np.zeros((coord_map1.shape[0], num_time), dtype=np.float32) prog_bar = ptime.progressBar(maxValue=num_time) + model_name = 'Stable' + if model_name == 'Stable': + model = gs.Stable(dim=2) + params = extractModelParams(model, logger) + model_params = np.array([np.full(num_time, param) for param in params], dtype=np.float32) + else: + # TODO: Other models to be added later + error_msg = f"Model {model_name} not implemented yet." + logger.error(msg=error_msg) + raise ValueError(error_msg) + for i in range(num_time): field = residuals[:, i].astype(np.float32) @@ -97,9 +217,9 @@ def launchSpatialFiltering(parameters: tuple): bin_center, vario = gs.vario_estimate(pos=[x, y], field=field, bin_edges=bins) # 2) fit model to empirical variogram - model = gs.Stable(dim=2) try: model.fit_variogram(x_data=bin_center, y_data=vario, nugget=True, max_eval=1500) + model_params[:, i] = extractModelParams(model=model, logger=logger) except RuntimeError as err: logger.error(msg="\nIMAGE {}: Not able to fit variogram! {}".format(idx_range[i], err)) if bool_plot: @@ -125,47 +245,12 @@ def launchSpatialFiltering(parameters: tuple): fld_sk, _ = sk((x, y), return_var=True) aps1[:, i] = fld_sk - # 5) evaluate the kriging model at NEW locations - fld_sk_new, var_sk_new = sk((x_new, y_new), return_var=True) - aps2[:, i] = fld_sk_new - prog_bar.update(value=i + 1, every=1, suffix='{}/{} images'.format(i + 1, num_time)) - # 5) show results - if bool_plot: - min_val = np.min(field) - max_val = np.max(field) - - fig, ax = plt.subplots(2, 2, figsize=[10, 5]) - - cur_ax = ax[0, 0] - sca1 = cur_ax.scatter(x, y, c=field, vmin=min_val, vmax=max_val) - plt.colorbar(sca1, ax=cur_ax, pad=0.03, shrink=0.5) - cur_ax.set_title("PS1 residuals") - - cur_ax = ax[0, 1] - cur_ax = model.plot(x_max=bin_center[-1], ax=cur_ax) - cur_ax.scatter(bin_center, vario) - cur_ax.set_xlabel("distance in [m]") - cur_ax.set_ylabel("semi-variogram") - - if coord_utm2 is not None: - cur_ax = ax[1, 0] - sca2 = cur_ax.scatter(x_new, y_new, c=fld_sk_new, vmin=min_val, vmax=max_val) - plt.colorbar(sca2, ax=cur_ax, pad=0.03, shrink=0.5) - cur_ax.set_title("PS2 prediction of atmospheric effect") - - cur_ax = ax[0, 1] - sca4 = cur_ax.scatter(x_new, y_new, c=var_sk_new) - plt.colorbar(sca4, ax=cur_ax, pad=0.03, shrink=0.5) - cur_ax.set_title("Variance of predicted atmospheric effect") + return idx_range, aps1, model_params - plt.close(fig) - return idx_range, aps1, aps2 - - -def estimateAtmosphericPhaseScreen(*, residuals: np.ndarray, coord_utm1: np.ndarray, coord_utm2: np.ndarray, +def estimateAtmosphericPhaseScreen(*, residuals: np.ndarray, coord_map1: np.ndarray, num_cores: int = 1, bool_plot: bool = False, logger: Logger) -> tuple[np.ndarray, np.ndarray]: """Estimate_atmospheric_phase_screen. @@ -177,10 +262,8 @@ def estimateAtmosphericPhaseScreen(*, residuals: np.ndarray, coord_utm1: np.ndar ---------- residuals: np.ndarray residual phase (size: num_points1 x num_images) - coord_utm1: np.ndarray + coord_map1: np.ndarray coordinates in UTM of the points for which the residuals are given (size: num_points1 x 2) - coord_utm2: np.ndarray - coordinates in UTM of the new points which shall be interpolated (size: num_points2 x 2) num_cores: int Number of cores bool_plot: bool @@ -192,8 +275,8 @@ def estimateAtmosphericPhaseScreen(*, residuals: np.ndarray, coord_utm1: np.ndar ------- aps1: np.ndarray atmospheric phase screen for the known points (size: num_points1 x num_images) - aps2: np.ndarray - atmospheric phase screen for the new points (size: num_points2 x num_images) + model_params: np.ndarray + model parameters for the variogram model (size: num_model_parameters x num_images) """ msg = "#" * 10 msg += " ESTIMATE ATMOSPHERIC PHASE SCREEN (KRIGING) " @@ -203,21 +286,20 @@ def estimateAtmosphericPhaseScreen(*, residuals: np.ndarray, coord_utm1: np.ndar start_time = time.time() num_points1 = residuals.shape[0] - num_points2 = coord_utm2.shape[0] num_time = residuals.shape[1] # can be either num_ifgs or num_images - bins = gs.variogram.standard_bins(pos=(coord_utm1[:, 1], coord_utm1[:, 0]), + bins = gs.variogram.standard_bins(pos=(coord_map1[:, 1], coord_map1[:, 0]), dim=2, latlon=False, mesh_type='unstructured', bin_no=30, max_dist=None) if num_cores == 1: - args = (np.arange(0, num_time), num_time, residuals, coord_utm1, coord_utm2, bins, bool_plot, logger) - _, aps1, aps2 = launchSpatialFiltering(parameters=args) + args = (np.arange(0, num_time), num_time, residuals, coord_map1, bins, bool_plot, logger) + _, aps1, model_params = launchSpatialFiltering(parameters=args) else: logger.info(msg="start parallel processing with {} cores.".format(num_cores)) pool = multiprocessing.Pool(processes=num_cores) aps1 = np.zeros((num_points1, num_time), dtype=np.float32) - aps2 = np.zeros((num_points2, num_time), dtype=np.float32) + model_params = np.zeros((4, num_time), dtype=np.float32) # 4 parameters for the Stable variogram model. num_cores = num_time if num_cores > num_time else num_cores # avoids having more samples than cores idx = ut.splitDatasetForParallelProcessing(num_samples=num_time, num_cores=num_cores) @@ -226,8 +308,7 @@ def estimateAtmosphericPhaseScreen(*, residuals: np.ndarray, coord_utm1: np.ndar idx_range, idx_range.shape[0], residuals[:, idx_range], - coord_utm1, - coord_utm2, + coord_map1, bins, False, logger) for idx_range in idx] @@ -235,18 +316,17 @@ def estimateAtmosphericPhaseScreen(*, residuals: np.ndarray, coord_utm1: np.ndar results = pool.map(func=launchSpatialFiltering, iterable=args) # retrieve results - for i, aps1_i, aps2_i in results: + for i, aps1_i, model_params_i in results: aps1[:, i] = aps1_i - aps2[:, i] = aps2_i + model_params[:, i] = model_params_i m, s = divmod(time.time() - start_time, 60) logger.debug(msg='time used: {:02.0f} mins {:02.1f} secs.\n'.format(m, s)) - return aps1, aps2 + return aps1, model_params -def simpleInterpolation(*, residuals: np.ndarray, coord_utm1: np.ndarray, coord_utm2: np.ndarray, - interp_method: str = "linear"): +def simpleInterpolation(*, residuals: np.ndarray, coord_map1: np.ndarray, interp_method: str = "linear"): """SimpleInterpolation. Simple interpolation of atmospheric phase screen using scipy's griddata function with options "linear" or "cubic". @@ -256,10 +336,8 @@ def simpleInterpolation(*, residuals: np.ndarray, coord_utm1: np.ndarray, coord_ ---------- residuals: np.ndarray residual phase (size: num_points x num_ifgs) - coord_utm1: np.ndarray - coordinates in UTM of the points for which the residuals are given (size: num_points_p1 x 2) - coord_utm2: np.ndarray - coordinates in UTM of the new points which shall be interpolated (size: num_points_p2 x 2) + coord_map1: np.ndarray + map coordinates of the points for which the residuals are given (size: num_points_p1 x 2) interp_method: str interpolation method (default: "linear"; options: "linear", "cubic") @@ -267,25 +345,78 @@ def simpleInterpolation(*, residuals: np.ndarray, coord_utm1: np.ndarray, coord_ ------- aps1: np.ndarray atmospheric phase screen for the known points (size: num_points_p1 x num_images) - aps2: np.ndarray - atmospheric phase screen for the new points (size: num_points_p2 x num_images) """ - num_points2 = coord_utm2.shape[0] num_images = residuals.shape[1] aps1 = np.zeros_like(residuals, dtype=np.float32) - aps2 = np.zeros((num_points2, num_images), dtype=np.float32) for i in range(num_images): - aps1[:, i] = griddata(coord_utm1, residuals[:, i], coord_utm1, method=interp_method) - aps2[:, i] = griddata(coord_utm1, residuals[:, i], coord_utm2, method=interp_method) - # interpolation with 'linear' or 'cubic' yields nan values for pixel that need to be extrapolated. - # interpolation with 'knn' solves this problem. - mask_extrapolate = np.isnan(aps2[:, i]) - aps2[mask_extrapolate, i] = griddata( - coord_utm1, - residuals[:, i], - coord_utm2[mask_extrapolate, :], - method='nearest' - ) + aps1[:, i] = griddata(coord_map1, residuals[:, i], coord_map1, method=interp_method) + + return aps1 - return aps1, aps2 + +def extractModelParams(model: gs.covmodel.models, logger: Logger): + """ + Extract parameters from the given gstools model. + + Parameters + ---------- + model : gs.covmodel.models + The model from which to extract parameters. + logger : Logger + Logging handler. + + Returns + ------- + tuple + A tuple containing the model parameters. + + Raises + ------ + ValueError + If the model is not implemented. + """ + if model.name == 'Stable': + params = (model.var, model.len_scale, model.nugget, model.alpha) + logger.debug(msg=f"Extract {len(params)} parameters from gs model {model.name}.") + # elif model.name == 'Gaussian': + # params = (model.var, model.len_scale, model.nugget) + else: + error_msg = f"Model {model.name} not implemented yet." + logger.error(msg=error_msg) + raise ValueError(error_msg) + return params + + +def applyModelParams(model: gs.covmodel.models, params: tuple, logger: Logger): + """ + Apply parameters to the given gstools model. + + Parameters + ---------- + model : gs.covmodel.models + The model to which parameters will be applied. + params : tuple + A tuple containing the model parameters. + logger : Logger + Logging handler. + + Returns + ------- + model + The model with applied parameters. + + Raises + ------ + ValueError + If the model is not implemented. + """ + if model.name == 'Stable': + model.var, model.len_scale, model.nugget, model.alpha = params + # elif model.name == 'Gaussian': # TODO: implement other models + # model.var, model.len_scale, model.nugget = params + else: + error_msg = f"Model {model.name} not implemented yet." + logger.error(msg=error_msg) + raise ValueError(error_msg) + return model diff --git a/sarvey/objects.py b/sarvey/objects.py index f2028c1a..aea428bd 100644 --- a/sarvey/objects.py +++ b/sarvey/objects.py @@ -35,8 +35,6 @@ import matplotlib.pyplot as plt import numpy as np from pyproj import Proj, CRS -from pyproj.aoi import AreaOfInterest -from pyproj.database import query_utm_crs_info from logging import Logger from cmcrameri import cm @@ -141,8 +139,8 @@ def plot(self, *, ax: plt.Axes = None, logger: Logger): return ax -class CoordinatesUTM: - """Coordinates in UTM for all pixels in the radar image.""" +class CoordinatesMap: + """Map Coordinates for all pixels in the radar image.""" def __init__(self, *, file_path: str, logger: Logger): """Init. @@ -155,7 +153,9 @@ def __init__(self, *, file_path: str, logger: Logger): Logging handler. """ self.file_path = file_path - self.coord_utm = None + self.coord_map = None + self.lat_0 = None + self.lon_0 = None self.logger = logger def prepare(self, *, input_path: str): @@ -170,19 +170,38 @@ def prepare(self, *, input_path: str): lat = readfile.read(input_path, datasetName='latitude')[0] lon = readfile.read(input_path, datasetName='longitude')[0] - log.info(msg="Transform coordinates from latitude and longitude (WGS84) to North and East (UTM).") + log.info(msg="Transform coordinates from latitude and longitude (WGS84) to North and East.") # noinspection PyTypeChecker - utm_crs_list = query_utm_crs_info( - datum_name="WGS 84", - area_of_interest=AreaOfInterest( - west_lon_degree=np.nanmin(lon.ravel()), - south_lat_degree=np.nanmin(lat.ravel()), - east_lon_degree=np.nanmax(lon.ravel()), - north_lat_degree=np.nanmax(lat.ravel())), - contains=True) - utm_crs = CRS.from_epsg(utm_crs_list[0].code) - lola2utm = Proj(utm_crs) - self.coord_utm = np.array(lola2utm(lon, lat)) + + log.info(msg="Transform coordinates from lat/lon (WGS84) to North/East (Transverse Mercator).") + # noinspection PyTypeChecker + + min_lat = np.nanmin(lat.ravel()) + max_lat = np.nanmax(lat.ravel()) + min_lon = np.nanmin(lon.ravel()) + max_lon = np.nanmax(lon.ravel()) + log.debug(f"WGS84 Longitude range: {min_lon} - {max_lon} Latitude range: {min_lat} - {max_lat}") + + self.lon_0 = (max_lon + min_lon) / 2 # Central_meridian + self.lat_0 = (max_lat + min_lat) / 2 # Latitude of origin + log.debug(f"Central meridian: {self.lon_0} Latitude of origin: {self.lat_0}") + + log.debug("Construct Transverse Mercator projection.") + map_crs = CRS.from_dict({ + "proj": "tmerc", # Transverse Mercator + "lat_0": self.lat_0, + "lon_0": self.lon_0, + "k": 1, + "x_0": 0, # False Easting + "y_0": 0, + "datum": "WGS84", + "units": "m", + "no_defs": True + }) + lola2map = Proj(map_crs) + self.coord_map = np.array(lola2map(lon, lat)) + log.debug(f"North coordinate range: {np.nanmin(self.coord_map[0, :])} - {np.nanmax(self.coord_map[0, :])}") + log.debug(f"East coordinate range: {np.nanmin(self.coord_map[1, :])} - {np.nanmax(self.coord_map[1, :])}") log.info(msg="write data to {}...".format(self.file_path)) @@ -190,12 +209,16 @@ def prepare(self, *, input_path: str): os.remove(self.file_path) with h5py.File(self.file_path, 'w') as f: - f.create_dataset('coord_utm', data=self.coord_utm) + f.create_dataset('coord_map', data=self.coord_map) + f.attrs["lon_0"] = self.lon_0 + f.attrs["lat_0"] = self.lat_0 def open(self): """Open.""" with h5py.File(self.file_path, 'r') as f: - self.coord_utm = f["coord_utm"][:] + self.coord_map = f["coord_map"][:] + self.lon_0 = f.attrs["lon_0"] + self.lat_0 = f.attrs["lat_0"] class BaseStack: @@ -449,7 +472,7 @@ def __init__(self, *, file_path: str, logger: Logger): Logging handler. """ self.ifg_net_obj = IfgNetwork() # use parent class here which doesn't know and care about 'star' or 'sb' - self.coord_utm = None + self.coord_map = None self.coord_lalo = None self.height = None self.slant_range = None @@ -461,7 +484,7 @@ def prepare(self, *, point_id: np.ndarray, coord_xy: np.ndarray, input_path: str """Assign point_id and radar coordinates to the object. Store the point_id and radar coordinates of the scatterers in the object (not file) and read further - attributes from external files (ifg_network.h5, slcStack.h5, geometryRadar.h5, coordinates_utm.h5). + attributes from external files (ifg_network.h5, slcStack.h5, geometryRadar.h5, coordinates_map.h5). Parameters ---------- @@ -521,7 +544,7 @@ def open(self, input_path: str, other_file_path: str = None): self.openExternalData(input_path=input_path) def openExternalData(self, *, input_path: str): - """Load data which is stored in slcStack.h5, geometryRadar.h5, ifg_network.h5 and coordinates_utm.h5.""" + """Load data which is stored in slcStack.h5, geometryRadar.h5, ifg_network.h5 and coordinates_map.h5.""" # 1) read IfgNetwork self.ifg_net_obj.open(path=join(dirname(self.file_path), "ifg_network.h5")) @@ -550,11 +573,11 @@ def openExternalData(self, *, input_path: str): self.height = height[mask].ravel() self.coord_lalo = np.array([lat[mask].ravel(), lon[mask].ravel()]).transpose() - # 4) read UTM coordinates - coord_utm_obj = CoordinatesUTM(file_path=join(dirname(self.file_path), "coordinates_utm.h5"), + # 4) read Map coordinates + coord_map_obj = CoordinatesMap(file_path=join(dirname(self.file_path), "coordinates_map.h5"), logger=self.logger) - coord_utm_obj.open() - self.coord_utm = coord_utm_obj.coord_utm[:, mask].transpose() + coord_map_obj.open() + self.coord_map = coord_map_obj.coord_map[:, mask].transpose() def createMask(self): """Create a mask. @@ -794,3 +817,96 @@ def open(self, *, input_path: str): self.demerr = f["demerr"][:] self.vel = f["vel"][:] self.gamma = f["gamma"][:] + + +class ApsParameters: + """APS model parameters.""" + + def __init__(self, *, file_path: str, logger: Logger): + """Init. + + Parameters + ---------- + file_path: str + path to filename + logger: Logger + Logging handler + """ + self.file_path = file_path + self.method = None + self.model_name = None + self.model_params = None + self.num_params = None + self.phase = None + self.logger = logger + self.logger.debug(f"ApsParameters initialized with file path: {file_path}") + + def prepare(self, *, method: str = "kriging", model_name: str = "stable", model_params: np.ndarray, + phase: np.ndarray): + """Prepare APS model parameters and store it into a file. + + Parameters + ---------- + method: str + method used for interpolation (kriging/simple) + model_name: str + name of the variogram model/interpolation type for kriging/simple + model_params: np.ndarray + Model parameters + phase: np.ndarray + Residual phase of P1 used for APS estimation + """ + self.logger.info(f"Preparing ApsParameters with method {method} and model {model_name}") + + if model_params is None or model_params.size == 0: + self.logger.debug("Model parameters are empty; defaulting to an empty list.") + model_params = np.array([]) + + self.method = method + self.model_name = model_name + self.model_params = model_params + self.num_params = model_params.shape[0] if model_params is not None and model_params.size > 0 else 0 + self.phase = phase + + self.logger.debug(f"Model parameters set with {self.num_params} parameters.") + self.logger.debug(f"Phase data shape: {phase.shape}") + + self.writeToFile() + + def open(self): + """Read data from file. + + Read stored information from already existing .h5 file. + """ + self.logger.info(f"Opening file: {self.file_path}") + + try: + with h5py.File(self.file_path, 'r') as f: + self.model_params = f["model_params"][:] + self.phase = f["phase"][:] + self.method = f.attrs["method"] + self.model_name = f.attrs["model_name"] + self.num_params = f.attrs["num_params"] + self.logger.debug(f"Loaded model name: {self.model_name}, number of parameters: {self.num_params}") + except Exception as e: + self.logger.error(f"Failed to open file {self.file_path}: {e}") + + def writeToFile(self): + """Write data to .h5 file.""" + self.logger.info(msg="write data to {}...".format(self.file_path)) + + try: + if exists(self.file_path): + self.logger.warning(f"File {self.file_path} exists and will be overwritten.") + os.remove(self.file_path) + + with h5py.File(self.file_path, 'w') as f: + f.create_dataset('model_params', data=self.model_params) + f.create_dataset('phase', data=self.phase) + f.attrs["method"] = self.method + f.attrs["model_name"] = self.model_name + f.attrs["num_params"] = self.num_params + + self.logger.info(f"Data successfully written to {self.file_path}.") + except Exception as e: + self.logger.error(f"Error writing data to {self.file_path}: {e}") diff --git a/sarvey/preparation.py b/sarvey/preparation.py index a926c4c8..cc2fbd19 100644 --- a/sarvey/preparation.py +++ b/sarvey/preparation.py @@ -38,7 +38,7 @@ from sarvey import viewer import sarvey.utils as ut -from sarvey.objects import CoordinatesUTM, AmplitudeImage, BaseStack, Points +from sarvey.objects import CoordinatesMap, AmplitudeImage, BaseStack, Points from sarvey.triangulation import PointNetworkTriangulation @@ -211,11 +211,11 @@ def selectPixels(*, path: str, selection_method: str, thrsh: float, # cmap = "lajolla" if grid_size is not None: # -> sparse pixel selection - coord_utm_obj = CoordinatesUTM(file_path=join(path, "coordinates_utm.h5"), logger=logger) - coord_utm_obj.open() - box_list = ut.createSpatialGrid(coord_utm_img=coord_utm_obj.coord_utm, - length=coord_utm_obj.coord_utm.shape[1], - width=coord_utm_obj.coord_utm.shape[2], + coord_map_obj = CoordinatesMap(file_path=join(path, "coordinates_map.h5"), logger=logger) + coord_map_obj.open() + box_list = ut.createSpatialGrid(coord_map_img=coord_map_obj.coord_map, + length=coord_map_obj.coord_map.shape[1], + width=coord_map_obj.coord_map.shape[2], grid_size=grid_size, logger=logger)[0] cand_mask_sparse = ut.selectBestPointsInGrid(box_list=box_list, quality=quality, sel_min=grid_min_val) @@ -260,7 +260,7 @@ def createArcsBetweenPoints(*, point_obj: Points, knn: int = None, max_arc_lengt arcs: np.ndarray Arcs of the triangulation containing the indices of the points for each arc. """ - triang_obj = PointNetworkTriangulation(coord_xy=point_obj.coord_xy, coord_utmxy=point_obj.coord_utm, logger=logger) + triang_obj = PointNetworkTriangulation(coord_xy=point_obj.coord_xy, coord_map_xy=point_obj.coord_map, logger=logger) if knn is not None: triang_obj.triangulateKnn(k=knn) diff --git a/sarvey/processing.py b/sarvey/processing.py index 47d66675..66c22473 100644 --- a/sarvey/processing.py +++ b/sarvey/processing.py @@ -2,7 +2,7 @@ # SARvey - A multitemporal InSAR time series tool for the derivation of displacements. # -# Copyright (C) 2021-2025 Andreas Piter (IPI Hannover, piter@ipi.uni-hannover.de) +# Copyright (C) 2021-2024 Andreas Piter (IPI Hannover, piter@ipi.uni-hannover.de) # # This software was developed together with FERN.Lab (fernlab@gfz-potsdam.de) in the context # of the SAR4Infra project with funds of the German Federal Ministry for Digital and @@ -39,11 +39,11 @@ from mintpy.utils.plot import auto_flip_direction from sarvey import viewer -from sarvey.densification import densifyNetwork +from sarvey.densification import densifyNetwork, selectP2 from sarvey.filtering import estimateAtmosphericPhaseScreen, simpleInterpolation from sarvey.ifg_network import (DelaunayNetwork, SmallBaselineYearlyNetwork, SmallTemporalBaselinesNetwork, SmallBaselineNetwork, StarNetwork) -from sarvey.objects import Network, Points, AmplitudeImage, CoordinatesUTM, NetworkParameter, BaseStack +from sarvey.objects import Network, Points, AmplitudeImage, CoordinatesMap, NetworkParameter, BaseStack, ApsParameters from sarvey.unwrapping import spatialParameterIntegration, \ parameterBasedNoisyPointRemoval, temporalUnwrapping, spatialUnwrapping, removeGrossOutliers from sarvey.preparation import createArcsBetweenPoints, selectPixels, createTimeMaskFromDates @@ -184,9 +184,9 @@ def runPreparation(self): logger=log) # store auxilliary datasets for faster access during processing - coord_utm_obj = CoordinatesUTM(file_path=join(self.path, "coordinates_utm.h5"), logger=self.logger) - coord_utm_obj.prepare(input_path=join(self.config.general.input_path, "geometryRadar.h5")) - del coord_utm_obj + coord_map_obj = CoordinatesMap(file_path=join(self.path, "coordinates_map.h5"), logger=self.logger) + coord_map_obj.prepare(input_path=join(self.config.general.input_path, "geometryRadar.h5")) + del coord_map_obj bmap_obj = AmplitudeImage(file_path=join(self.path, "background_map.h5")) bmap_obj.prepare(slc_stack_obj=slc_stack_obj, img=mean_amp_img, logger=self.logger) @@ -439,7 +439,7 @@ def runUnwrappingTimeAndSpace(self): if self.config.unwrapping.use_arcs_from_temporal_unwrapping: arcs = net_par_obj.arcs # use this to avoid unreliable connections. Takes a bit longer. else: - triang_obj = PointNetworkTriangulation(coord_xy=point_obj.coord_xy, coord_utmxy=point_obj.coord_utm, + triang_obj = PointNetworkTriangulation(coord_xy=point_obj.coord_xy, coord_map_xy=point_obj.coord_map, logger=self.logger) triang_obj.triangulateGlobal() arcs = triang_obj.getArcsFromAdjMat() @@ -497,7 +497,7 @@ def runUnwrappingSpace(self): arcs = net_par_obj.arcs # use this to avoid unreliable connections. Takes a bit longer. else: # re-triangulate with delaunay to make PUMA faster - triang_obj = PointNetworkTriangulation(coord_xy=point_obj.coord_xy, coord_utmxy=point_obj.coord_utm, + triang_obj = PointNetworkTriangulation(coord_xy=point_obj.coord_xy, coord_map_xy=point_obj.coord_map, logger=self.logger) triang_obj.triangulateGlobal() arcs = triang_obj.getArcsFromAdjMat() @@ -550,15 +550,12 @@ def runUnwrappingSpace(self): def runFiltering(self): """RunFiltering.""" - coh_value = int(self.config.filtering.coherence_p2 * 100) - # create output file which contains filtered phase time series point1_obj = Points(file_path=join(self.path, "p1_ts_filt.h5"), logger=self.logger) point1_obj.open( other_file_path=join(self.path, "p1_ts.h5"), input_path=self.config.general.input_path ) - p1_mask = point1_obj.createMask() # used later for selecting psCand2 when a spatial mask AOI is given. # select only pixels which have low phase noise and are well distributed mask = point1_obj.createMask() @@ -589,14 +586,14 @@ def runFiltering(self): plt.close(fig) # create grid - coord_utm_obj = CoordinatesUTM(file_path=join(self.path, "coordinates_utm.h5"), logger=self.logger) - coord_utm_obj.open() + coord_map_obj = CoordinatesMap(file_path=join(self.path, "coordinates_map.h5"), logger=self.logger) + coord_map_obj.open() # remove points based on threshold mask_thrsh = auto_corr_img <= self.config.filtering.max_temporal_autocorrelation auto_corr_img[~mask_thrsh] = np.inf - box_list, num_box = ut.createSpatialGrid(coord_utm_img=coord_utm_obj.coord_utm, length=point1_obj.length, + box_list, num_box = ut.createSpatialGrid(coord_map_img=coord_map_obj.coord_map, length=point1_obj.length, width=point1_obj.width, grid_size=self.config.filtering.grid_size, logger=self.logger) @@ -608,9 +605,6 @@ def runFiltering(self): self.logger.warning(msg=f"Only {num_p1_points_for_filtering} points for APS filtering selected. Filtering " f"results are probably not reliable. You can e.g. increase 'max_auto_corr' or try " f"to increase the number of first-order points during step 1 and 2.") - if num_p1_points_for_filtering == 0: - self.logger.error("No points selected for APS filtering.") - raise ValueError point_id_img = np.arange(0, point1_obj.length * point1_obj.width).reshape( (point1_obj.length, point1_obj.width)) @@ -640,165 +634,27 @@ def runFiltering(self): input_path=self.config.general.input_path ) - # select second-order points - cand_mask2 = selectPixels( - path=self.path, selection_method="temp_coh", - thrsh=self.config.filtering.coherence_p2, - grid_size=None, bool_plot=True, - logger=self.logger - ) # first-order points are included in second-order points - - if self.config.phase_linking.use_phase_linking_results: - # read PL results - pl_coh = readfile.read(join(self.config.phase_linking.inverted_path, "phase_series.h5"), - datasetName='temporalCoherence')[0] - pl_coh = pl_coh[1, :, :] - siblings = readfile.read(join(self.config.phase_linking.inverted_path, "phase_series.h5"), - datasetName='shp')[0] - - if self.config.phase_linking.use_ps: - mask_ps = readfile.read(self.config.phase_linking.mask_ps_file, - datasetName='mask')[0].astype(np.bool_) - cand_mask_pl = (pl_coh > self.config.filtering.coherence_p2) | mask_ps - else: - cand_mask_pl = (pl_coh > self.config.filtering.coherence_p2) - # remove ps candidates, because the ps detection strategy in miaplpy seems to be biased. - cand_mask_pl[siblings <= self.config.phase_linking.num_siblings] = False - - if self.config.phase_linking.mask_phase_linking_file is not None: - path_mask_pl_aoi = join(self.config.phase_linking.mask_phase_linking_file) - self.logger.info(msg="load mask for area of interest from: {}.".format(path_mask_pl_aoi)) - mask_pl_aoi = readfile.read(path_mask_pl_aoi, datasetName='mask')[0].astype(np.bool_) - - fig = plt.figure(figsize=(15, 5)) - ax = fig.add_subplot() - ax.imshow(mask_pl_aoi, cmap=cmc.cm.cmaps["grayC"], alpha=0.5, zorder=10, vmin=0, vmax=1) - bmap_obj.plot(ax=ax, logger=self.logger) - coord_xy = np.array(np.where(cand_mask_pl)).transpose() - val = np.ones_like(cand_mask_pl) - sc = ax.scatter(coord_xy[:, 1], coord_xy[:, 0], c=val[cand_mask_pl], s=0.5, - cmap=cmc.cm.cmaps["lajolla_r"], - vmin=1, vmax=2) # set min, max to ensure that points are yellow - cbar = plt.colorbar(sc, pad=0.03, shrink=0.5) - cbar.ax.set_visible(False) # make size of axis consistent with all others - plt.tight_layout() - plt.title("Mask for phase linking points") - fig.savefig(join(self.path, "pic", "step_3_mask_p2_coh{}_phase_linking.png".format(coh_value)), dpi=300) - plt.close(fig) - - # mask points after plotting, so that available coherent points are visible in figure - cand_mask_pl[~mask_pl_aoi] = False - - # combine phase linking coherence with TPC cand_mask2 - cand_mask2 = cand_mask2 | cand_mask_pl - - mask_valid_area = ut.detectValidAreas(bmap_obj=bmap_obj, logger=self.logger) - - if self.config.filtering.mask_p2_file is not None: - path_mask_aoi = join(self.config.filtering.mask_p2_file) - self.logger.info(msg="load mask for area of interest from: {}.".format(path_mask_aoi)) - mask_aoi = readfile.read(path_mask_aoi, datasetName='mask')[0].astype(np.bool_) - mask_valid_area &= mask_aoi - # todo: add unstable points from p1 for densification - else: - self.logger.info(msg="No mask for area of interest given.") - - cand_mask2[p1_mask] = True # add all selected 1.order points to avoid spatial gaps in 2D unwrapping - # cand_mask2[cand_mask_sparse] = True # add only stable points from 1.order points - - cand_mask2 &= mask_valid_area - - fig = plt.figure(figsize=(15, 5)) - ax = fig.add_subplot() - ax.imshow(mask_valid_area, cmap=cmc.cm.cmaps["grayC"], alpha=0.5, zorder=10, vmin=0, vmax=1) - bmap_obj.plot(ax=ax, logger=self.logger) - coord_xy = np.array(np.where(cand_mask2)).transpose() - val = np.ones_like(cand_mask2) - sc = ax.scatter(coord_xy[:, 1], coord_xy[:, 0], c=val[cand_mask2], s=0.5, cmap=cmc.cm.cmaps["lajolla_r"], - vmin=0, vmax=10) # set min, max to ensure that points are yellow - cbar = plt.colorbar(sc, pad=0.03, shrink=0.5) - cbar.ax.set_visible(False) # make size of axis consistent with all others - plt.tight_layout() - plt.title("Mask for dense point set") - fig.savefig(join(self.path, "pic", "step_3_mask_p2_coh{}.png".format(coh_value)), dpi=300) - plt.close(fig) - - point2_obj = Points(file_path=join(self.path, "p2_coh{}_ifg_wr.h5".format(coh_value)), logger=self.logger) - coord_xy = np.array(np.where(cand_mask2)).transpose() - point_id2 = point_id_img[cand_mask2] - point2_obj.prepare( - point_id=point_id2, - coord_xy=coord_xy, - input_path=self.config.general.input_path - ) - - ifg_stack_obj = BaseStack(file=join(self.path, "ifg_stack.h5"), logger=self.logger) - - point2_obj.phase = ut.readPhasePatchwise(stack_obj=ifg_stack_obj, dataset_name="ifgs", - num_patches=self.config.general.num_patches, cand_mask=cand_mask2, - point_id_img=point_id_img, logger=self.logger) - - if self.config.phase_linking.use_phase_linking_results: - self.logger.info(msg="read phase from MiaplPy results...") - phase_linking_obj = BaseStack( - file=join(self.config.phase_linking.inverted_path, "phase_series.h5"), - logger=self.logger - ) - - pl_phase = ut.readPhasePatchwise( - stack_obj=phase_linking_obj, dataset_name="phase", - num_patches=self.config.general.num_patches, - cand_mask=cand_mask2, - point_id_img=point_id_img, logger=self.logger - ) - - # subset to time span - slc_stack_obj = slcStack(join(self.config.general.input_path, "slcStack.h5")) - slc_stack_obj.open() - time_mask = createTimeMaskFromDates( - start_date=self.config.preparation.start_date, - stop_date=self.config.preparation.end_date, - date_list=slc_stack_obj.dateList, - logger=self.logger - )[0] - pl_phase = pl_phase[:, time_mask] - - pl_ifgs = np.zeros((point2_obj.num_points, point2_obj.ifg_net_obj.num_ifgs), dtype=np.float32) - - c = 0 - for i, j in np.asarray(point1_obj.ifg_net_obj.ifg_list): - pl_ifgs[:, c] = np.angle(np.exp(1j * pl_phase[:, i]) * np.conjugate(np.exp(1j * pl_phase[:, j]))) - c += 1 - - # change only phase for good phase linking pixels and keep original phase for good tpc pixels - mask_pl = cand_mask_pl[cand_mask2] - point2_obj.phase[mask_pl] = pl_ifgs[mask_pl] - - point2_obj.writeToFile() - del point2_obj, ifg_stack_obj - - aps2_obj = Points(file_path=join(self.path, "p2_coh{}_aps.h5".format(coh_value)), logger=self.logger) - aps2_obj.open( - other_file_path=join(self.path, "p2_coh{}_ifg_wr.h5".format(coh_value)), - input_path=self.config.general.input_path - ) + # select second-order points moved to selectP2() if self.config.filtering.apply_aps_filtering: # spatial filtering of points with linear motion only (no non-linear motion) if self.config.filtering.interpolation_method == "kriging": - aps1_phase, aps2_phase = estimateAtmosphericPhaseScreen( + method = "kriging" + model_name = "stable" # To be integrated in the config? + aps1_phase, aps_model_params = estimateAtmosphericPhaseScreen( residuals=phase_for_aps_filtering, - coord_utm1=point1_obj.coord_utm, - coord_utm2=aps2_obj.coord_utm, + coord_map1=point1_obj.coord_map, num_cores=self.config.general.num_cores, bool_plot=False, logger=self.logger ) else: - aps1_phase, aps2_phase = simpleInterpolation( + method = "simple" + model_name = "linear" # To be integrated in the config? + aps_model_params = None + aps1_phase = simpleInterpolation( residuals=phase_for_aps_filtering, - coord_utm1=point1_obj.coord_utm, - coord_utm2=aps2_obj.coord_utm, + coord_map1=point1_obj.coord_map2, interp_method=self.config.filtering.interpolation_method ) else: @@ -806,23 +662,33 @@ def runFiltering(self): msg += " SKIP ATMOSPHERIC FILTERING! " msg += "#" * 10 self.logger.info(msg=msg) + method = "None" + model_name = "None" + aps_model_params = None num_points1 = phase_for_aps_filtering.shape[0] - num_points2 = aps2_obj.coord_utm.shape[0] num_time = phase_for_aps_filtering.shape[1] aps1_phase = np.zeros((num_points1, num_time), dtype=np.float32) - aps2_phase = np.zeros((num_points2, num_time), dtype=np.float32) + + # save APS parameters for later use in step 4 + aps_file = join(self.path, "aps_parameters.h5") + self.logger.debug(f"Prepare Aps Parameters file {aps_file}") + aps_params_obj = ApsParameters(file_path=aps_file, logger=self.logger) + aps_params_obj.prepare(method=method, model_name=model_name, model_params=aps_model_params, + phase=phase_for_aps_filtering) point1_obj.phase -= aps1_phase point1_obj.writeToFile() aps1_obj.phase = aps1_phase - aps2_obj.phase = aps2_phase aps1_obj.writeToFile() - aps2_obj.writeToFile() def runDensificationTimeAndSpace(self): """RunDensificationTimeAndSpace.""" - coh_value = int(self.config.filtering.coherence_p2 * 100) + coherence_p2 = self.config.densification.coherence_p2 + coh_value = int(self.config.densification.coherence_p2 * 100) + self.logger.info(f"Select second-order points using coherence threshold {coherence_p2:.2f}.") + + _, aps2_phase = selectP2(output_path=self.path, config=self.config, logger=self.logger) point2_obj = Points(file_path=join(self.path, "p2_coh{}_ifg_unw.h5".format(coh_value)), logger=self.logger) point2_obj.open( @@ -842,10 +708,7 @@ def runDensificationTimeAndSpace(self): aps1_obj = Points(file_path=join(self.path, "p1_aps.h5"), logger=self.logger) aps1_obj.open(input_path=self.config.general.input_path) - aps2_obj = Points(file_path=join(self.path, "p2_coh{}_aps.h5".format(coh_value)), logger=self.logger) - aps2_obj.open(input_path=self.config.general.input_path) - - if self.config.filtering.mask_p2_file is None: + if self.config.densification.mask_p2_file is None: """ overview of points contained in the *_obj (unstable p1 means: p1 which were not used in atmospheric filtering) @@ -862,16 +725,16 @@ def runDensificationTimeAndSpace(self): mask_unstable_p1 = p1_mask & (~aps1_mask) unstable_p1_id = point_id_img[np.where(mask_unstable_p1)] - mask_unstable_p1_in_p2 = np.ones((aps2_obj.num_points,), dtype=np.bool_) - for p in aps2_obj.point_id: + mask_unstable_p1_in_p2 = np.ones((point2_obj.num_points,), dtype=np.bool_) + for p in point2_obj.point_id: if p not in unstable_p1_id: - mask_unstable_p1_in_p2[aps2_obj.point_id == p] = False + mask_unstable_p1_in_p2[point2_obj.point_id == p] = False # add unstable p1 from aps2 to aps1 aps1_obj.addPointsFromObj( - new_point_id=aps2_obj.point_id[mask_unstable_p1_in_p2], - new_coord_xy=aps2_obj.coord_xy[mask_unstable_p1_in_p2, :], - new_phase=aps2_obj.phase[mask_unstable_p1_in_p2, :], + new_point_id=point2_obj.point_id[mask_unstable_p1_in_p2], + new_coord_xy=point2_obj.coord_xy[mask_unstable_p1_in_p2, :], + new_phase=aps2_phase[mask_unstable_p1_in_p2, :], new_num_points=mask_unstable_p1_in_p2[mask_unstable_p1_in_p2].shape[0], input_path=self.config.general.input_path ) @@ -880,8 +743,10 @@ def runDensificationTimeAndSpace(self): p2_mask = point2_obj.createMask() mask_only_p2 = p2_mask & (~p1_mask) keep_id = point_id_img[np.where(mask_only_p2)] + + mask = np.isin(point2_obj.point_id, keep_id) + aps2_phase = aps2_phase[mask, :] point2_obj.removePoints(keep_id=keep_id, input_path=self.config.general.input_path) - aps2_obj.removePoints(keep_id=keep_id, input_path=self.config.general.input_path) else: """ @@ -917,13 +782,15 @@ def runDensificationTimeAndSpace(self): p2_mask = point2_obj.createMask() mask_p2 = ~(p1_mask & p2_mask) & p2_mask p2_id = point_id_img[np.where(mask_p2)] + + mask = np.isin(point2_obj.point_id, p2_id) + aps2_phase = aps2_phase[mask, :] point2_obj.removePoints(keep_id=p2_id, input_path=self.config.general.input_path) - aps2_obj.removePoints(keep_id=p2_id, input_path=self.config.general.input_path) # return to ifg-space a_ifg = point2_obj.ifg_net_obj.getDesignMatrix() aps1_ifg_phase = np.matmul(a_ifg, aps1_obj.phase.T).T - aps2_ifg_phase = np.matmul(a_ifg, aps2_obj.phase.T).T + aps2_ifg_phase = np.matmul(a_ifg, aps2_phase.T).T # correct for APS point2_obj.phase = np.angle(np.exp(1j * point2_obj.phase) * np.conjugate(np.exp(1j * aps2_ifg_phase))) @@ -1012,7 +879,7 @@ def runDensificationTimeAndSpace(self): wr_phase = point2_obj.phase wr_res_phase = np.angle(np.exp(1j * wr_phase) * np.conjugate(np.exp(1j * pred_phase))) - triang_obj = PointNetworkTriangulation(coord_xy=point2_obj.coord_xy, coord_utmxy=None, logger=self.logger) + triang_obj = PointNetworkTriangulation(coord_xy=point2_obj.coord_xy, coord_map_xy=None, logger=self.logger) triang_obj.triangulateGlobal() arcs = triang_obj.getArcsFromAdjMat() @@ -1052,7 +919,11 @@ def runDensificationTimeAndSpace(self): def runDensificationSpace(self): """RunDensification.""" - coh_value = int(self.config.filtering.coherence_p2 * 100) + coherence_p2 = self.config.densification.coherence_p2 + coh_value = int(coherence_p2 * 100) + + self.logger.info(f"Select second-order points using coherence threshold {coherence_p2:.2f}.") + _, aps2_phase = selectP2(output_path=self.path, config=self.config, logger=self.logger) point_obj = Points(file_path=join(self.path, "p2_coh{}_ifg_unw.h5".format(coh_value)), logger=self.logger) point_obj.open( @@ -1060,18 +931,15 @@ def runDensificationSpace(self): input_path=self.config.general.input_path ) # open wr phase - aps2_obj = Points(file_path=join(self.path, "p2_coh{}_aps.h5".format(coh_value)), logger=self.logger) - aps2_obj.open(input_path=self.config.general.input_path) - # return to ifg-space a_ifg = point_obj.ifg_net_obj.getDesignMatrix() - aps2_ifg_phase = np.matmul(a_ifg, aps2_obj.phase.T).T + aps2_ifg_phase = np.matmul(a_ifg, aps2_phase.T).T # correct for APS point_obj.phase = np.angle(np.exp(1j * point_obj.phase) * np.conjugate(np.exp(1j * aps2_ifg_phase))) - triang_obj = PointNetworkTriangulation(coord_xy=point_obj.coord_xy, coord_utmxy=None, logger=self.logger) - triang_obj.triangulateGlobal() # if coord_utm is not given, only global delaunay and knn can be calculated + triang_obj = PointNetworkTriangulation(coord_xy=point_obj.coord_xy, coord_map_xy=None, logger=self.logger) + triang_obj.triangulateGlobal() # if coord_map is not given, only global delaunay and knn can be calculated arcs = triang_obj.getArcsFromAdjMat() unw_phase = spatialUnwrapping(num_ifgs=point_obj.ifg_net_obj.num_ifgs, diff --git a/sarvey/sarvey_export.py b/sarvey/sarvey_export.py index 83144c0d..25d9fdeb 100755 --- a/sarvey/sarvey_export.py +++ b/sarvey/sarvey_export.py @@ -39,15 +39,13 @@ import numpy as np import pandas as pd import geopandas as gpd -from pyproj import CRS -from pyproj.aoi import AreaOfInterest -from pyproj.database import query_utm_crs_info +from pyproj import CRS, Transformer from shapely import Point from shapely.errors import ShapelyDeprecationWarning from sarvey.config import loadConfiguration from sarvey.console import showLogoSARvey -from sarvey.objects import Points +from sarvey.objects import Points, CoordinatesMap import sarvey.utils as ut from sarvey.geolocation import calculateGeolocationCorrection @@ -83,7 +81,7 @@ def exportDataToGisFormat(*, file_path: str, output_path: str, input_path: str, vel, demerr, _, coherence, omega, _ = ut.estimateParameters(obj=point_obj, ifg_space=False) - stc = ut.spatiotemporalConsistency(coord_utm=point_obj.coord_utm, phase=point_obj.phase, + stc = ut.spatiotemporalConsistency(coord_map=point_obj.coord_map, phase=point_obj.phase, wavelength=point_obj.wavelength) point_obj.phase *= point_obj.wavelength / (4 * np.pi) # in [m] @@ -98,19 +96,6 @@ def exportDataToGisFormat(*, file_path: str, output_path: str, input_path: str, # transform into meters defo_ts *= 1000 # in [mm] - utm_crs_list = query_utm_crs_info( - datum_name="WGS 84", - area_of_interest=AreaOfInterest( - west_lon_degree=point_obj.coord_lalo[:, 1].min(), - south_lat_degree=point_obj.coord_lalo[:, 0].min(), - east_lon_degree=point_obj.coord_lalo[:, 1].max(), - north_lat_degree=point_obj.coord_lalo[:, 0].max(), - ), - contains=True - ) - - utm_epsg = utm_crs_list[0].code - dates = ["D{}".format(date).replace("-", "") for date in point_obj.ifg_net_obj.dates] dates = dates[:point_obj.phase.shape[1]] # remove dates which were not processed @@ -134,9 +119,38 @@ def exportDataToGisFormat(*, file_path: str, output_path: str, input_path: str, coord_correction = 0 logger.info("geolocation correction skipped.") - coord_utm = point_obj.coord_utm - coord_utm += coord_correction - df_points['coord'] = (coord_utm).tolist() + coord_map = point_obj.coord_map + coord_map += coord_correction + + # reconstruct Transverse Mercator projection + coord_map_obj = CoordinatesMap(file_path=join(dirname(file_path), "coordinates_map.h5"), logger=logger) + coord_map_obj.open() + lon_0 = coord_map_obj.lon_0 + lat_0 = coord_map_obj.lat_0 + + map_crs = CRS.from_dict({ + "proj": "tmerc", + "lat_0": lat_0, + "lon_0": lon_0, + "k": 1, + "x_0": 0, + "y_0": 0, + "datum": "WGS84", + "units": "m", + "no_defs": True + }) + + # construct Transverse Mercator to WGS84 transformer + wgs_epsg = "4326" + transformer = Transformer.from_crs(map_crs, wgs_epsg) + + lat, lon = transformer.transform(coord_map[:, 0], coord_map[:, 1]) + logger.debug(f"WGS84 Lon range: {np.min(lon):.6f} - {np.max(lon):.6f}") + logger.debug(f"WGS84 Lat range: {np.min(lat):.6f} - {np.max(lat):.6f}") + lonlat = np.vstack((lon, lat)).T + + logger.info('Construct output dataframe.') + df_points['coord'] = (lonlat).tolist() df_points['coord'] = df_points['coord'].apply(Point) df_points.insert(0, 'point_id', point_obj.point_id.tolist()) df_points.insert(1, 'velocity', vel * 1000) # in [mm] @@ -153,7 +167,7 @@ def exportDataToGisFormat(*, file_path: str, output_path: str, input_path: str, df_points[date] = defo_ts[:, i] gdf_points = gpd.GeoDataFrame(df_points, geometry='coord') - gdf_points = gdf_points.set_crs(CRS.from_epsg(utm_epsg)) + gdf_points = gdf_points.set_crs(CRS.from_epsg(wgs_epsg)) logger.info(msg="write to file.") gdf_points.to_file(output_path) diff --git a/sarvey/sarvey_mti.py b/sarvey/sarvey_mti.py index b770b918..5e1e810c 100755 --- a/sarvey/sarvey_mti.py +++ b/sarvey/sarvey_mti.py @@ -104,7 +104,7 @@ def run(*, config: Config, args: argparse.Namespace, logger: Logger): config_section_default=config_default_dict["preparation"], logger=logger) proc_obj.runPreparation() - required_files = ["background_map.h5", "coordinates_utm.h5", "ifg_network.h5", "ifg_stack.h5", + required_files = ["background_map.h5", "coordinates_map.h5", "ifg_network.h5", "ifg_stack.h5", "temporal_coherence.h5"] if 1 in steps: @@ -147,8 +147,7 @@ def run(*, config: Config, args: argparse.Namespace, logger: Logger): config_section_default=config_default_dict["filtering"], logger=logger) proc_obj.runFiltering() - coh_value = int(config.filtering.coherence_p2 * 100) - required_files.extend(["p1_aps.h5", f"p2_coh{coh_value}_ifg_wr.h5", f"p2_coh{coh_value}_aps.h5"]) + required_files.extend(["p1_aps.h5", "aps_parameters.h5"]) if 4 in steps: checkIfRequiredFilesExist( @@ -295,9 +294,9 @@ def main(iargs=None): if config.consistency_check.mask_p1_file is not None: config.consistency_check.mask_p1_file = os.path.abspath( join(args.workdir, config.consistency_check.mask_p1_file)) - if config.filtering.mask_p2_file is not None: - config.filtering.mask_p2_file = os.path.abspath( - join(args.workdir, config.filtering.mask_p2_file)) + if config.densification.mask_p2_file is not None: + config.densification.mask_p2_file = os.path.abspath( + join(args.workdir, config.densification.mask_p2_file)) # create all necessary directories if not os.path.exists(config.general.output_path): diff --git a/sarvey/sarvey_plot.py b/sarvey/sarvey_plot.py index 0f0f75d1..65c866da 100755 --- a/sarvey/sarvey_plot.py +++ b/sarvey/sarvey_plot.py @@ -142,7 +142,7 @@ def plotMap(*, obj_name: str, save_path: str, interactive: bool = False, input_p else: plt.close(plt.gcf()) - stc = ut.spatiotemporalConsistency(coord_utm=point_obj.coord_utm, phase=point_obj.phase, + stc = ut.spatiotemporalConsistency(coord_map=point_obj.coord_map, phase=point_obj.phase, wavelength=point_obj.wavelength, min_dist=50, max_dist=np.inf, knn=40) diff --git a/sarvey/triangulation.py b/sarvey/triangulation.py index c70f62e3..e59d08d9 100644 --- a/sarvey/triangulation.py +++ b/sarvey/triangulation.py @@ -42,15 +42,15 @@ class PointNetworkTriangulation: """PointNetworkTriangulation.""" - def __init__(self, *, coord_xy: np.ndarray, coord_utmxy: Optional[np.ndarray], logger: Logger): + def __init__(self, *, coord_xy: np.ndarray, coord_map_xy: Optional[np.ndarray], logger: Logger): """Triangulate points in space based on distance. Parameters ---------- coord_xy: np.ndarray Radar coordinates of the points. - coord_utmxy: np.ndarray - UTM coordinates of the points. + coord_map_xy: np.ndarray + map coordinates of the points. logger: Logger Logging handler. """ @@ -62,9 +62,9 @@ def __init__(self, *, coord_xy: np.ndarray, coord_utmxy: Optional[np.ndarray], l # create network afterwards once. reduces time. self.adj_mat = lil_matrix((num_points, num_points), dtype=np.bool_) - if coord_utmxy is not None: + if coord_map_xy is not None: logger.info(msg="create distance matrix between all points...") - self.dist_mat = distance_matrix(coord_utmxy, coord_utmxy) + self.dist_mat = distance_matrix(coord_map_xy, coord_map_xy) # todo: check out alternatives: # scipy.spatial.KDTree.sparse_distance_matrix else: # if only global delaunay shall be computed without memory issues diff --git a/sarvey/utils.py b/sarvey/utils.py index b811719b..18261139 100644 --- a/sarvey/utils.py +++ b/sarvey/utils.py @@ -430,8 +430,8 @@ def splitImageIntoBoxesRngAz(*, length: int, width: int, num_box_az: int, num_bo num_box: int number of boxes """ - y_step = int(np.rint((length / num_box_rng) / 10) * 10) - x_step = int(np.rint((width / num_box_az) / 10) * 10) + y_step = int(length / num_box_rng) + x_step = int(width / num_box_az) box_list = [] y0 = 0 @@ -548,13 +548,13 @@ def splitDatasetForParallelProcessing(*, num_samples: int, num_cores: int): return idx -def createSpatialGrid(*, coord_utm_img: np.ndarray, length: int, width: int, grid_size: int, logger: Logger): +def createSpatialGrid(*, coord_map_img: np.ndarray, length: int, width: int, grid_size: int, logger: Logger): """Create a spatial grid over the image. Parameters ---------- - coord_utm_img: np.ndarray - coordinates of all image pixels in UTM + coord_map_img: np.ndarray + map coordinates of all image pixels length: int number of pixels in length of the image width: int @@ -571,9 +571,9 @@ def createSpatialGrid(*, coord_utm_img: np.ndarray, length: int, width: int, gri num_box: int actual number of boxes created by the function """ - p0 = coord_utm_img[:, 0, 0] - p1 = coord_utm_img[:, 0, -1] - p2 = coord_utm_img[:, -1, 0] + p0 = coord_map_img[:, 0, 0] + p1 = coord_map_img[:, 0, -1] + p2 = coord_map_img[:, -1, 0] dist_width = np.linalg.norm(p0 - p1) dist_length = np.linalg.norm(p0 - p2) @@ -632,14 +632,14 @@ def selectBestPointsInGrid(*, box_list: list, quality: np.ndarray, sel_min: bool return cand_mask_sparse -def spatiotemporalConsistency(*, coord_utm: np.ndarray, phase: np.ndarray, wavelength: float, min_dist: int = 15, +def spatiotemporalConsistency(*, coord_map: np.ndarray, phase: np.ndarray, wavelength: float, min_dist: int = 15, max_dist: float = np.inf, knn: int = 50): """Spatiotemporal consistency proposed by Hanssen et al. (2008) and implemented in DePSI (van Leijen, 2014). Parameters ---------- - coord_utm: np.ndarray - UTM coordinates of the points + coord_map: np.ndarray + map coordinates of the points phase: np.ndarray phase time series of the points wavelength: float @@ -659,12 +659,12 @@ def spatiotemporalConsistency(*, coord_utm: np.ndarray, phase: np.ndarray, wavel from scipy.spatial import KDTree num_samples, num_time = phase.shape - tree = KDTree(data=coord_utm) + tree = KDTree(data=coord_map) stc = np.zeros((num_samples,), np.float64) for p in range(num_samples): - dist, idx = tree.query([coord_utm[p, 0], coord_utm[p, 1]], k=knn) + dist, idx = tree.query([coord_map[p, 0], coord_map[p, 1]], k=knn) mask = (dist < max_dist) & (dist > min_dist) & (dist != 0) rho = list() for i in idx[mask]: diff --git a/setup.py b/setup.py index 31f02e78..f483d483 100644 --- a/setup.py +++ b/setup.py @@ -43,7 +43,7 @@ exec(version_file.read(), version) req = [ - "cython", "numpy", "pyproj", "matplotlib", "numba", "scipy", + "cython", "numpy<=1.26", "pyproj", "matplotlib", "numba", "scipy", "mintpy", "h5py", "overpy", "miaplpy", "gstools", "shapely", "pandas", "geopandas", "pymaxflow", "pillow", "pydantic<=1.10.10", "importlib_resources", "kamui", "json5", "cmcrameri" ] @@ -63,6 +63,8 @@ req_dev = ['twine'] + req_setup + req_test + req_doc + req_lint +extra_req = ["gdal"] + setup( author="Andreas Piter", author_email='piter@ipi.uni-hannover.de', diff --git a/tests/CI_docker/build_sarvey_testsuite_image.sh b/tests/CI_docker/build_sarvey_testsuite_image.sh index c08ad9e2..04485541 100755 --- a/tests/CI_docker/build_sarvey_testsuite_image.sh +++ b/tests/CI_docker/build_sarvey_testsuite_image.sh @@ -12,7 +12,6 @@ print(version["__version__"]) ' version=`python -c "$python_script"` tag="sarvey_ci:$version" -gitlab_runner="sarvey_gitlab_CI_runner" echo "#### Build runner docker image" if [[ "$(docker images ${tag} | grep ${tag} 2> /dev/null)" != "" ]]; then diff --git a/tests/CI_docker/context/entrypoint.sh b/tests/CI_docker/context/entrypoint.sh new file mode 100644 index 00000000..94325667 --- /dev/null +++ b/tests/CI_docker/context/entrypoint.sh @@ -0,0 +1,26 @@ +#!/bin/bash +set -e + +# Check for necessary environment variables +if [ -z "$RUNNER_TOKEN" ] || [ -z "$RUNNER_REPO_URL" ]; then + echo "Error: RUNNER_TOKEN and RUNNER_REPO_URL environment variables must be set." + exit 1 +fi + +# Configure the runner if it hasn’t been configured already +if [ ! -f ".runner" ]; then + echo "Configuring the GitHub Actions runner..." + ./config.sh --unattended \ + --replace \ + --labels self-hosted \ + --url "$RUNNER_REPO_URL" \ + --token "$RUNNER_TOKEN" \ + --name "${RUNNER_NAME:-docker-runner}" || exit 1 +fi + +# Trap SIGTERM and deregister the runner on shutdown if needed +trap './config.sh remove --token "$RUNNER_TOKEN"; exit 0' SIGTERM + +echo "Starting the GitHub Actions runner..." +exec ./run.sh + diff --git a/tests/CI_docker/context/sarvey_ci.docker b/tests/CI_docker/context/sarvey_ci.docker index 51c4a85b..d34dcc45 100644 --- a/tests/CI_docker/context/sarvey_ci.docker +++ b/tests/CI_docker/context/sarvey_ci.docker @@ -1,5 +1,8 @@ FROM condaforge/miniforge3:latest +# Set Mamba root prefix +ENV MAMBA_ROOT_PREFIX="/opt/conda" + # update base environment RUN --mount=type=cache,target=/opt/conda/pkgs \ mamba update --all -y && \ @@ -8,22 +11,16 @@ RUN --mount=type=cache,target=/opt/conda/pkgs \ ARG DEBIAN_FRONTEND=noninteractive RUN mkdir actions-runner; cd actions-runner && \ - apt-get update && apt-get install -y curl + apt-get update && apt-get install -y curl gfortran build-essential openssh-client WORKDIR /actions-runner -RUN curl -o actions-runner-linux-x64-2.317.0.tar.gz -L https://github.com/actions/runner/releases/download/v2.317.0/actions-runner-linux-x64-2.317.0.tar.gz && \ - tar xzf ./actions-runner-linux-x64-2.317.0.tar.gz && \ +RUN curl -o actions-runner-linux-x64.tar.gz -L https://github.com/actions/runner/releases/download/v2.322.0/actions-runner-linux-x64-2.322.0.tar.gz&& \ + tar xzf ./actions-runner-linux-x64.tar.gz && \ ./bin/installdependencies.sh && \ useradd -m runneruser && \ chown -R runneruser:runneruser /actions-runner -USER runneruser - -RUN ./config.sh --url https://github.com/luhipi/sarvey --token --unattended --replace --name mefe2_sarvey_ci_1.0.0 --labels self-hosted - -USER root - # install some needed packages RUN --mount=type=cache,target=/opt/conda/pkgs \ mamba install -y bzip2 fish gcc gdb git ipython make nano pip tree wget unzip @@ -31,22 +28,29 @@ RUN --mount=type=cache,target=/opt/conda/pkgs \ # use bash shell instead of sh shell for all docker commands SHELL ["/bin/bash", "-c"] -# copy some needed stuff to /root -COPY ../environment_sarvey.yml . -# create ci_env environment +# Create ci_env environment with pip installed RUN --mount=type=cache,target=/opt/conda/pkgs \ - pip install conda-merge && \ - wget https://raw.githubusercontent.com/insarlab/MiaplPy/main/conda-env.yml && \ - conda-merge conda-env.yml environment_sarvey.yml > env.yml && \ - mamba env create -n ci_env -f env.yml && \ - source /opt/conda/bin/activate ci_env && \ - pip install git+https://github.com/insarlab/MiaplPy.git && \ - conda list && \ + conda create -n ci_env python=3.10 pip -y && \ conda clean -afy -USER runneruser +# Install additional packages using mamba and pip +RUN --mount=type=cache,target=/opt/conda/pkgs \ + conda install -n ci_env conda-forge::pysolid -y && \ + conda install -n ci_env conda-forge::gdal && \ + conda run -n ci_env pip install git+https://github.com/insarlab/MiaplPy.git && \ + # TODO: replace the followong with the main branch + #conda run -n ci_env pip install git+https://github.com/luhipi/sarvey.git@main && \ + conda run -n ci_env pip install git+https://github.com/mahmud1/sarvey.git@update-runner2 && \ + conda run -n ci_env pip install sarvey[dev] && \ + conda run -n ci_env pip install sphinx_rtd_theme && \ + conda clean -afy +COPY ../entrypoint.sh . +RUN chown runneruser:runneruser entrypoint.sh +RUN chmod +x entrypoint.sh + +USER runneruser RUN chmod +x /actions-runner/run.sh -ENTRYPOINT ["/actions-runner/run.sh"] +ENTRYPOINT ["./entrypoint.sh"] diff --git a/tests/test_processing.py b/tests/test_processing.py index 3fdfa270..4f15ffdf 100644 --- a/tests/test_processing.py +++ b/tests/test_processing.py @@ -103,8 +103,8 @@ def testUnwrappingTimeAndSpace(self): run(config=config, args=args, logger=self.logger) assert glob(os.path.join(config.general.output_path, "background_map.h5")), \ 'Processing failed (background_map.h5 not created).' - assert glob(os.path.join(config.general.output_path, "coordinates_utm.h5")), \ - 'Processing failed (coordinates_utm.h5 not created).' + assert glob(os.path.join(config.general.output_path, "coordinates_map.h5")), \ + 'Processing failed (coordinates_map.h5 not created).' assert glob(os.path.join(config.general.output_path, "ifg_network.h5")), \ 'Processing failed (ifg_network.h5 not created).' assert glob(os.path.join(config.general.output_path, "ifg_stack.h5")), \ @@ -143,18 +143,17 @@ def testUnwrappingTimeAndSpace(self): 'Processing failed (p1_ts_filt.h5 not created).' assert glob(os.path.join(config.general.output_path, "p1_aps.h5")), \ 'Processing failed (p1_aps.h5 not created).' - coh_value = int(config.filtering.coherence_p2 * 100) - assert glob(os.path.join(config.general.output_path, f"p2_coh{coh_value}_ifg_wr.h5")), \ - f'Processing failed (p2_coh{coh_value}_ifg_wr.h5 not created).' - assert glob(os.path.join(config.general.output_path, f"p2_coh{coh_value}_aps.h5")), \ - f'Processing failed (p2_coh{coh_value}_aps.h5 not created).' + assert glob(os.path.join(config.general.output_path, "aps_parameters.h5")), \ + 'Processing failed (aps_parameters.h5 not created).' # densification args.start = 4 args.stop = 4 run(config=config, args=args, logger=self.logger) - coh_value = int(config.filtering.coherence_p2 * 100) + coh_value = int(config.densification.coherence_p2 * 100) + assert glob(os.path.join(config.general.output_path, f"p2_coh{coh_value}_ifg_wr.h5")), \ + f'Processing failed (p2_coh{coh_value}_ifg_wr.h5 not created).' assert glob(os.path.join(config.general.output_path, f"p2_coh{coh_value}_ifg_unw.h5")), \ f'Processing failed (p2_coh{coh_value}_ifg_unw.h5 not created).' assert glob(os.path.join(config.general.output_path, f"p2_coh{coh_value}_ts.h5")), \ @@ -176,8 +175,8 @@ def testUnwrappingSpace(self): run(config=config, args=args, logger=self.logger) assert glob(os.path.join(config.general.output_path, "background_map.h5")), \ 'Processing failed (background_map.h5 not created).' - assert glob(os.path.join(config.general.output_path, "coordinates_utm.h5")), \ - 'Processing failed (coordinates_utm.h5 not created).' + assert glob(os.path.join(config.general.output_path, "coordinates_map.h5")), \ + 'Processing failed (coordinates_map.h5 not created).' assert glob(os.path.join(config.general.output_path, "ifg_network.h5")), \ 'Processing failed (ifg_network.h5 not created).' assert glob(os.path.join(config.general.output_path, "ifg_stack.h5")), \ @@ -216,18 +215,17 @@ def testUnwrappingSpace(self): 'Processing failed (p1_ts_filt.h5 not created).' assert glob(os.path.join(config.general.output_path, "p1_aps.h5")), \ 'Processing failed (p1_aps.h5 not created).' - coh_value = int(config.filtering.coherence_p2 * 100) - assert glob(os.path.join(config.general.output_path, f"p2_coh{coh_value}_ifg_wr.h5")), \ - f'Processing failed (p2_coh{coh_value}_ifg_wr.h5 not created).' - assert glob(os.path.join(config.general.output_path, f"p2_coh{coh_value}_aps.h5")), \ - f'Processing failed (p2_coh{coh_value}_aps.h5 not created).' + assert glob(os.path.join(config.general.output_path, "aps_parameters.h5")), \ + 'Processing failed (aps_parameters.h5 not created).' # densification args.start = 4 args.stop = 4 run(config=config, args=args, logger=self.logger) - coh_value = int(config.filtering.coherence_p2 * 100) + coh_value = int(config.densification.coherence_p2 * 100) + assert glob(os.path.join(config.general.output_path, f"p2_coh{coh_value}_ifg_wr.h5")), \ + f'Processing failed (p2_coh{coh_value}_ifg_wr.h5 not created).' assert glob(os.path.join(config.general.output_path, f"p2_coh{coh_value}_ifg_unw.h5")), \ f'Processing failed (p2_coh{coh_value}_ifg_unw.h5 not created).' assert glob(os.path.join(config.general.output_path, f"p2_coh{coh_value}_ts.h5")), \ @@ -242,7 +240,7 @@ def testPhaseLinking(self): config.phase_linking.use_phase_linking_results = True config.phase_linking.use_ps = True config.general.apply_temporal_unwrapping = True - config.filtering.coherence_p2 = 0.75 + config.densification.coherence_p2 = 0.75 # preparation args.start = 0 @@ -251,8 +249,8 @@ def testPhaseLinking(self): run(config=config, args=args, logger=self.logger) assert glob(os.path.join(config.general.output_path, "background_map.h5")), \ 'Processing failed (background_map.h5 not created).' - assert glob(os.path.join(config.general.output_path, "coordinates_utm.h5")), \ - 'Processing failed (coordinates_utm.h5 not created).' + assert glob(os.path.join(config.general.output_path, "coordinates_map.h5")), \ + 'Processing failed (coordinates_map.h5 not created).' assert glob(os.path.join(config.general.output_path, "ifg_network.h5")), \ 'Processing failed (ifg_network.h5 not created).' assert glob(os.path.join(config.general.output_path, "ifg_stack.h5")), \ @@ -291,18 +289,17 @@ def testPhaseLinking(self): 'Processing failed (p1_ts_filt.h5 not created).' assert glob(os.path.join(config.general.output_path, "p1_aps.h5")), \ 'Processing failed (p1_aps.h5 not created).' - coh_value = int(config.filtering.coherence_p2 * 100) - assert glob(os.path.join(config.general.output_path, f"p2_coh{coh_value}_ifg_wr.h5")), \ - f'Processing failed (p2_coh{coh_value}_ifg_wr.h5 not created).' - assert glob(os.path.join(config.general.output_path, f"p2_coh{coh_value}_aps.h5")), \ - f'Processing failed (p2_coh{coh_value}_aps.h5 not created).' + assert glob(os.path.join(config.general.output_path, "aps_parameters.h5")), \ + 'Processing failed (aps_parameters.h5 not created).' # densification args.start = 4 args.stop = 4 run(config=config, args=args, logger=self.logger) - coh_value = int(config.filtering.coherence_p2 * 100) + coh_value = int(config.densification.coherence_p2 * 100) + assert glob(os.path.join(config.general.output_path, f"p2_coh{coh_value}_ifg_wr.h5")), \ + f'Processing failed (p2_coh{coh_value}_ifg_wr.h5 not created).' assert glob(os.path.join(config.general.output_path, f"p2_coh{coh_value}_ifg_unw.h5")), \ f'Processing failed (p2_coh{coh_value}_ifg_unw.h5 not created).' assert glob(os.path.join(config.general.output_path, f"p2_coh{coh_value}_ts.h5")), \ @@ -317,7 +314,7 @@ def testMasking(self): config.general.apply_temporal_unwrapping = False config.preparation.ifg_network_type = "sb" config.consistency_check.mask_p1_file = "tests/testdata/aoi_mask.h5" - config.filtering.mask_p2_file = "tests/testdata/aoi_mask.h5" + config.densification.mask_p2_file = "tests/testdata/aoi_mask.h5" # preparation args.start = 0 @@ -326,8 +323,8 @@ def testMasking(self): run(config=config, args=args, logger=self.logger) assert glob(os.path.join(config.general.output_path, "background_map.h5")), \ 'Processing failed (background_map.h5 not created).' - assert glob(os.path.join(config.general.output_path, "coordinates_utm.h5")), \ - 'Processing failed (coordinates_utm.h5 not created).' + assert glob(os.path.join(config.general.output_path, "coordinates_map.h5")), \ + 'Processing failed (coordinates_map.h5 not created).' assert glob(os.path.join(config.general.output_path, "ifg_network.h5")), \ 'Processing failed (ifg_network.h5 not created).' assert glob(os.path.join(config.general.output_path, "ifg_stack.h5")), \ @@ -366,18 +363,17 @@ def testMasking(self): 'Processing failed (p1_ts_filt.h5 not created).' assert glob(os.path.join(config.general.output_path, "p1_aps.h5")), \ 'Processing failed (p1_aps.h5 not created).' - coh_value = int(config.filtering.coherence_p2 * 100) - assert glob(os.path.join(config.general.output_path, f"p2_coh{coh_value}_ifg_wr.h5")), \ - f'Processing failed (p2_coh{coh_value}_ifg_wr.h5 not created).' - assert glob(os.path.join(config.general.output_path, f"p2_coh{coh_value}_aps.h5")), \ - f'Processing failed (p2_coh{coh_value}_aps.h5 not created).' + assert glob(os.path.join(config.general.output_path, "aps_parameters.h5")), \ + 'Processing failed (aps_parameters.h5 not created).' # densification args.start = 4 args.stop = 4 run(config=config, args=args, logger=self.logger) - coh_value = int(config.filtering.coherence_p2 * 100) + coh_value = int(config.densification.coherence_p2 * 100) + assert glob(os.path.join(config.general.output_path, f"p2_coh{coh_value}_ifg_wr.h5")), \ + f'Processing failed (p2_coh{coh_value}_ifg_wr.h5 not created).' assert glob(os.path.join(config.general.output_path, f"p2_coh{coh_value}_ifg_unw.h5")), \ f'Processing failed (p2_coh{coh_value}_ifg_unw.h5 not created).' assert glob(os.path.join(config.general.output_path, f"p2_coh{coh_value}_ts.h5")), \ diff --git a/tests/testdata/config_test.json b/tests/testdata/config_test.json index 6cfa1abd..4a5b916a 100644 --- a/tests/testdata/config_test.json +++ b/tests/testdata/config_test.json @@ -41,16 +41,16 @@ "use_arcs_from_temporal_unwrapping": true }, "filtering": { - "coherence_p2": 0.9, "apply_aps_filtering": true, "interpolation_method": "kriging", "grid_size": 1000, - "mask_p2_file": null, "use_moving_points": true, "max_temporal_autocorrelation": 0.3 }, "densification": { + "coherence_p2": 0.9, "num_connections_to_p1": 5, + "mask_p2_file": null, "max_distance_to_p1": 2000, "velocity_bound": 0.15, "dem_error_bound": 100.0,