diff --git a/.circleci/config.yml b/.circleci/config.yml index 21ee62ba7..4e11186c6 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -1,63 +1,91 @@ version: 2.1 + jobs: build: + parameters: + python-version: + type: string docker: - image: cimg/base:current user: root steps: - checkout + - run: - name: Setup micromamba + name: Install micromamba + shell: /bin/bash -l command: | apt update --yes && apt-get upgrade --yes apt install -y --no-install-recommends wget ca-certificates git - cd $HOME + + cd ${HOME} curl -Ls https://micro.mamba.pm/api/micromamba/linux-64/latest | tar -xvj bin/micromamba + eval "$(${HOME}/bin/micromamba shell hook -s posix)" - run: - name: Setup environment + name: Create RAiDER environment + shell: /bin/bash -l command: | - eval "$($HOME/bin/micromamba shell hook -s posix)" + eval "$(${HOME}/bin/micromamba shell hook -s posix)" + + PYTHON_VERSION="<< parameters.python-version >>" + sed -i "/python>=/c\ - python=${PYTHON_VERSION}" environment.yml micromamba create -f environment.yml + + - run: + name: Install raider and check environment + shell: /bin/bash -l + command: | + eval "$(${HOME}/bin/micromamba shell hook -s posix)" micromamba activate RAiDER - pip install coveralls - echo url: https://cds.climate.copernicus.eu/api/v2 > $HOME/.cdsapirc - echo key: $cdsak >> $HOME/.cdsapirc - - echo { > $HOME/.ecmwfapirc - echo ' "url": "https://api.ecmwf.int/v1",' >> $HOME/.ecmwfapirc - echo ' "email": "'$ecmwfu'",' >> $HOME/.ecmwfapirc - echo ' "key": "'$ecmwfk'"' >> $HOME/.ecmwfapirc - echo } >> $HOME/.ecmwfapirc - - echo url: $NCUMloc > $HOME/.ncmrlogin - echo username: $NCUMu >> $HOME/.ncmrlogin - echo password: $NCUMp >> $HOME/.ncmrlogin + + python -m pip install --no-deps . + + python -c "import RAiDER; from RAiDER.delay import tropo_delay" + python -c "import RAiDER; from RAiDER.interpolator import interp_along_axis" python --version python -c "import numpy; print(numpy.__version__)" python -c "import pyproj; print(pyproj.__version__)" + - run: - name: Install RAiDER and test the install + name: Setup data stores + shell: /bin/bash -l command: | - eval "$($HOME/bin/micromamba shell hook -s posix)" + eval "$(${HOME}/bin/micromamba shell hook -s posix)" micromamba activate RAiDER - python -m pip install . - python -c "import RAiDER; from RAiDER.delay import tropo_delay" - python -c "import RAiDER; from RAiDER.interpolator import interp_along_axis" + + python -c 'from RAiDER.models.credentials import setup_from_env; setup_from_env()' + - run: name: Run unit tests shell: /bin/bash -l command: | - eval "$($HOME/bin/micromamba shell hook -s posix)" + eval "$(${HOME}/bin/micromamba shell hook -s posix)" micromamba activate RAiDER + COV_OPTIONS=`python -c "import importlib;print(*(' --cov='+p for p in importlib.util.find_spec('RAiDER').submodule_search_locations))"` - pytest -m "not long" test/ $COV_OPTIONS --cov-report= + pytest -m "not long" test/ ${COV_OPTIONS} --cov-report= + - run: name: Report coverage + shell: /bin/bash -l command: | - eval "$($HOME/bin/micromamba shell hook -s posix)" - micromamba activate RAiDER - python .circleci/fix_coverage_paths.py .coverage $(pwd)/tools/RAiDER/ - coverage report -mi - coveralls + PYTHON_VERSION="<< parameters.python-version >>" + if [ "${PYTHON_VERSION}" == "3.12" ]; then + eval "$(${HOME}/bin/micromamba shell hook -s posix)" + micromamba activate RAiDER + + python -m pip install coveralls + python .circleci/fix_coverage_paths.py .coverage ${PWD}/tools/RAiDER/ + coverage report -mi + coveralls + fi + +workflows: + all-tests: + jobs: + - build: + matrix: + parameters: + python-version: ["3.9", "3.10", "3.11", "3.12"] diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 671db476e..fc4d13e32 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -12,13 +12,13 @@ on: jobs: call-version-info-workflow: - uses: ASFHyP3/actions/.github/workflows/reusable-version-info.yml@v0.8.2 + uses: ASFHyP3/actions/.github/workflows/reusable-version-info.yml@v0.8.3 with: python_version: '3.10' call-docker-ghcr-workflow: needs: call-version-info-workflow - uses: ASFHyP3/actions/.github/workflows/reusable-docker-ghcr.yml@v0.8.2 + uses: ASFHyP3/actions/.github/workflows/reusable-docker-ghcr.yml@v0.8.3 with: version_tag: ${{ needs.call-version-info-workflow.outputs.version_tag }} release_branch: main diff --git a/.github/workflows/changelog.yml b/.github/workflows/changelog.yml index 0a50ce126..369b918a5 100644 --- a/.github/workflows/changelog.yml +++ b/.github/workflows/changelog.yml @@ -13,6 +13,6 @@ on: jobs: call-changelog-check-workflow: - uses: ASFHyP3/actions/.github/workflows/reusable-changelog-check.yml@v0.8.2 + uses: ASFHyP3/actions/.github/workflows/reusable-changelog-check.yml@v0.8.3 secrets: USER_TOKEN: ${{ secrets.GITHUB_TOKEN }} diff --git a/.github/workflows/labeled-pr.yml b/.github/workflows/labeled-pr.yml index ee8a5f5ea..103ae29a7 100644 --- a/.github/workflows/labeled-pr.yml +++ b/.github/workflows/labeled-pr.yml @@ -12,4 +12,4 @@ on: jobs: call-labeled-pr-check-workflow: - uses: ASFHyP3/actions/.github/workflows/reusable-labeled-pr-check.yml@v0.8.2 + uses: ASFHyP3/actions/.github/workflows/reusable-labeled-pr-check.yml@v0.8.3 diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index bb7e27ce1..310ea6db4 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -7,7 +7,7 @@ on: jobs: call-release-workflow: - uses: ASFHyP3/actions/.github/workflows/reusable-release.yml@v0.8.2 + uses: ASFHyP3/actions/.github/workflows/reusable-release.yml@v0.8.3 with: release_prefix: RAiDER develop_branch: dev diff --git a/.github/workflows/tag.yml b/.github/workflows/tag.yml index ad0761eb0..f2328e9c7 100644 --- a/.github/workflows/tag.yml +++ b/.github/workflows/tag.yml @@ -7,7 +7,7 @@ on: jobs: call-bump-version-workflow: - uses: ASFHyP3/actions/.github/workflows/reusable-bump-version.yml@v0.8.2 + uses: ASFHyP3/actions/.github/workflows/reusable-bump-version.yml@v0.8.3 with: user: dbekaert email: bekaertdavid@gmail.com diff --git a/CHANGELOG.md b/CHANGELOG.md index e70ad1834..598687cb8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,18 +6,43 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [PEP 440](https://www.python.org/dev/peps/pep-0440/) and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html). + +## [0.4.6] + +### Added +* Adds an `s1_orbits.py` module which includes: + * `get_orbits_from_slc_ids` to download the associated orbit files for a list of Sentinel-1 SLC IDs + * `ensure_orbit_credentials` to ensure ESA CSDE credentials have been provides to download orbit files. This should be called before `sentineleof` is used to download orbits. +* Adds a `setup_from_env` function to `models/credentials.py` which will pull *all* credentials needed for acquiring weather model data from environment variables and ensure the correct config file is written. This makes setting up credentials in CI pipelines significantly easier + +### Changed +* `sentineleof` upgraded to version 0.9.5 or later to (a) fetch orbits from ESA CDSE and (b) ensure that if CDSE fetch fails, code resorts to ASF orbit repository + +### Fixes +* RAiDER is now tested on Python version 3.9-3.12 +* All typehints are now Python 3.9 compatible +* [607](https://github.com/dbekaert/RAiDER/issues/607): Python entrypoint loading is now compatible with Python 3.12 +* [610](https://github.com/dbekaert/RAiDER/issues/610): Sentinel-1 orbit availability due to ESA migrating Sentinel-1 orbit files from Copernicus Open Access Hub (Scihub) to the new Copernicus Data Space Ecosystem (CDSE) +* make weather file directory when it doesn't exist +* Ensures the `models/data/alaska.geojson.zip` file is packaged when building from the source tarball +* Make ISCE3 an optional dependency in `s1_azimuth_timing.py` ++ Added unit tests and removed unused and depracated functions + +### Removed +* `hyp3lib`, which was only used for downloading orbit fies, has been removed in favor of `sentineleof` + ## [0.4.5] -## Fixes +### Fixes * [#583](https://github.com/dbekaert/RAiDER/issues/583): it appears that since the issues with geo2rdr cropped up during our processing campaign, there has been a new release of ISCE3 that resolves these failures with `geo2rdr` and the time interpolation that uses this ISCE3 routine. * [#584](https://github.com/dbekaert/RAiDER/issues/584): failed Raider step function in hyp3 job submission when HRRR model times are not available (even within the valid model range) - to resolve, we check availability of files when delay workflow called with a) azimuth_grid_interpolation and b) input to workflow is GUNW. If weather model files are unavailable and the GUNW is on s3, do nothing to GUNW (i.e. do not add tropo delay) and exit successfully. If weather model files are unavailable and the GUNW is on local disk, raise `ValueError` * [#587](https://github.com/dbekaert/RAiDER/issues/587): similar to 584 except added here to the mix is control flow in RAiDER.py passes over numerous exceptions in workflow. This is fixed identically as above. * [#596](https://github.com/dbekaert/RAiDER/issues/596): the "prefix" for aws does not include the final netcdf file name, just the sub-directories in the bucket and therefore extra logic must be added to determine the GUNW netcdf file name (and the assocaited reference/secondary dates). We proceed by downloading the data which is needed regardless. Tests are updated. -## Removed +### Removed * Removes `update` option (either `True` or `False`) from calcGUNW workflow which asks whether the GUNW should be updated or not. In existing code, it was not being used/applied, i.e. previous workflow always updated GUNW. Removed input arguments related from respective functions so that it can be updated later. -## Added +### Added * Allow for Hyp3 GUNW workflow with HRRR (i.e. specifying a gunw path in s3) to successfully exit if any of the HRRR model times required for `azimuth-time-grid` interpolation (which is default interpolatin method) are not available when using bucket inputs (i.e. only on the cloud) * For GUNW workflow, when model is HRRR, azimuth_time_grid interpolation used, and using a local GUNW, if requisite weather model files are not available for raise a ValueError (before processing) * Raise a value error if non-unique dates are given in the function `get_inverse_weights_for_dates` in `s1_azimuth_timing.py` @@ -30,7 +55,7 @@ and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html). * Ensures ISCE3 is `>=0.15.0` * Uses correct HyP3 S3 prefix conventions and filename suffix within test patches to improve readability of what tests are mocking (see comments in #597). -## Changed +### Changed * Get only 2 or 3 model times required for azimuth-time-interpolation (previously obtained all 3 as it was easier to implement) - this ensures slightly less failures associated with HRRR availability. Importantly, if a acquisition time occurs during a model time, then we order by distance to the reference time and how early it occurs (so earlier times come first if two times are equidistant to the aquisition time). * Made test names in `test_GUNW.py` more descriptive * Numpy docstrings and general linting to modified function including removing variables that were not being accessed @@ -40,11 +65,10 @@ and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html). * Fixes tests for checking availability of HRRR due Issue #596 (above). ## [0.4.4] - -## Fixes * For s1-azimuth-time interpolation, overlapping orbits when one orbit does not cover entire GUNW product errors out. We now ensure state-vectors are both unique and in order before creating a orbit object in ISCE3. ## [0.4.3] ++ Bug fixes, unit tests, docstrings + Prevent ray tracing integration from occuring at exactly top of weather model + Properly expose z_ref (max integration height) parameter, and dont allow higher than weather model + Min version for sentineleof for obtaining restituted orbits. diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index b29736270..eb2173730 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -13,6 +13,34 @@ If you get stuck at any point you can create an [issue on GitHub](https://github For more information on contributing to open source projects, [GitHub's own guide](https://guides.github.com/activities/contributing-to-open-source/) is a great starting point if you are new to version control. +## Optional Dependencies + +In order to better support the NISAR SDS (see: [#533](https://github.com/dbekaert/RAiDER/issues/533)), RAiDER has some optional dependencies: + +* ISCE3 +* Pandas +* Rasterio +* Progressbar + +RAiDER distributes two conda packages, `raider-base` a lighter-weight package that does depend on the optional dependencies, and `raider` which includes all dependencies. When using, or adding new, optional dependenices in RAiDER, please follow this pattern: +1. When you import the optional dependency, handle import errors like: + ```python + try: + import optional_dependency + except ImportError: + optional_dependency = None + ``` + Note: you *do not* need to delay imports until use with this pattern. +2. At the top of any function/method that uses the optional dependency, throw if it's missing like: + ```python + if optional_dependency is None: + raise ImportError('optional_dependency is required for this function. Use conda to install optional_dependency') + ``` +3. If you want to add type hints for objects in the optional_dependency, use a forward declaration like: + ```python + def function_that_uses_optional_dependency(obj: 'optional_dependency.obj'): + ``` + Note: the typehint is a string here. ## Git workflows ## diff --git a/environment.yml b/environment.yml index 464a9ee91..e7de60279 100644 --- a/environment.yml +++ b/environment.yml @@ -13,18 +13,15 @@ dependencies: # For running - asf_search - boto3 + - cartopy - cdsapi - cfgrib - - cmake - - cxx-compiler - - cython - dask - dem_stitcher>=2.3.1 - ecmwf-api-client + - h5netcdf - h5py - herbie-data - - h5netcdf - - hyp3lib - isce3>=0.15.0 - jsonschema==3.2.0 # this is for ASF DAAC ingest schema validation - lxml @@ -33,36 +30,39 @@ dependencies: - numpy - pandas - progressbar - - pybind11 - pydap>3.2.2 - pyproj>=2.2.0 - pyyaml - rasterio>=1.3.0 - - rioxarray - requests + - rioxarray - s3fs - scipy>1.10.0 + - sentineleof>=0.9.5 - shapely - - sysroot_linux-64 - tqdm - xarray # For packaging and testing - autopep8 + - cmake + - cxx-compiler + - cython + - pybind11 - pytest + - pytest-console-scripts - pytest-cov - - pytest-timeout - pytest-mock - - pytest-console-scripts - - setuptools_scm >=6.2 + - pytest-timeout + - setuptools_scm>=6.2 + - sysroot_linux-64 # For docs website - mkdocs - - mkdocstrings - - mkdocstrings-python - mkdocs-macros-plugin - mkdocs-material - mkdocs-material-extensions - - sentineleof>=0.8.1 + - mkdocstrings + - mkdocstrings-python # For RAiDER-docs - - jupyterlab - jupyter_contrib_nbextensions + - jupyterlab - wand diff --git a/pyproject.toml b/pyproject.toml index f01a75a1c..d5815cdb2 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -58,7 +58,7 @@ zip-safe = false where = ["tools"] [tool.setuptools.package-data] -"*" = ["*.yml", "*.yaml"] +"*" = ["*.yml", "*.yaml", "*.zip"] [tool.isort] known_first_party = "RAiDER" diff --git a/setup.py b/setup.py index 9f163fc77..a8a76bc1e 100755 --- a/setup.py +++ b/setup.py @@ -44,5 +44,4 @@ setup( ext_modules=cython_extensions + pybind_extensions, cmdclass={"build_ext": build_ext}, - package_data={'tools': ['RAiDER/models/*.zip']} ) diff --git a/test/scenario_0/small_dem.tif b/test/scenario_0/small_dem.tif new file mode 100644 index 000000000..bab74e11a Binary files /dev/null and b/test/scenario_0/small_dem.tif differ diff --git a/test/test_GUNW.py b/test/test_GUNW.py index 71c18db90..a57152d2f 100644 --- a/test/test_GUNW.py +++ b/test/test_GUNW.py @@ -2,12 +2,10 @@ import json import os import shutil -import subprocess import unittest from pathlib import Path import eof.download -import hyp3lib import jsonschema import numpy as np import pandas as pd @@ -17,12 +15,14 @@ import RAiDER import RAiDER.cli.raider as raider +import RAiDER.s1_azimuth_timing from RAiDER import aws from RAiDER.aria.prepFromGUNW import ( check_hrrr_dataset_availablity_for_s1_azimuth_time_interpolation, - check_weather_model_availability + check_weather_model_availability, ) from RAiDER.cli.raider import calcDelaysGUNW +from RAiDER.models.customExceptions import * def compute_transform(lats, lons): @@ -206,24 +206,13 @@ def test_azimuth_timing_interp_against_center_time_interp(weather_model_name: st # https://github.com/dbekaert/RAiDER/blob/ # f77af9ce2d3875b00730603305c0e92d6c83adc2/tools/RAiDER/aria/prepFromGUNW.py#L151-L200 - mocker.patch('hyp3lib.get_orb.downloadSentinelOrbitFile', - # Hyp3 Lib returns 2 values - side_effect=[ - # For azimuth time - (orbit_dict_for_azimuth_time_test['reference'], ''), - (orbit_dict_for_azimuth_time_test['secondary'], ''), - (orbit_dict_for_azimuth_time_test['secondary'], ''), - ] - ) - # For prepGUNW side_effect = [ - # center-time - orbit_dict_for_azimuth_time_test['reference'], - # azimuth-time - orbit_dict_for_azimuth_time_test['reference'], - ] - side_effect = list(map(str, side_effect)) + # center-time + [Path(orbit_dict_for_azimuth_time_test['reference'])], + # azimuth-time + [Path(orbit_dict_for_azimuth_time_test['reference'])], + ] mocker.patch('eof.download.download_eofs', side_effect=side_effect) @@ -237,6 +226,15 @@ def test_azimuth_timing_interp_against_center_time_interp(weather_model_name: st ['secondary_slc_id', 'secondary_slc_id'], ]) + mocker.patch( + 'RAiDER.s1_azimuth_timing.get_orbits_from_slc_ids', + side_effect=[ + # For azimuth time + [Path(orbit_dict_for_azimuth_time_test['reference'])], + [Path(orbit_dict_for_azimuth_time_test['secondary']), Path(orbit_dict_for_azimuth_time_test['secondary'])], + ] + ) + side_effect = (weather_model_dict_for_center_time_test[weather_model_name] + weather_model_dict_for_azimuth_time_test[weather_model_name]) # RAiDER needs strings for paths @@ -259,9 +257,8 @@ def test_azimuth_timing_interp_against_center_time_interp(weather_model_name: st # Calls 4 times for azimuth time and 4 times for center time assert RAiDER.processWM.prepareWeatherModel.call_count == 8 - # Only calls for azimuth timing for reference and secondary; - # There is 1 slc for ref and 2 for secondary, totalling 3 calls - assert hyp3lib.get_orb.downloadSentinelOrbitFile.call_count == 3 + # Only calls once each ref and sec list of slcs + assert RAiDER.s1_azimuth_timing.get_orbits_from_slc_ids.call_count == 2 # Only calls for azimuth timing: once for ref and sec assert RAiDER.s1_azimuth_timing.get_slc_id_from_point_and_time.call_count == 2 # Once for center-time and azimuth-time each @@ -371,24 +368,13 @@ def test_provenance_metadata_for_tropo_group(weather_model_name: str, out_path = shutil.copy(gunw_azimuth_test, tmp_path / out) if interp_method == 'azimuth_time_grid': - mocker.patch('hyp3lib.get_orb.downloadSentinelOrbitFile', - # Hyp3 Lib returns 2 values - side_effect=[ - # For azimuth time - (orbit_dict_for_azimuth_time_test['reference'], ''), - (orbit_dict_for_azimuth_time_test['secondary'], ''), - (orbit_dict_for_azimuth_time_test['secondary'], ''), - ] - ) - # For prepGUNW side_effect = [ - # center-time - orbit_dict_for_azimuth_time_test['reference'], - # azimuth-time - orbit_dict_for_azimuth_time_test['reference'], - ] - side_effect = list(map(str, side_effect)) + # center-time + [Path(orbit_dict_for_azimuth_time_test['reference'])], + # azimuth-time + [Path(orbit_dict_for_azimuth_time_test['reference'])], + ] mocker.patch('eof.download.download_eofs', side_effect=side_effect) @@ -402,6 +388,14 @@ def test_provenance_metadata_for_tropo_group(weather_model_name: str, ['secondary_slc_id', 'secondary_slc_id'], ]) + mocker.patch( + 'RAiDER.s1_azimuth_timing.get_orbits_from_slc_ids', + side_effect=[ + # For azimuth time + [Path(orbit_dict_for_azimuth_time_test['reference'])], + [Path(orbit_dict_for_azimuth_time_test['secondary']), Path(orbit_dict_for_azimuth_time_test['secondary'])], + ] + ) weather_model_path_dict = (weather_model_dict_for_center_time_test if interp_method == 'center_time' else weather_model_dict_for_azimuth_time_test) @@ -475,24 +469,14 @@ def test_GUNW_workflow_fails_if_a_download_fails(gunw_azimuth_test, orbit_dict_f # The first part is the same mock up as done in test_azimuth_timing_interp_against_center_time_interp # Maybe better mocks could be done - but this is sufficient or simply a factory for this test given # This is reused so many times. - mocker.patch('hyp3lib.get_orb.downloadSentinelOrbitFile', - # Hyp3 Lib returns 2 values - side_effect=[ - # For azimuth time - (orbit_dict_for_azimuth_time_test['reference'], ''), - (orbit_dict_for_azimuth_time_test['secondary'], ''), - (orbit_dict_for_azimuth_time_test['secondary'], ''), - ] - ) # For prepGUNW side_effect = [ - # center-time - orbit_dict_for_azimuth_time_test['reference'], - # azimuth-time - orbit_dict_for_azimuth_time_test['reference'], - ] - side_effect = list(map(str, side_effect)) + # center-time + [Path(orbit_dict_for_azimuth_time_test['reference'])], + # azimuth-time + [Path(orbit_dict_for_azimuth_time_test['reference'])], + ] mocker.patch('eof.download.download_eofs', side_effect=side_effect) @@ -506,6 +490,15 @@ def test_GUNW_workflow_fails_if_a_download_fails(gunw_azimuth_test, orbit_dict_f ['secondary_slc_id', 'secondary_slc_id'], ]) + mocker.patch( + 'RAiDER.s1_azimuth_timing.get_orbits_from_slc_ids', + side_effect=[ + # For azimuth time + [Path(orbit_dict_for_azimuth_time_test['reference'])], + [Path(orbit_dict_for_azimuth_time_test['secondary']), Path(orbit_dict_for_azimuth_time_test['secondary'])], + ] + ) + # These are the important parts of this test # Makes sure that a value error is raised if a download fails via a Runtime Error # There are two weather model files required for this particular mock up. First, one fails. @@ -517,7 +510,7 @@ def test_GUNW_workflow_fails_if_a_download_fails(gunw_azimuth_test, orbit_dict_f '-interp', 'azimuth_time_grid' ] - with pytest.raises(ValueError): + with pytest.raises(RuntimeError): calcDelaysGUNW(iargs_1) RAiDER.s1_azimuth_timing.get_s1_azimuth_time_grid.assert_not_called() @@ -536,6 +529,6 @@ def test_value_error_for_file_inputs_when_no_data_available(mocker): '-interp', 'azimuth_time_grid' ] - with pytest.raises(ValueError): + with pytest.raises(NoWeatherModelData): calcDelaysGUNW(iargs) RAiDER.aria.prepFromGUNW.main.assert_not_called() diff --git a/test/test_checkArgs.py b/test/test_checkArgs.py index 5e4604c98..9f210ee80 100644 --- a/test/test_checkArgs.py +++ b/test/test_checkArgs.py @@ -10,11 +10,12 @@ from test import TEST_DIR, pushd from RAiDER.cli import DEFAULT_DICT -from RAiDER.checkArgs import checkArgs, makeDelayFileNames +from RAiDER.checkArgs import checkArgs, makeDelayFileNames, get_raster_ext from RAiDER.llreader import BoundingBox, StationFile, RasterRDR from RAiDER.losreader import Zenith, Conventional, Raytracing from RAiDER.models.gmao import GMAO + SCENARIO_1 = os.path.join(TEST_DIR, "scenario_1") SCENARIO_2 = os.path.join(TEST_DIR, "scenario_2") @@ -180,4 +181,7 @@ def test_makeDelayFileNames_4(): ) +def test_get_raster_ext(): + with pytest.raises(ValueError): + get_raster_ext('dummy_format') diff --git a/test/test_credentials.py b/test/test_credentials.py index fa9fd88ec..caf0c7a60 100644 --- a/test/test_credentials.py +++ b/test/test_credentials.py @@ -5,26 +5,26 @@ # Test checking/creating ECMWF_RC CAPI file def test_ecmwfApi_createFile(): import ecmwfapi - + #Check extension for hidden files hidden_ext = '_' if system()=="Windows" else '.' # Test creation of ~/.ecmwfapirc file - ecmwf_file = os.path.expanduser('./') + hidden_ext + credentials.API_FILENAME['ERAI'] - credentials.check_api('ERAI', 'dummy', 'dummy', './', update_flag=True) - assert os.path.exists(ecmwf_file) == True,f'{ecmwf_file} does not exist' + ecmwf_file = os.path.expanduser('./') + hidden_ext + credentials.API_FILENAME['HRES'] + credentials.check_api('HRES', 'dummy', 'dummy', './', update_flag=True) + assert os.path.exists(ecmwf_file) == True,f'{ecmwf_file} does not exist' # Get existing ECMWF_API_RC env if exist default_ecmwf_file = os.getenv("ECMWF_API_RC_FILE") if default_ecmwf_file is None: default_ecmwf_file = ecmwfapi.api.DEFAULT_RCFILE_PATH - #Set it to current dir to avoid overwriting ~/.ecmwfapirc file + #Set it to current dir to avoid overwriting ~/.ecmwfapirc file os.environ["ECMWF_API_RC_FILE"] = ecmwf_file key, url, uid = ecmwfapi.api.get_apikey_values() - + # Return to default_ecmwf_file and remove local API file - os.environ["ECMWF_API_RC_FILE"] = default_ecmwf_file + os.environ["ECMWF_API_RC_FILE"] = default_ecmwf_file os.remove(ecmwf_file) #Check if API is written correctly @@ -33,17 +33,17 @@ def test_ecmwfApi_createFile(): # Test checking/creating Copernicus Climate Data Store -# CDS_RC CAPI file +# CDS_RC CAPI file def test_cdsApi_createFile(): import cdsapi #Check extension for hidden files hidden_ext = '_' if system()=="Windows" else '.' - + # Test creation of .cdsapirc file in current dir cds_file = os.path.expanduser('./') + hidden_ext + credentials.API_FILENAME['ERA5'] credentials.check_api('ERA5', 'dummy', 'dummy', './', update_flag=True) - assert os.path.exists(cds_file) == True,f'{cds_file} does not exist' + assert os.path.exists(cds_file) == True,f'{cds_file} does not exist' # Check the content cds_credentials = cdsapi.api.read_config(cds_file) @@ -51,21 +51,21 @@ def test_cdsApi_createFile(): # Remove local API file os.remove(cds_file) - + assert uid == 'dummy', f'{cds_file}: UID was not written correctly' assert key == 'dummy', f'{cds_file}: KEY was not written correctly' -# Test checking/creating EARTHDATA_RC API file +# Test checking/creating EARTHDATA_RC API file def test_netrcApi_createFile(): import netrc - + #Check extension for hidden files hidden_ext = '_' if system()=="Windows" else '.' # Test creation of ~/.cdsapirc file netrc_file = os.path.expanduser('./') + hidden_ext + credentials.API_FILENAME['GMAO'] credentials.check_api('GMAO', 'dummy', 'dummy', './', update_flag=True) - assert os.path.exists(netrc_file) == True,f'{netrc_file} does not exist' + assert os.path.exists(netrc_file) == True,f'{netrc_file} does not exist' # Check the content host = 'urs.earthdata.nasa.gov' @@ -74,6 +74,6 @@ def test_netrcApi_createFile(): # Remove local API file os.remove(netrc_file) - + assert uid == 'dummy', f'{netrc_file}: UID was not written correctly' - assert key == 'dummy', f'{netrc_file}: KEY was not written correctly' \ No newline at end of file + assert key == 'dummy', f'{netrc_file}: KEY was not written correctly' diff --git a/test/test_dem.py b/test/test_dem.py index c005ef3e4..a008c7e56 100644 --- a/test/test_dem.py +++ b/test/test_dem.py @@ -1,7 +1,38 @@ import os import pytest -import numpy as np +from test import TEST_DIR, pushd +from RAiDER.dem import download_dem + + +def test_download_dem_1(): + SCENARIO_1 = os.path.join(TEST_DIR, "scenario_4") + hts, meta = download_dem( + demName=os.path.join(SCENARIO_1,'warpedDEM.dem'), + overwrite=False + ) + assert hts.shape == (45,226) + assert meta['crs'] is None + + +def test_download_dem_2(): + with pytest.raises(ValueError): + download_dem() + + +def test_download_dem_3(tmp_path): + with pushd(tmp_path): + fname = os.path.join(tmp_path, 'tmp_file.nc') + with pytest.raises(ValueError): + download_dem(demName=fname) + + +@pytest.mark.long +def test_download_dem_4(tmp_path): + with pushd(tmp_path): + fname = os.path.join(tmp_path, 'tmp_file.nc') + z,m = download_dem(demName=fname, overwrite=True, ll_bounds=[37.9,38.,-91.8,-91.7], writeDEM=True) + assert len(z.shape) == 2 + assert 'crs' in m.keys() -from test import pushd diff --git a/test/test_gnss.py b/test/test_gnss.py index 3b9a621ee..b626a84a7 100644 --- a/test/test_gnss.py +++ b/test/test_gnss.py @@ -4,13 +4,20 @@ import pandas as pd -from test import pushd +from test import pushd, TEST_DIR +SCENARIO2_DIR = os.path.join(TEST_DIR, "scenario_2") + from RAiDER.gnss.processDelayFiles import ( addDateTimeToFiles, getDateTime, concatDelayFiles ) +from RAiDER.gnss.downloadGNSSDelays import ( + get_stats_by_llh, get_station_list, download_tropo_delays, + filterToBBox, +) +from RAiDER.models.customExceptions import NoStationDataFoundError def file_len(fname): @@ -99,3 +106,57 @@ def test_concatDelayFiles(tmp_path, temp_file): outName=out_name ) assert file_len(out_name) == file_length + +def test_get_stats_by_llh2(): + stations = get_stats_by_llh(llhBox=[10,18,360-93,360-88]) + assert isinstance(stations, pd.DataFrame) + + +def test_get_stats_by_llh3(): + with pytest.raises(ValueError): + get_stats_by_llh(llhBox=[10,18,-93,-88]) + + + +def test_get_station_list(): + stations, output_file = get_station_list(stationFile=os.path.join(SCENARIO2_DIR, 'stations.csv'), writeStationFile=False) + assert isinstance(stations,list) + assert isinstance(output_file,pd.DataFrame) + + +def test_download_tropo_delays1(): + with pytest.raises(NotImplementedError): + download_tropo_delays(stats=['GUAT', 'SLAC', 'CRSE'], years=[2022], gps_repo='dummy_repo') + + +def test_download_tropo_delays2(): + with pytest.raises(NoStationDataFoundError): + download_tropo_delays(stats=['dummy_station'], years=[2022]) + + +def test_download_tropo_delays2(tmp_path): + stations, output_file = get_station_list(stationFile=os.path.join(SCENARIO2_DIR, 'stations.csv')) + + # spot check a couple of stations + assert 'CAPE' in stations + assert 'FGNW' in stations + assert isinstance(output_file, str) + + # try downloading the delays + download_tropo_delays(stats=stations, years=[2022], writeDir=tmp_path) + assert True + + +def test_filterByBBox1(): + _, station_data = get_station_list(stationFile=os.path.join(SCENARIO2_DIR, 'stations.csv'), writeStationFile=False) + with pytest.raises(ValueError): + filterToBBox(station_data, llhBox=[34,38,-120,-115]) + + +def test_filterByBBox2(): + _, station_data = get_station_list(stationFile=os.path.join(SCENARIO2_DIR, 'stations.csv'), writeStationFile=False) + new_data = filterToBBox(station_data, llhBox=[34,38,240,245]) + for stat in ['CAPE','MHMS','NVCO']: + assert stat not in new_data['ID'].to_list() + for stat in ['FGNW', 'JPLT', 'NVTP', 'WLHG', 'WORG']: + assert stat in new_data['ID'].to_list() diff --git a/test/test_llreader.py b/test/test_llreader.py index 564829ae9..c8b584014 100644 --- a/test/test_llreader.py +++ b/test/test_llreader.py @@ -5,6 +5,7 @@ import pandas as pd from test import GEOM_DIR, TEST_DIR +from pyproj import CRS from RAiDER.cli.raider import calcDelays @@ -14,8 +15,9 @@ bounds_from_latlon_rasters, bounds_from_csv ) -SCENARIO2_DIR = os.path.join(TEST_DIR, "scenario_2") +SCENARIO0_DIR = os.path.join(TEST_DIR, "scenario_0") SCENARIO1_DIR = os.path.join(TEST_DIR, "scenario_1", "geom") +SCENARIO2_DIR = os.path.join(TEST_DIR, "scenario_2") @pytest.fixture @@ -41,6 +43,31 @@ def test_latlon_reader_2(): RasterRDR(lat_file='doesnotexist.rdr', lon_file='doesnotexist.rdr') +def test_aoi_epsg(): + bbox = [20, 27, -115, -104] + r = BoundingBox(bbox) + r.set_output_spacing(ll_res=0.05) + test = r.get_output_spacing(4978) + assert test == 0.05 * 1e5 + + +def test_set_output_dir(): + bbox = [20, 27, -115, -104] + r = BoundingBox(bbox) + r.set_output_directory('dummy_directory') + assert r._output_directory == 'dummy_directory' + + +def test_set_xygrid(): + bbox = [20, 27, -115, -104] + crs = CRS.from_epsg(4326) + r = BoundingBox(bbox) + r.set_output_spacing(ll_res=0.1) + r.set_output_xygrid(dst_crs=4978) + r.set_output_xygrid(dst_crs=crs) + assert True + + def test_latlon_reader(): latfile = os.path.join(GEOM_DIR, 'lat.rdr') lonfile = os.path.join(GEOM_DIR, 'lon.rdr') @@ -60,6 +87,17 @@ def test_latlon_reader(): assert all([np.allclose(b, t, rtol=1e-4) for b, t in zip(query.bounds(), bounds_true)]) +def test_badllfiles(station_file): + latfile = os.path.join(GEOM_DIR, 'lat.rdr') + lonfile = os.path.join(GEOM_DIR, 'lon_dummy.rdr') + station_file = station_file + with pytest.raises(ValueError): + RasterRDR(lat_file=latfile, lon_file=lonfile) + with pytest.raises(ValueError): + RasterRDR(lat_file=latfile, lon_file=station_file) + with pytest.raises(ValueError): + RasterRDR(lat_file=station_file, lon_file=lonfile) + def test_read_bbox(): bbox = [20, 27, -115, -104] query = BoundingBox(bbox) @@ -102,3 +140,12 @@ def test_readZ_sf(station_file): aoi = StationFile(station_file) assert np.allclose(aoi.readZ(), .1) + +def test_GeocodedFile(): + aoi = GeocodedFile(os.path.join(SCENARIO0_DIR, 'small_dem.tif'), is_dem=True) + z = aoi.readZ() + x,y = aoi.readLL() + assert z.shape == (569,558) + assert x.shape == z.shape + assert True + diff --git a/test/test_losreader.py b/test/test_losreader.py index 3def26308..2a7125354 100644 --- a/test/test_losreader.py +++ b/test/test_losreader.py @@ -147,7 +147,7 @@ def test_get_sv_3(svs): def test_get_sv_4(svs): true_svs = svs filename = os.path.join(ORB_DIR, 'no_exist.txt') - with pytest.raises(FileNotFoundError): + with pytest.raises(ValueError): get_sv(filename, true_svs[0][0], pad=3*60) diff --git a/test/test_s1_orbits.py b/test/test_s1_orbits.py new file mode 100644 index 000000000..2c7c8c8e6 --- /dev/null +++ b/test/test_s1_orbits.py @@ -0,0 +1,110 @@ +import netrc +from pathlib import Path + +import eof.download +import pytest + +from RAiDER import s1_orbits + + +def test_ensure_orbit_credentials(monkeypatch): + class EmptyNetrc(): + def __init__(self, netrc_file): + self.netrc_file = netrc_file + self.hosts = {} + def __str__(self): + return str(self.hosts) + + # No .netrc, no ESA CDSE env variables + with monkeypatch.context() as mp: + mp.setattr(netrc, 'netrc', EmptyNetrc, raising=False) + mp.delenv('ESA_USERNAME', raising=False) + mp.delenv('ESA_PASSWORD', raising=False) + with pytest.raises(ValueError): + s1_orbits.ensure_orbit_credentials() + + # No .netrc, set ESA CDSE env variables + with monkeypatch.context() as mp: + mp.setattr(netrc, 'netrc', EmptyNetrc, raising=False) + mp.setenv('ESA_USERNAME', 'foo') + mp.setenv('ESA_PASSWORD', 'bar') + mp.setattr(Path, 'write_text', lambda self, write_text: write_text) + written_credentials = s1_orbits.ensure_orbit_credentials() + assert written_credentials == str({s1_orbits.ESA_CDSE_HOST: ('foo', None, 'bar')}) + + class NoCDSENetrc(): + def __init__(self, netrc_file): + self.netrc_file = netrc_file + self.hosts = {'fizz.buzz.org': ('foo', None, 'bar')} + def __str__(self): + return str(self.hosts) + + # No CDSE in .netrc, no ESA CDSE env variables + with monkeypatch.context() as mp: + mp.setattr(netrc, 'netrc', NoCDSENetrc, raising=False) + mp.delenv('ESA_USERNAME', raising=False) + mp.delenv('ESA_PASSWORD', raising=False) + with pytest.raises(ValueError): + s1_orbits.ensure_orbit_credentials() + + # No CDSE in .netrc, set ESA CDSE env variables + with monkeypatch.context() as mp: + mp.setattr(netrc, 'netrc', NoCDSENetrc, raising=False) + mp.setenv('ESA_USERNAME', 'foo') + mp.setenv('ESA_PASSWORD', 'bar') + mp.setattr(Path, 'write_text', lambda self, write_text: write_text) + written_credentials = s1_orbits.ensure_orbit_credentials() + assert written_credentials == str({'fizz.buzz.org': ('foo', None, 'bar'), s1_orbits.ESA_CDSE_HOST: ('foo', None, 'bar')}) + + class CDSENetrc(): + def __init__(self, netrc_file): + self.netrc_file = netrc_file + self.hosts = {s1_orbits.ESA_CDSE_HOST: ('foo', None, 'bar')} + def __str__(self): + return str(self.hosts) + + # cdse in .netrc, no ESA CDSE env variables + with monkeypatch.context() as mp: + mp.setattr(netrc, 'netrc', CDSENetrc, raising=False) + mp.delenv('ESA_USERNAME', raising=False) + mp.delenv('ESA_PASSWORD', raising=False) + written_credentials = s1_orbits.ensure_orbit_credentials() + assert written_credentials is None + + # cdse in .netrc, set ESA CDSE env variables + with monkeypatch.context() as mp: + mp.setattr(netrc, 'netrc', CDSENetrc, raising=False) + mp.setenv('ESA_USERNAME', 'foo') + mp.setenv('ESA_PASSWORD', 'bar') + written_credentials = s1_orbits.ensure_orbit_credentials() + assert written_credentials is None + + +def test_get_orbits_from_slc_ids(mocker): + side_effect = [ + [Path('foo.txt')], + [Path('bar.txt'), Path('fiz.txt')], + ] + mocker.patch('eof.download.download_eofs', + side_effect=side_effect) + + orbit_files = s1_orbits.get_orbits_from_slc_ids( + ['S1A_IW_SLC__1SSV_20150621T120220_20150621T120232_006471_008934_72D8'] + ) + assert orbit_files == [Path('foo.txt')] + assert eof.download.download_eofs.call_count == 1 + eof.download.download_eofs.assert_called_with( + [ '20150621T120220', '20150621T120232'], ['S1A'], save_dir=str(Path.cwd()) + ) + + orbit_files = s1_orbits.get_orbits_from_slc_ids( + ['S1B_IW_SLC__1SDV_20201115T162313_20201115T162340_024278_02E29D_5C54', + 'S1A_IW_SLC__1SDV_20201203T162353_20201203T162420_035524_042744_6D5C'] + ) + assert orbit_files == [Path('bar.txt'), Path('fiz.txt')] + assert eof.download.download_eofs.call_count == 2 + eof.download.download_eofs.assert_called_with( + ['20201115T162313', '20201115T162340', '20201203T162353', '20201203T162420'], + ['S1A', 'S1B'], + save_dir=str(Path.cwd()) + ) diff --git a/test/test_s1_time_grid.py b/test/test_s1_time_grid.py index 0b4102b75..152bad866 100644 --- a/test/test_s1_time_grid.py +++ b/test/test_s1_time_grid.py @@ -1,7 +1,7 @@ import datetime from pathlib import Path +from typing import Union -import hyp3lib import numpy as np import pandas as pd import pytest @@ -87,9 +87,13 @@ def test_s1_timing_array_wrt_slc_center_time(gunw_azimuth_test: Path, slc_start_time = get_start_time_from_slc_id(slc_ids[n // 2]).to_pydatetime() # Azimuth time grid - mocker.patch('hyp3lib.get_orb.downloadSentinelOrbitFile', - # Hyp3 Lib returns 2 values - return_value=(orbit_dict_for_azimuth_time_test[ifg_type], '')) + mocker.patch( + 'RAiDER.s1_azimuth_timing.get_orbits_from_slc_ids', + side_effect=[ + [Path(orbit_dict_for_azimuth_time_test[ifg_type])], + ] + ) + mocker.patch('RAiDER.s1_azimuth_timing._asf_query', return_value=slc_id_dict_for_azimuth_time_test[ifg_type]) time_grid = get_s1_azimuth_time_grid(lon, lat, hgt, slc_start_time) @@ -101,9 +105,7 @@ def test_s1_timing_array_wrt_slc_center_time(gunw_azimuth_test: Path, assert np.all(abs_diff < 40) assert RAiDER.s1_azimuth_timing._asf_query.call_count == 1 - # There are 2 slc ids so it downloads orbits twice - call_count = 1 if ifg_type == 'reference' else 2 - assert hyp3lib.get_orb.downloadSentinelOrbitFile.call_count == call_count + assert RAiDER.s1_azimuth_timing.get_orbits_from_slc_ids.call_count == 1 @pytest.mark.parametrize('ifg_type', ['reference', 'secondary']) @@ -132,10 +134,13 @@ def test_s1_timing_array_wrt_variance(gunw_azimuth_test: Path, slc_start_time = get_start_time_from_slc_id(slc_ids[0]).to_pydatetime() # Azimuth time grid - # Azimuth time grid - mocker.patch('hyp3lib.get_orb.downloadSentinelOrbitFile', - # Hyp3 Lib returns 2 values - return_value=(orbit_dict_for_azimuth_time_test[ifg_type], '')) + mocker.patch( + 'RAiDER.s1_azimuth_timing.get_orbits_from_slc_ids', + side_effect=[ + [Path(orbit_dict_for_azimuth_time_test[ifg_type])], + ] + ) + mocker.patch('RAiDER.s1_azimuth_timing._asf_query', return_value=slc_id_dict_for_azimuth_time_test[ifg_type]) X = get_s1_azimuth_time_grid(lon, lat, hgt, slc_start_time) @@ -146,9 +151,7 @@ def test_s1_timing_array_wrt_variance(gunw_azimuth_test: Path, assert np.all(std_hgt < 2e-3) assert RAiDER.s1_azimuth_timing._asf_query.call_count == 1 - # There are 2 slc ids so it downloads orbits twice - call_count = 1 if ifg_type == 'reference' else 2 - assert hyp3lib.get_orb.downloadSentinelOrbitFile.call_count == call_count + assert RAiDER.s1_azimuth_timing.get_orbits_from_slc_ids.call_count == 1 def test_n_closest_dts(): @@ -216,7 +219,7 @@ def test_n_closest_dts(): @pytest.mark.parametrize('input_time, temporal_window, expected_weights', zip(input_times, windows, expected_weights_list)) def test_inverse_weighting(input_time: np.datetime64, - temporal_window: int | float, + temporal_window: Union[int, float], expected_weights: list[float]): """The test is designed to determine valid inverse weighting @@ -337,18 +340,19 @@ def test_duplicate_orbits(mocker, orbit_paths_for_duplicate_orbit_xml_test): mocker.patch('RAiDER.s1_azimuth_timing.get_slc_id_from_point_and_time', side_effect=[['slc_id_0', 'slc_id_1', 'slc_id_2', 'slc_id_3']]) - # Hyp3 Lib returns 2 values - side_effect = [(o_path, '') for o_path in orbit_paths_for_duplicate_orbit_xml_test] - mocker.patch('hyp3lib.get_orb.downloadSentinelOrbitFile', - side_effect=side_effect - ) + mocker.patch( + 'RAiDER.s1_azimuth_timing.get_orbits_from_slc_ids', + side_effect=[ + [Path(o_path) for o_path in orbit_paths_for_duplicate_orbit_xml_test], + ] + ) time_grid = get_s1_azimuth_time_grid(lon, lat, hgt, t) assert time_grid.shape == (len(hgt), len(lat), len(lon)) assert RAiDER.s1_azimuth_timing.get_slc_id_from_point_and_time.call_count == 1 - assert hyp3lib.get_orb.downloadSentinelOrbitFile.call_count == 4 + assert RAiDER.s1_azimuth_timing.get_orbits_from_slc_ids.call_count == 1 def test_get_times_for_az(): diff --git a/test/test_temporal_interpolate.py b/test/test_temporal_interpolate.py index 94f8218ab..15b2041c2 100644 --- a/test/test_temporal_interpolate.py +++ b/test/test_temporal_interpolate.py @@ -10,6 +10,7 @@ wm = 'ERA5' if WM == 'ERA-5' else WM +@pytest.mark.long def test_cube_timemean(): """ Test the mean interpolation by computing cube delays at 1:30PM vs mean of 12 PM / 3PM for GMAO """ SCENARIO_DIR = os.path.join(TEST_DIR, "INTERP_TIME") @@ -70,6 +71,7 @@ def test_cube_timemean(): return +@pytest.mark.long def test_cube_weighting(): """ Test the weighting by comparing a small crop with numpy directly """ from datetime import datetime @@ -142,4 +144,3 @@ def test_cube_weighting(): os.remove('temp.yaml') return - diff --git a/test/test_util.py b/test/test_util.py index 5539fbb70..8bc086587 100644 --- a/test/test_util.py +++ b/test/test_util.py @@ -7,13 +7,15 @@ import pyproj import rasterio -from test import TEST_DIR +from test import TEST_DIR, pushd from RAiDER.utilFcns import ( _least_nonzero, cosd, rio_open, sind, writeArrayToRaster, rio_profile, rio_extents, getTimeFromFile, enu2ecef, ecef2enu, transform_bbox, clip_bbox, get_nearest_wmtimes, + robmax,robmin,padLower,convertLons, + projectDelays,floorish, ) @@ -229,11 +231,6 @@ def test_rio_extent(): os.remove("test.tif") -def test_rio_extent2(): - with pytest.raises(AttributeError): - rio_profile(os.path.join(TEST_DIR, "test_geom", "lat.rdr")) - - def test_getTimeFromFile(): name1 = 'abcd_2020_01_01_T00_00_00jijk.xyz' assert getTimeFromFile(name1) == datetime.datetime(2020, 1, 1, 0, 0, 0) @@ -486,3 +483,99 @@ def test_get_nearest_wmtimes_4(): test_out = get_nearest_wmtimes(t0, 1) true_out = [datetime.datetime(2020, 1, 1, 11, 0), datetime.datetime(2020, 1, 1, 12, 0)] assert [t == t0 for t, t0 in zip(test_out, true_out)] + + +def test_rio(): + SCENARIO0_DIR = os.path.join(TEST_DIR, "scenario_0") + geotif = os.path.join(SCENARIO0_DIR, 'small_dem.tif') + profile = rio_profile(geotif) + assert profile['crs'] is not None + + +def test_rio_2(): + SCENARIO0_DIR = os.path.join(TEST_DIR, "scenario_0") + geotif = os.path.join(SCENARIO0_DIR, 'small_dem.tif') + prof = rio_profile(geotif) + del prof['transform'] + with pytest.raises(KeyError): + rio_extents(prof) + + +def test_rio_3(): + SCENARIO0_DIR = os.path.join(TEST_DIR, "scenario_0") + geotif = os.path.join(SCENARIO0_DIR, 'small_dem.tif') + data = rio_open(geotif, returnProj=False, userNDV=None, band=1) + assert data.shape == (569,558) + + +def test_rio_4(): + SCENARIO_DIR = os.path.join(TEST_DIR, "scenario_4") + los = os.path.join(SCENARIO_DIR, 'los.rdr') + inc, hd = rio_open(los, returnProj=False) + assert len(inc.shape) == 2 + assert len(hd.shape) == 2 + + +def test_writeArrayToRaster_2(): + test = np.random.randn(10,10,10) + with pytest.raises(RuntimeError): + writeArrayToRaster(test, 'dummy_file') + + +def test_writeArrayToRaster_3(tmp_path): + test = np.random.randn(10,10) + test = test + test * 1j + with pushd(tmp_path): + fname = os.path.join(tmp_path, 'tmp_file.tif') + writeArrayToRaster(test, fname) + tmp = rio_profile(fname) + assert tmp['dtype'] == np.complex64 + + +def test_writeArrayToRaster_3(tmp_path): + SCENARIO0_DIR = os.path.join(TEST_DIR, "scenario_0") + geotif = os.path.join(SCENARIO0_DIR, 'small_dem.tif') + profile = rio_profile(geotif) + data = rio_open(geotif) + with pushd(tmp_path): + fname = os.path.join(tmp_path, 'tmp_file.nc') + writeArrayToRaster( + data, + fname, + proj=profile['crs'], + gt=profile['transform'], + fmt='nc', + ) + new_fname = os.path.join(tmp_path, 'tmp_file.tif') + prof = rio_profile(new_fname) + assert prof['driver'] == 'GTiff' + + +def test_robs(): + assert robmin([1, 2, 3, np.nan])==1 + assert robmin([1,2,3])==1 + assert robmax([1, 2, 3, np.nan])==3 + assert robmax([1,2,3])==3 + + +def test_floorish1(): + assert np.isclose(floorish(5.6,0.2), 5.4) +def test_floorish2(): + assert np.isclose(floorish(5.71,0.2),5.6) +def test_floorish3(): + assert np.isclose(floorish(5.71,1),5) + +def test_projectDelays1(): + assert np.allclose(projectDelays(10,45),14.1421312) + + +def test_padLower(): + test = np.random.randn(2,3,4) + val = test[1,2,1] + test[1,2,0] = np.nan + out = padLower(test) + assert out[1,2,0] == val + + +def test_convertLons(): + assert np.allclose(convertLons(np.array([0, 10, -10, 190, 360])), np.array([0, 10, -10, -170, 0])) diff --git a/test/test_weather_model.py b/test/test_weather_model.py index dc2d0d1c2..6ead8cfc2 100644 --- a/test/test_weather_model.py +++ b/test/test_weather_model.py @@ -23,6 +23,7 @@ from RAiDER.models.gmao import GMAO from RAiDER.models.merra2 import MERRA2 from RAiDER.models.ncmr import NCMR +from RAiDER.models.customExceptions import * _LON0 = 0 @@ -156,12 +157,12 @@ def test_weatherModel_basic1(model: MockWeatherModel): wm.setTime('19720229', fmt='%Y%m%d') # test a leap year assert wm._time == datetime.datetime(1972, 2, 29, 0, 0, 0) - with pytest.raises(RuntimeError): + with pytest.raises(DatetimeOutsideRange): wm.checkTime(datetime.datetime(1950, 1, 1)) wm.checkTime(datetime.datetime(2000, 1, 1)) - with pytest.raises(RuntimeError): + with pytest.raises(DatetimeOutsideRange): wm.checkTime(datetime.datetime.now()) @@ -309,7 +310,7 @@ def test_hrrr(hrrr: HRRR): assert wm._Name == 'HRRR' assert wm._valid_range[0] == datetime.datetime(2016, 7, 15) assert wm._proj.to_epsg() is None - with pytest.raises(RuntimeError): + with pytest.raises(DatetimeOutsideRange): wm.checkTime(datetime.datetime(2010, 7, 15)) wm.checkTime(datetime.datetime(2018, 7, 12)) @@ -329,7 +330,7 @@ def test_hrrrak(hrrrak: HRRRAK): with pytest.raises(ValueError): wm.checkValidBounds([15, 20, 265, 270]) - with pytest.raises(RuntimeError): + with pytest.raises(DatetimeOutsideRange): wm.checkTime(datetime.datetime(2018, 7, 12)) wm.checkTime(datetime.datetime(2018, 7, 15)) @@ -431,7 +432,6 @@ def test_hrrr_badloc(wm:hrrr=HRRR): with pytest.raises(ValueError): wm._fetch('dummy_filename') - def test_hrrrak_dl(tmp_path: Path, wm:hrrrak=HRRRAK): wm = wm() d = tmp_path / "files" @@ -443,7 +443,6 @@ def test_hrrrak_dl(tmp_path: Path, wm:hrrrak=HRRRAK): wm._fetch(fname) assert True - def test_hrrrak_dl2(tmp_path: Path, wm:hrrrak=HRRRAK): # test the international date line crossing wm = wm() diff --git a/tools/RAiDER/aria/prepFromGUNW.py b/tools/RAiDER/aria/prepFromGUNW.py index c2b4c984e..4561abe4f 100644 --- a/tools/RAiDER/aria/prepFromGUNW.py +++ b/tools/RAiDER/aria/prepFromGUNW.py @@ -11,21 +11,19 @@ import eof.download import xarray as xr import rasterio -import geopandas as gpd import pandas as pd import yaml import shapely.wkt from dataclasses import dataclass import sys from shapely.geometry import box -from rasterio.crs import CRS import RAiDER -from RAiDER.utilFcns import rio_open, writeArrayToRaster from RAiDER.logger import logger from RAiDER.models import credentials from RAiDER.models.hrrr import HRRR_CONUS_COVERAGE_POLYGON, AK_GEO, check_hrrr_dataset_availability from RAiDER.s1_azimuth_timing import get_times_for_azimuth_interpolation +from RAiDER.s1_orbits import ensure_orbit_credentials ## cube spacing in degrees for each model DCT_POSTING = {'HRRR': 0.05, 'HRES': 0.10, 'GMAO': 0.10, 'ERA5': 0.10, 'ERA5T': 0.10} @@ -278,9 +276,10 @@ def get_orbit_file(self): sat = slc.split('_')[0] dt = datetime.strptime(f'{self.dates[0]}T{self.mid_time}', '%Y%m%dT%H:%M:%S') + ensure_orbit_credentials() path_orb = eof.download.download_eofs([dt], [sat], save_dir=orbit_dir) - return path_orb + return [str(o) for o in path_orb] ## ------ methods below are not used @@ -381,7 +380,7 @@ def update_yaml(dct_cfg:dict, dst:str='GUNW.yaml'): params = yaml.safe_load(f) except yaml.YAMLError as exc: print(exc) - raise ValueError(f'Something is wrong with the yaml file {example_yaml}') + raise ValueError(f'Something is wrong with the yaml file {template_file}') params = {**params, **dct_cfg} diff --git a/tools/RAiDER/cli/__main__.py b/tools/RAiDER/cli/__main__.py index 447b82bd4..45a656212 100644 --- a/tools/RAiDER/cli/__main__.py +++ b/tools/RAiDER/cli/__main__.py @@ -29,11 +29,9 @@ def main(): sys.argv = [args.process, *unknowns] - from RAiDER.cli.raider import calcDelays, downloadGNSS, calcDelaysGUNW - try: # python >=3.10 interface - process_entry_point = entry_points(group='console_scripts', name=f'{args.process}.py')[0] + (process_entry_point,) = entry_points(group='console_scripts', name=f'{args.process}.py') except TypeError: # python 3.8 and 3.9 interface scripts = entry_points()['console_scripts'] diff --git a/tools/RAiDER/cli/raider.py b/tools/RAiDER/cli/raider.py index ea90c3363..4034d8da4 100644 --- a/tools/RAiDER/cli/raider.py +++ b/tools/RAiDER/cli/raider.py @@ -20,12 +20,15 @@ from RAiDER.cli.parser import add_out, add_cpus, add_verbose from RAiDER.cli.validators import DateListAction, date_type from RAiDER.models.allowed import ALLOWED_MODELS +from RAiDER.models.customExceptions import * from RAiDER.utilFcns import get_dt from RAiDER.s1_azimuth_timing import get_s1_azimuth_time_grid, get_inverse_weights_for_dates, get_times_for_azimuth_interpolation import traceback +TIME_INTERPOLATION_METHODS = ['none', 'center_time', 'azimuth_time_grid'] + HELP_MESSAGE = """ Command line options for RAiDER processing. Default options can be found by running @@ -241,7 +244,8 @@ def calcDelays(iargs=None): ): ########################################################### - # weather model calculation + # Weather model calculation + ########################################################### logger.debug('Starting to run the weather model calculation') logger.debug(f'Requested date,time: {t.strftime("%Y%m%d, %H:%M")}') logger.debug('Beginning weather model pre-processing') @@ -252,13 +256,11 @@ def calcDelays(iargs=None): logger.warning('interp_method is not specified, defaulting to \'none\', i.e. nearest datetime for delay ' 'calculation') - # Grab the closest two times unless the user specifies 'nearest' via 'none' or None. - # If the model time_delta is not specified then use 6 - # The two datetimes will be combined to a single file and processed - # TODO: make more transparent control flow for GUNW and non-GUNW workflow - if (interp_method in ['none', 'center_time']): - times = get_nearest_wmtimes(t, [model.dtime() if \ - model.dtime() is not None else 6][0]) if interp_method == 'center_time' else [t] + if (interp_method != 'azimuth_time_grid'): + times = get_nearest_wmtimes( + t, [model.dtime() if model.dtime() is not None else 6][0] + ) if interp_method == 'center_time' else [t] + elif interp_method == 'azimuth_time_grid': step = model.dtime() time_step_hours = step if step is not None else 6 @@ -270,129 +272,37 @@ def calcDelays(iargs=None): wfiles = [] for tt in times: try: - wfile = RAiDER.processWM.prepareWeatherModel(model, - tt, - aoi.bounds(), - makePlots=params['verbose']) + wfile = RAiDER.processWM.prepareWeatherModel( + model, + tt, + aoi.bounds(), + makePlots=params['verbose'] + ) wfiles.append(wfile) - # catch when requested datetime fails - except RuntimeError as re: - continue + except TryToKeepGoingError: + if interp_method in ['azimuth_time_grid', 'none']: + raise DatetimeFailed(model.Model(), tt) + else: + continue - # catch when something else within weather model class fails + # log when something else happens and then re-raise the error except Exception as e: S, N, W, E = wm_bounds logger.info(f'Weather model point bounds are {S:.2f}/{N:.2f}/{W:.2f}/{E:.2f}') logger.info(f'Query datetime: {tt}') msg = f'Downloading and/or preparation of {model._Name} failed.' logger.error(e) + logger.error('Weather model files are: {}'.format(wfiles)) logger.error(msg) - if interp_method == 'azimuth_time_grid': - break + raise + # dont process the delays for download only if dl_only: continue - if (len(wfiles) == 0) and (interp_method != 'azimuth_time_grid'): - logger.error('No weather model data was successfully processed.') - if len(params['date_list']) == 1: - raise RuntimeError - # skip date and continue processing if multiple dates are requested - else: - continue - - # nearest weather model time via 'none' is specified - # When interp_method is 'none' only 1 weather model file and one relevant time - elif (interp_method == 'none') and (len(wfiles) == 1) and (len(times) == 1): - weather_model_file = wfiles[0] - - # only one time in temporal interpolation worked - # TODO: this seems problematic - unexpected behavior possibly for 'center_time' - elif (len(wfiles) == 1) and (len(times) == 2) and (interp_method != 'azimuth_time_grid'): - logger.warning('Time interpolation did not succeed, defaulting to nearest available date') - weather_model_file = wfiles[0] - - # TODO: ensure this additional conditional is appropriate; assuming wfiles == 2 ONLY for 'center_time' - # value of 'interp_method' parameter - elif (interp_method == 'center_time') and (len(wfiles) == 2): - ds1 = xr.open_dataset(wfiles[0]) - ds2 = xr.open_dataset(wfiles[1]) - - # calculate relative weights of each dataset - date1 = datetime.datetime.strptime(ds1.attrs['datetime'], '%Y_%m_%dT%H_%M_%S') - date2 = datetime.datetime.strptime(ds2.attrs['datetime'], '%Y_%m_%dT%H_%M_%S') - wgts = [ 1 - get_dt(t, date1) / get_dt(date2, date1), 1 - get_dt(date2, t) / get_dt(date2, date1)] - try: - assert np.isclose(np.sum(wgts), 1) - except AssertionError: - logger.error('Time interpolation weights do not sum to one; something is off with query datetime: %s', t) - continue - - # combine datasets - ds = ds1 - for var in ['wet', 'hydro', 'wet_total', 'hydro_total']: - ds[var] = (wgts[0] * ds1[var]) + (wgts[1] * ds2[var]) - ds.attrs['Date1'] = 0 - ds.attrs['Date2'] = 0 - weather_model_file = os.path.join( - os.path.dirname(wfiles[0]), - os.path.basename(wfiles[0]).split('_')[0] + '_' + t.strftime('%Y_%m_%dT%H_%M_%S') + '_timeInterp_' + '_'.join(wfiles[0].split('_')[-4:]), - ) - ds.to_netcdf(weather_model_file) - elif (interp_method == 'azimuth_time_grid'): - n_files = len(wfiles) - n_times = len(times) - if n_files != n_times: - raise ValueError('The model files for the datetimes for requisite azimuth interpolation were not ' - 'succesfully downloaded or processed') - datasets = [xr.open_dataset(f) for f in wfiles] - - # Each model will require some inspection here - # the subsequent s1 azimuth time grid requires dimension - # inputs to all have same dimensions and either be - # 1d or 3d. - if model._dataset in ['hrrr', 'hrrrak']: - lat_2d = datasets[0].latitude.data - lon_2d = datasets[0].longitude.data - z_1d = datasets[0].z.data - m, n, p = z_1d.shape[0], lat_2d.shape[0], lat_2d.shape[1] - - lat = np.broadcast_to(lat_2d, (m, n, p)) - lon = np.broadcast_to(lon_2d, (m, n, p)) - hgt = np.broadcast_to(z_1d[:, None, None], (m, n, p)) - - else: - raise NotImplementedError('Azimuth Time is currently only implemented for HRRR') - - time_grid = RAiDER.s1_azimuth_timing.get_s1_azimuth_time_grid(lon, - lat, - hgt, - # This is the acq time from loop - t) - - if np.any(np.isnan(time_grid)): - raise ValueError('The azimuth time grid return nans meaning no orbit was downloaded.') - wgts = get_inverse_weights_for_dates(time_grid, - times, - temporal_window_hours=model._time_res) - # combine datasets - ds_out = datasets[0] - for var in ['wet', 'hydro', 'wet_total', 'hydro_total']: - ds_out[var] = sum([wgt * ds[var] for (wgt, ds) in zip(wgts, datasets)]) - ds_out.attrs['Date1'] = 0 - ds_out.attrs['Date2'] = 0 - weather_model_file = os.path.join( - os.path.dirname(wfiles[0]), - # TODO: clean up - os.path.basename(wfiles[0]).split('_')[0] + '_' + t.strftime('%Y_%m_%dT%H_%M_%S') + '_timeInterpAziGrid_' + '_'.join(wfiles[0].split('_')[-4:]), - ) - ds_out.to_netcdf(weather_model_file) - # TODO: test to ensure this error is caught - else: - n = len(wfiles) - raise NotImplementedError(f'The {interp_method} with {n} retrieved weather model files was not well posed ' - 'for the the delay workflow.') + # Get the weather model file + weather_model_file = getWeatherFile(wfiles, times, t, model._Name, interp_method) # Now process the delays try: @@ -576,7 +486,7 @@ def calcDelaysGUNW(iargs: list[str] = None) -> xr.Dataset: p.add_argument( '-interp', '--interpolate-time', default='azimuth_time_grid', type=str, - choices=['none', 'center_time', 'azimuth_time_grid'], + choices=TIME_INTERPOLATION_METHODS, help=('How to interpolate across model time steps. Possible options are: ' '[\'none\', \'center_time\', \'azimuth_time_grid\'] ' 'None: means nearest model time; center_time: linearly across center time; ' @@ -606,7 +516,7 @@ def calcDelaysGUNW(iargs: list[str] = None) -> xr.Dataset: file_name = iargs.file.split('/')[-1] gunw_id = file_name.replace('.nc', '') if not RAiDER.aria.prepFromGUNW.check_hrrr_dataset_availablity_for_s1_azimuth_time_interpolation(gunw_id): - raise ValueError('The required HRRR data for time-grid interpolation is not available') + raise NoWeatherModelData('The required HRRR data for time-grid interpolation is not available') if not iargs.file and iargs.bucket: # only use GUNW ID for checking if HRRR available @@ -682,3 +592,198 @@ def combineZTDFiles(): outName=args.out_name, localTime=args.local_time ) + + +def getWeatherFile(wfiles, times, t, model, interp_method='none'): + ''' + # Time interpolation + # + # Need to handle various cases, including if the exact weather model time is + # requested, or if one or more datetimes are not available from the weather + # model data provider + ''' + + # time interpolation method: number of expected files + EXPECTED_NUM_FILES = {'none': 1, 'center_time': 2, 'azimuth_time_grid': 3} + + Nfiles = len(wfiles) + Ntimes = len(times) + + try: + Nfiles_expected = EXPECTED_NUM_FILES[interp_method] + except KeyError: + raise ValueError('getWeatherFile: interp_method {} is not known'.format(interp_method)) + + Nmatch = (Nfiles_expected == Nfiles) + Tmatch = (Nfiles == Ntimes) + + # Case 1: no files downloaded + if Nfiles==0: + logger.error('No weather model data was successfully processed.') + return None + + # Case 2 - nearest weather model time is requested and retrieved + if (interp_method == 'none'): + weather_model_file = wfiles[0] + + + elif (interp_method == 'center_time'): + + if Nmatch: # Case 3: two weather files downloaded + weather_model_file = combine_weather_files( + wfiles, + t, + model, + interp_method='center_time' + ) + elif Tmatch: # Case 4: Exact time is available without interpolation + logger.warning('Time interpolation is not needed as exact time is available') + weather_model_file = wfiles[0] + elif Nfiles == 1: # Case 5: one file does not download for some reason + logger.warning('getWeatherFile: One datetime is not available to download, defaulting to nearest available date') + weather_model_file = wfiles[0] + else: + raise WrongNumberOfFiles(Nfiles_expected, Nfiles) + + elif (interp_method) == 'azimuth_time_grid': + + if Nmatch or Tmatch: # Case 6: all files downloaded + weather_model_file = combine_weather_files( + wfiles, + t, + model, + interp_method='azimuth_time_grid' + ) + else: + raise WrongNumberOfFiles(Nfiles_expected, Nfiles) + + # Case 7 - Anything else errors out + else: + N = len(wfiles) + raise NotImplementedError(f'The {interp_method} with {N} retrieved weather model files was not well posed ' + 'for the current workflow.') + + return weather_model_file + + +def combine_weather_files(wfiles, t, model, interp_method='center_time'): + '''Interpolate downloaded weather files and save to a single file''' + + STYLE = {'center_time': '_timeInterp_', 'azimuth_time_grid': '_timeInterpAziGrid_'} + + # read the individual datetime datasets + datasets = [xr.open_dataset(f) for f in wfiles] + + # Pull the datetimes from the datasets + times = [] + for ds in datasets: + times.append(datetime.datetime.strptime(ds.attrs['datetime'], '%Y_%m_%dT%H_%M_%S')) + + if len(times)==0: + raise NoWeatherModelData() + + # calculate relative weights of each dataset + if interp_method == 'center_time': + wgts = get_weights_time_interp(times, t) + elif interp_method == 'azimuth_time_grid': + time_grid = get_time_grid_for_aztime_interp(datasets, t, model) + wgts = get_inverse_weights_for_dates(time_grid, times) + + # combine datasets + ds_out = datasets[0] + for var in ['wet', 'hydro', 'wet_total', 'hydro_total']: + ds_out[var] = sum([wgt * ds[var] for (wgt, ds) in zip(wgts, datasets)]) + ds_out.attrs['Date1'] = 0 + ds_out.attrs['Date2'] = 0 + + # Give the weighted combination a new file name + weather_model_file = os.path.join( + os.path.dirname(wfiles[0]), + os.path.basename(wfiles[0]).split('_')[0] + '_' + + t.strftime('%Y_%m_%dT%H_%M_%S') + STYLE[interp_method] + + '_'.join(wfiles[0].split('_')[-4:]), + ) + + # write the combined results to disk + ds_out.to_netcdf(weather_model_file) + + return weather_model_file + + +def combine_files_using_azimuth_time(wfiles, t, times): + '''Combine files using azimuth time interpolation''' + + # read the individual datetime datasets + datasets = [xr.open_dataset(f) for f in wfiles] + + # Pull the datetimes from the datasets + times = [] + for ds in datasets: + times.append(datetime.datetime.strptime(ds.attrs['datetime'], '%Y_%m_%dT%H_%M_%S')) + + model = datasets[0].attrs['model_name'] + + time_grid = get_time_grid_for_aztime_interp(datasets, times, t, model) + + wgts = get_inverse_weights_for_dates(time_grid, times) + + # combine datasets + ds_out = datasets[0] + for var in ['wet', 'hydro', 'wet_total', 'hydro_total']: + ds_out[var] = sum([wgt * ds[var] for (wgt, ds) in zip(wgts, datasets)]) + ds_out.attrs['Date1'] = 0 + ds_out.attrs['Date2'] = 0 + + # Give the weighted combination a new file name + weather_model_file = os.path.join( + os.path.dirname(wfiles[0]), + os.path.basename(wfiles[0]).split('_')[0] + '_' + t.strftime('%Y_%m_%dT%H_%M_%S') + '_timeInterpAziGrid_' + '_'.join(wfiles[0].split('_')[-4:]), + ) + + # write the combined results to disk + ds_out.to_netcdf(weather_model_file) + + return weather_model_file + + +def get_weights_time_interp(times, t): + '''Calculate weights for time interpolation using simple inverse linear weighting''' + date1,date2 = times + wgts = [ 1 - get_dt(t, date1) / get_dt(date2, date1), 1 - get_dt(date2, t) / get_dt(date2, date1)] + + try: + assert np.isclose(np.sum(wgts), 1) + except AssertionError: + logger.error('Time interpolation weights do not sum to one; something is off with query datetime: %s', t) + return None + + return wgts + + +def get_time_grid_for_aztime_interp(datasets, t, model): + '''Calculate the time-varying grid for use with azimuth time interpolation''' + + # Each model will require some inspection here + # the subsequent s1 azimuth time grid requires dimension + # inputs to all have same dimensions and either be + # 1d or 3d. + AZ_TIME_ALLOWED_MODELS = ['hrrr', 'hrrrak'] + + if model.lower() in AZ_TIME_ALLOWED_MODELS: + lat_2d = datasets[0].latitude.data + lon_2d = datasets[0].longitude.data + z_1d = datasets[0].z.data + m, n, p = z_1d.shape[0], lat_2d.shape[0], lat_2d.shape[1] + + lat = np.broadcast_to(lat_2d, (m, n, p)) + lon = np.broadcast_to(lon_2d, (m, n, p)) + hgt = np.broadcast_to(z_1d[:, None, None], (m, n, p)) + + else: + raise NotImplementedError('Azimuth Time is currently only implemented for HRRR') + + time_grid = get_s1_azimuth_time_grid(lon, lat, hgt, t) # This is the acq time from loop + if np.any(np.isnan(time_grid)): + raise ValueError('The Time Grid return nans meaning no orbit was downloaded.') + + return time_grid diff --git a/tools/RAiDER/cli/validators.py b/tools/RAiDER/cli/validators.py index c5d8847d4..78c470248 100755 --- a/tools/RAiDER/cli/validators.py +++ b/tools/RAiDER/cli/validators.py @@ -165,7 +165,6 @@ def get_query_region(args): # TODO: Need to incorporate the cube raise ValueError('No valid query points or bounding box found in the configuration file') - return query diff --git a/tools/RAiDER/delayFcns.py b/tools/RAiDER/delayFcns.py index c56a2f01c..a71de3b21 100755 --- a/tools/RAiDER/delayFcns.py +++ b/tools/RAiDER/delayFcns.py @@ -79,12 +79,3 @@ def make_shared_raw(inarr): return shared_arr_np -def interpolate2(fun, x, y, z): - ''' - helper function to make the interpolation step cleaner - ''' - in_shape = x.shape - out = fun((y.ravel(), x.ravel(), z.ravel())) # note that this re-ordering is on purpose to match the weather model - outData = out.reshape(in_shape) - return outData - diff --git a/tools/RAiDER/dem.py b/tools/RAiDER/dem.py index c21338ea7..7d6099c46 100644 --- a/tools/RAiDER/dem.py +++ b/tools/RAiDER/dem.py @@ -14,47 +14,17 @@ import rasterio from dem_stitcher.stitcher import stitch_dem -from RAiDER.interpolator import interpolateDEM from RAiDER.logger import logger -from RAiDER.utilFcns import rio_open, get_file_and_band - - -def getHeights(ll_bounds, dem_type, dem_file, lats=None, lons=None): - ''' - Fcn to return heights from a DEM, either one that already exists - or will download one if needed. - ''' - # height_type, height_data = heights - if dem_type == 'hgt': - htinfo = get_file_and_band(dem_file) - hts = rio_open(htinfo[0], band=htinfo[1]) - - elif dem_type == 'csv': - # Heights are in the .csv file - hts = pd.read_csv(dem_file)['Hgt_m'].values - - elif dem_type == 'interpolate': - # heights will be vertically interpolated to the heightlvs - hts = None - - elif (dem_type == 'download') or (dem_type == 'dem'): - zvals, metadata = download_dem(ll_bounds, writeDEM=True, outName=dem_file) - - #TODO: check this - lons, lats = np.meshgrid(lons, lats) - # Interpolate to the query points - hts = interpolateDEM(zvals, metadata['transform'], (lats, lons), method='nearest') - - return hts +from RAiDER.utilFcns import rio_open def download_dem( - ll_bounds, - writeDEM=False, - outName='warpedDEM', - buf=0.02, - overwrite=False, -): + ll_bounds=None, + demName='warpedDEM.dem', + overwrite=False, + writeDEM=False, + buf=0.02, + ): """ Download a DEM if one is not already present. Args: @@ -67,11 +37,22 @@ def download_dem( zvals: np.array -DEM heights metadata: -metadata for the DEM """ - if os.path.exists(outName) and not overwrite: - logger.info('Using existing DEM: %s', outName) - zvals, metadata = rio_open(outName, returnProj=True) + if os.path.exists(demName): + if overwrite: + download = True + else: + download = False + else: + download = True + + if download and (ll_bounds is None): + raise ValueError('download_dem: Either an existing file or lat/lon bounds must be passed') + if not download: + logger.info('Using existing DEM: %s', demName) + zvals, metadata = rio_open(demName, returnProj=True) else: + # download the dem # inExtent is SNWE # dem-stitcher wants WSEN bounds = [ @@ -86,9 +67,9 @@ def download_dem( dst_area_or_point='Area', ) if writeDEM: - with rasterio.open(outName, 'w', **metadata) as ds: + with rasterio.open(demName, 'w', **metadata) as ds: ds.write(zvals, 1) ds.update_tags(AREA_OR_POINT='Point') - logger.info('Wrote DEM: %s', outName) + logger.info('Wrote DEM: %s', demName) return zvals, metadata diff --git a/tools/RAiDER/getStationDelays.py b/tools/RAiDER/getStationDelays.py index 5b0cc4911..f61f8b8db 100644 --- a/tools/RAiDER/getStationDelays.py +++ b/tools/RAiDER/getStationDelays.py @@ -24,28 +24,42 @@ def get_delays_UNR(stationFile, filename, dateList, returnTime=None): Parses and returns a dictionary containing either (1) all the GPS delays, if returnTime is None, or (2) only the delay at the closest times to to returnTime. - Inputs: - stationFile - a .gz station delay file - returnTime - specified time of GPS delay - Outputs: - a dict and CSV file containing the times and delay information - (delay in mm, delay uncertainty, delay gradients) - *NOTE: Due to a formatting error in the tropo SINEX files, the two tropospheric gradient columns - (TGNTOT and TGETOT) are interchanged, as are the formal error columns (_SIG). + + Args: + stationFile: binary - a .gz station delay file + filename: ? - ? + dateList: list of datetime - ? + returnTime: datetime - specified time of GPS delay (default all times) + + Returns: + None + + The function writes a CSV file containing the times and delay information + (delay in mm, delay uncertainty, delay gradients) + + Refer to the following sites to interpret stationFile variable names: + ftp://igs.org/pub/data/format/sinex_tropo.txt + http://geodesy.unr.edu/gps_timeseries/README_trop2.txt + Wet and hydrostratic delays were derived as so: + Constants —> k1 = 0.704, k2 = 0.776, k3 = 3739.0, m = 18.0152/28.9644, + k2' = k2-(k1*m) = 0.33812796398337275, Rv = 461.5 J/(kg·K), ρl = 997 kg/m^3 + + *NOTE: wet delays passed here are computed using + PMV = precipitable water vapor, + P = total atm pressure, + Tm = mean temp of the column, as: + + Wet zenith delay = 10^-6 ρlRv(k2' + k3/Tm) PMV + Hydrostatic zenith delay = Total zenith delay - wet zenith delay = k1*(P/Tm) + + Source —> Hanssen, R. F. (2001) eqns. 6.2.7-10 + + *NOTE: Due to a formatting error in the tropo SINEX files, the two + tropospheric gradient columns (TGNTOT and TGETOT) are interchanged, + as are the formal error columns (_SIG). + Source —> http://geodesy.unr.edu/gps_timeseries/README_trop2.txt) ''' - # Refer to the following sites to interpret stationFile variable names: - # ftp://igs.org/pub/data/format/sinex_tropo.txt - # http://geodesy.unr.edu/gps_timeseries/README_trop2.txt - # Wet and hydrostratic delays were derived as so: - # Constants —> k1 = 0.704, k2 = 0.776, k3 = 3739.0, m = 18.0152/28.9644, - # k2' = k2-(k1*m) = 0.33812796398337275, Rv = 461.5 J/(kg·K), ρl = 997 kg/m^3 - # Note wet delays passed here may be computed as so - # where PMV = precipitable water vapor, P = total atm pressure, Tm = mean temp of the column —> - # Wet zenith delay = 10^-6 ρlRv(k2' + k3/Tm) PMV - # Hydrostatic zenith delay = Total zenith delay - wet zenith delay = k1*(P/Tm) - # Source —> Hanssen, R. F. (2001) eqns. 6.2.7-10 - # sort through station zip files allstationTarfiles = [] # if URL diff --git a/tools/RAiDER/gnss/downloadGNSSDelays.py b/tools/RAiDER/gnss/downloadGNSSDelays.py index 04f1ed283..c4792beb6 100755 --- a/tools/RAiDER/gnss/downloadGNSSDelays.py +++ b/tools/RAiDER/gnss/downloadGNSSDelays.py @@ -13,79 +13,99 @@ from RAiDER.logger import logger, logging from RAiDER.getStationDelays import get_station_data from RAiDER.utilFcns import requests_retry_session +from RAiDER.models.customExceptions import NoStationDataFoundError # base URL for UNR repository _UNR_URL = "http://geodesy.unr.edu/" -def get_station_list(bbox=None, writeLoc=None, userstatList=None, name_appendix=''): +def get_station_list( + bbox=None, + stationFile=None, + userstatList=None, + writeStationFile=True, + writeLoc=None, + name_appendix='' + ): ''' Creates a list of stations inside a lat/lon bounding box from a source - Inputs: - bbox - length-4 list of floats that describes a bounding box. Format is - S N W E + + Args: + bbox: list of float - length-4 list of floats that describes a bounding box. + Format is S N W E + writeLoc: string - Directory to write data + userstatList: list - list of specific IDs to access + name_appendix: str - name to append to output file + + Returns: + stations: list of strings - station IDs to access + output_file: string - file to write delays ''' - writeLoc = os.path.join(writeLoc or os.getcwd(), 'gnssStationList_overbbox' + name_appendix + '.csv') - - if userstatList: - userstatList = read_text_file(userstatList) - - statList = get_stats_by_llh(llhBox=bbox, userstatList=userstatList) + if stationFile: + station_data = pd.read_csv(stationFile) + elif userstatList: + station_data = read_text_file(userstatList) + elif bbox: + station_data = get_stats_by_llh(llhBox=bbox) # write to file and pass final stations list - statList.to_csv(writeLoc, index=False) - stations = list(statList['ID'].values) + if writeStationFile: + output_file = os.path.join(writeLoc or os.getcwd(), 'gnssStationList_overbbox' + name_appendix + '.csv') + station_data.to_csv(output_file, index=False) - return stations, writeLoc + return list(station_data['ID'].values), [output_file if writeStationFile else station_data][0] -def get_stats_by_llh(llhBox=None, baseURL=_UNR_URL, userstatList=None): +def get_stats_by_llh(llhBox=None, baseURL=_UNR_URL): ''' Function to pull lat, lon, height, beginning date, end date, and number of solutions for stations inside the bounding box llhBox. - llhBox should be a tuple with format (lat1, lat2, lon1, lon2), where lat1, lon1 define the lower left-hand corner and lat2, lon2 - define the upper right corner. + llhBox should be a tuple in SNWE format. ''' if llhBox is None: llhBox = [-90, 90, 0, 360] + S,N,W,E = llhBox + if (W < 0) or (E < 0): + raise ValueError('get_stats_by_llh: bounding box must use 0 < lon < 360') - stationHoldings = '{}NGLStationPages/DataHoldings.txt'.format(baseURL) + stationHoldings = '{}NGLStationPages/llh.out'.format(baseURL) # it's a file like object and works just like a file - session = requests_retry_session() - data = session.get(stationHoldings) - stations = [] - for ind, line in enumerate(data.text.splitlines()): # files are iterable - if ind == 0: - continue - statID, lat, lon, height = get_ID(line) - # Only pass if in bbox - # And if user list of stations specified, only pass info for stations within list - if in_box(lat, lon, llhBox) and (not userstatList or statID in userstatList): - # convert lon into range [-180,180] - lon = fix_lons(lon) - stations.append({'ID': statID, 'Lat': lat, 'Lon': lon, 'Hgt_m': height}) - - logger.info('%d stations were found in %s', len(stations), llhBox) - stations = pd.DataFrame(stations) - # Report stations from user's list that do not cover bbox - if userstatList: - userstatList = [ - i for i in userstatList if i not in stations['ID'].to_list()] - if userstatList: - logger.warning( - "The following user-input stations are not covered by the input " - "bounding box %s: %s", - str(llhBox).strip('[]'), str(userstatList).strip('[]') - ) + stations = pd.read_csv( + stationHoldings, + delim_whitespace=True, + names=['ID', 'Lat', 'Lon', 'Hgt_m'] + ) + stations = filterToBBox(stations, llhBox) + stations['Lon'] = ((stations['Lon'].values + 180) % 360) - 180 # convert lons to -180 - 180 + + stations = filterToBBox(stations, llhBox) return stations -def download_tropo_delays(stats, years, gps_repo=None, writeDir='.', numCPUs=8, download=False): +def download_tropo_delays( + stats, years, + gps_repo='UNR', + writeDir='.', + numCPUs=8, + download=False, + ): ''' - Check for and download GNSS tropospheric delays from an archive. If download is True then - files will be physically downloaded, which again is not necessary as data can be virtually accessed. + Check for and download GNSS tropospheric delays from an archive. If + download is True then files will be physically downloaded, but this + is not necessary as data can be virtually accessed. + + Args: + stats: stations - Stations to access + years: list of int - A list of years to be downloaded + gps_repo: string - SNWE bounds target area to ensure weather model contains them + writeDir: string - False if preprocessing weather model data + numCPUs: int - whether to write debug plots + download: bool - True if you want to download even when the weather model exists + + Returns: + None ''' # argument checking @@ -97,7 +117,9 @@ def download_tropo_delays(stats, years, gps_repo=None, writeDir='.', numCPUs=8, # Iterate over stations and years and check or download data stat_year_tup = itertools.product(stats, years) stat_year_tup = ((*tup, writeDir, download) for tup in stat_year_tup) + # Parallelize remote querying of station locations + results = [] with multiprocessing.Pool(numCPUs) as multipool: # only record valid path if gps_repo == 'UNR': @@ -105,8 +127,12 @@ def download_tropo_delays(stats, years, gps_repo=None, writeDir='.', numCPUs=8, fileurl for fileurl in multipool.starmap(download_UNR, stat_year_tup) if fileurl['path'] ] + else: + raise NotImplementedError('download_tropo_delays: gps_repo "{}" not yet implemented'.format(gps_repo)) # Write results to file + if len(results)==0: + raise NoStationDataFoundError(station_list=stats['ID'].to_list(), years=years) statDF = pd.DataFrame(results).set_index('ID') statDF.to_csv(os.path.join(writeDir, '{}gnssStationList_overbbox_withpaths.csv'.format(gps_repo))) @@ -224,52 +250,14 @@ def main(inps=None): # Setup bounding box if bounding_box: - if isinstance(bounding_box, str) and not os.path.isfile(bounding_box): - try: - bbox = [float(val) for val in bounding_box.split()] - except ValueError: - raise Exception( - 'Cannot understand the --bbox argument. String input is incorrect or path does not exist.') - elif isinstance(bounding_box, list): - bbox = bounding_box - - else: - raise Exception('Passing a file with a bounding box not yet supported.') - - long_cross_zero = 1 if bbox[2] * bbox[3] < 0 else 0 - - # if necessary, convert negative longitudes to positive - if bbox[2] < 0: - bbox[2] += 360 - - if bbox[3] < 0: - bbox[3] += 360 - + bbox, long_cross_zero = parse_bbox(bounding_box) # If bbox not specified, query stations across the entire globe else: bbox = [-90, 90, 0, 360] long_cross_zero = 1 # Handle station query - if long_cross_zero == 1: - bbox1 = bbox.copy() - bbox2 = bbox.copy() - bbox1[3] = 360.0 - bbox2[2] = 0.0 - stats1, origstatsFile1 = get_station_list(bbox=bbox1, writeLoc=out, userstatList=station_file, name_appendix='_a') - stats2, origstatsFile2 = get_station_list(bbox=bbox2, writeLoc=out, userstatList=station_file, name_appendix='_b') - stats = stats1 + stats2 - origstatsFile = origstatsFile1[:-6] + '.csv' - file_a = pd.read_csv(origstatsFile1) - file_b = pd.read_csv(origstatsFile2) - frames = [file_a, file_b] - result = pd.concat(frames, ignore_index=True) - result.to_csv(origstatsFile, index=False) - else: - if bbox[3] < bbox[2]: - bbox[3] = 360.0 - stats, origstatsFile = get_station_list( - bbox=bbox, writeLoc=out, userstatList=station_file) + stats = get_stats(bbox, long_cross_zero, out, station_file) # iterate over years years = list(set([i.year for i in dateList])) @@ -299,3 +287,100 @@ def main(inps=None): ) logger.debug('Completed processing') + + +def parse_bbox(bounding_box): + ''' + Parse bounding box arguments + ''' + if isinstance(bounding_box, str) and not os.path.isfile(bounding_box): + try: + bbox = [float(val) for val in bounding_box.split()] + except ValueError: + raise Exception( + 'Cannot understand the --bbox argument. String input is incorrect or path does not exist.') + elif isinstance(bounding_box, list): + bbox = bounding_box + + else: + raise Exception('Passing a file with a bounding box not yet supported.') + + long_cross_zero = 1 if bbox[2] * bbox[3] < 0 else 0 + + # if necessary, convert negative longitudes to positive + if bbox[2] < 0: + bbox[2] += 360 + + if bbox[3] < 0: + bbox[3] += 360 + + return bbox, long_cross_zero + + +def get_stats(bbox, long_cross_zero, out, station_file): + ''' + Pull the stations needed + ''' + if long_cross_zero == 1: + bbox1 = bbox.copy() + bbox2 = bbox.copy() + bbox1[3] = 360.0 + bbox2[2] = 0.0 + stats1, origstatsFile1 = get_station_list(bbox=bbox1, writeLoc=out, stationFile=station_file, name_appendix='_a') + stats2, origstatsFile2 = get_station_list(bbox=bbox2, writeLoc=out, stationFile=station_file, name_appendix='_b') + stats = stats1 + stats2 + origstatsFile = origstatsFile1[:-6] + '.csv' + file_a = pd.read_csv(origstatsFile1) + file_b = pd.read_csv(origstatsFile2) + frames = [file_a, file_b] + result = pd.concat(frames, ignore_index=True) + result.to_csv(origstatsFile, index=False) + else: + if bbox[3] < bbox[2]: + bbox[3] = 360.0 + stats, origstatsFile = get_station_list( + bbox=bbox, + writeLoc=out, + stationFile=station_file + ) + + return stats + + +def filterToBBox(stations, llhBox): + ''' + Filter a dataframe by lat/lon. + *NOTE: llhBox longitude format should be 0-360 + + Args: + stations: DataFrame - a pandas dataframe with "Lat" and "Lon" columns + llhBox: list of float - 4-element list: [S, N, W, E] + + Returns: + a Pandas Dataframe with stations removed that are not inside llhBox + ''' + S,N,W,E = llhBox + if (W<0) or (E<0): + raise ValueError('llhBox longitude format should 0-360') + + # For a user-provided file, need to check the column names + keys = stations.columns + lat_keys = ['lat', 'latitude', 'Lat', 'Latitude'] + lon_keys = ['lon', 'longitude', 'Lon', 'Longitude'] + index = None + for k,key in enumerate(lat_keys): + if key in list(keys): + index = k + break + if index is None: + raise KeyError('filterToBBox: No valid column names found for latitude and longitude') + lon_key = lon_keys[k] + lat_key = lat_keys[k] + + if stations[lon_key].min() < 0: + # convert lon format to -180 to 180 + W,E = [((D+ 180) % 360) - 180 for D in [W,E]] + + mask = (stations[lat_key] > S) & (stations[lat_key] < N) & (stations[lon_key] < E) & (stations[lon_key] > W) + return stations[mask] + diff --git a/tools/RAiDER/llreader.py b/tools/RAiDER/llreader.py index 769b3b32c..a92b59745 100644 --- a/tools/RAiDER/llreader.py +++ b/tools/RAiDER/llreader.py @@ -183,8 +183,12 @@ def set_output_xygrid(self, dst_crs=4326): try: out_proj = CRS.from_epsg(dst_crs.replace('EPSG:', '')) - except pyproj.exceptions.CRSError: - out_proj = dst_crs + except AttributeError: + try: + out_proj = CRS.from_epsg(dst_crs) + except pyproj.exceptions.CRSError: + out_proj = dst_crs + out_snwe = transform_bbox(self.bounds(), src_crs=4326, dest_crs=out_proj) logger.debug(f"Output SNWE: {out_snwe}") @@ -230,7 +234,7 @@ def readZ(self): _, _ = download_dem( self._bounding_box, writeDEM=True, - outName=demFile, + demName=demFile, ) ## interpolate the DEM to the query points @@ -303,7 +307,7 @@ def readZ(self): _, _ = download_dem( self._bounding_box, writeDEM=True, - outName=demFile, + demName=demFile, ) z_out = interpolateDEM(demFile, self.readLL()) @@ -331,6 +335,10 @@ def __init__(self, filename, is_dem=False): self._is_dem = is_dem _, self._proj, self._geotransform = rio_stats(filename) self._type = 'geocoded_file' + try: + self.crs = self.p['crs'] + except KeyError: + self.crs = None def readLL(self): @@ -354,7 +362,7 @@ def readZ(self): demFile = self._filename if self._is_dem else 'GLO30_fullres_dem.tif' bbox = self._bounding_box - _, _ = download_dem(bbox, writeDEM=True, outName=demFile) + _, _ = download_dem(bbox, writeDEM=True, demName=demFile) z_out = interpolateDEM(demFile, self.readLL()) return z_out @@ -375,7 +383,6 @@ def get_extent(self): W, E = ds.longitude.min().item(), ds.longitude.max().item() return [S, N, W, E] - ## untested def readLL(self): with xarray.open_dataset(self.path) as ds: @@ -424,7 +431,5 @@ def bounds_from_csv(station_file): and "Lon" columns, which should be EPSG: 4326 projection (i.e WGS84) ''' stats = pd.read_csv(station_file).drop_duplicates(subset=["Lat", "Lon"]) - if 'Hgt_m' in stats.columns: - use_csv_heights = True snwe = [stats['Lat'].min(), stats['Lat'].max(), stats['Lon'].min(), stats['Lon'].max()] return snwe diff --git a/tools/RAiDER/losreader.py b/tools/RAiDER/losreader.py index c11b9a84c..426305f0b 100644 --- a/tools/RAiDER/losreader.py +++ b/tools/RAiDER/losreader.py @@ -9,6 +9,7 @@ import os import datetime import shelve +from abc import ABC from typing import Union from pathlib import PosixPath @@ -17,8 +18,10 @@ import xml.etree.ElementTree as ET except ImportError: ET = None - -from abc import ABC +try: + import isce3.ext.isce3 as isce +except ImportError: + isce = None from RAiDER.constants import _ZREF from RAiDER.utilFcns import ( @@ -178,8 +181,12 @@ class Raytracing(LOS): >>> from RAiDER.losreader import Raytracing >>> import numpy as np """ + def __init__(self, filename=None, los_convention='isce', time=None, look_dir = 'right', pad=600): '''read in and parse a statevector file''' + if isce is None: + raise ImportError(f'isce3 is required for this class. Use conda to install isce3`') + super().__init__() self._ray_trace = True self._file = filename @@ -191,7 +198,6 @@ def __init__(self, filename=None, los_convention='isce', time=None, look_dir = ' raise NotImplementedError() # ISCE3 data structures - import isce3.ext.isce3 as isce if self._time is not None: # __call__ called in checkArgs; keep for modularity self._orbit = get_orbit(self._file, self._time, pad=pad) @@ -231,14 +237,15 @@ def getLookVectors(self, ht, llh, xyz, yy): ''' Calculate look vectors for raytracing ''' + if isce is None: + raise ImportError(f'isce3 is required for this method. Use conda to install isce3`') + # TODO - Modify when isce3 vectorization is available los = np.full(yy.shape + (3,), np.nan) llh = llh.copy() llh[0] = np.deg2rad(llh[0]) llh[1] = np.deg2rad(llh[1]) - import isce3.ext.isce3 as isce - for ii in range(yy.shape[0]): for jj in range(yy.shape[1]): inp = np.array([llh[0][ii, jj], llh[1][ii, jj], ht]) @@ -373,6 +380,9 @@ def filter_ESA_orbit_file_p(path: str) -> bool: raise ValueError( 'get_sv: I cannot parse the statevector file {}'.format(los_file) ) + except: + raise ValueError('get_sv: I cannot parse the statevector file {}'.format(los_file)) + if ref_time: idx = cut_times(svs[0], ref_time, pad=pad) @@ -595,6 +605,9 @@ def state_to_los(svs, llh_targets): >>> svs = losr.read_ESA_Orbit_file(esa_orbit_file) >>> LOS = losr.state_to_los(*svs, [lats, lons, heights], xyz) ''' + if isce is None: + raise ImportError(f'isce3 is required for this function. Use conda to install isce3`') + # check the inputs if np.min(svs.shape) < 4: raise RuntimeError( @@ -603,7 +616,6 @@ def state_to_los(svs, llh_targets): ) # Convert svs to isce3 orbit - import isce3.ext.isce3 as isce orb = isce.core.Orbit([ isce.core.StateVector( isce.core.DateTime(row[0]), @@ -657,6 +669,8 @@ def get_radar_pos(llh, orb): los: ndarray - Satellite incidence angle sr: ndarray - Slant range in meters ''' + if isce is None: + raise ImportError(f'isce3 is required for this function. Use conda to install isce3`') num_iteration = 30 residual_threshold = 1.0e-7 @@ -667,7 +681,6 @@ def get_radar_pos(llh, orb): ) # Get some isce3 constants for this inversion - import isce3.ext.isce3 as isce # TODO - Assuming right-looking for now elp = isce.core.Ellipsoid() dop = isce.core.LUT2d() @@ -760,9 +773,10 @@ def get_orbit(orbit_file: Union[list, str], requested time (should be about 600 seconds) ''' - # First load the state vectors into an isce orbit - import isce3.ext.isce3 as isce + if isce is None: + raise ImportError(f'isce3 is required for this function. Use conda to install isce3`') + # First load the state vectors into an isce orbit svs = np.stack(get_sv(orbit_file, ref_time, pad), axis=-1) sv_objs = [] # format for ISCE diff --git a/tools/RAiDER/models/credentials.py b/tools/RAiDER/models/credentials.py index e35d7883b..85f8eb2d7 100644 --- a/tools/RAiDER/models/credentials.py +++ b/tools/RAiDER/models/credentials.py @@ -1,13 +1,13 @@ ''' -API credential information and help url for downloading weather model data - saved in a hidden file in home directory +API credential information and help url for downloading weather model data + saved in a hidden file in home directory api filename weather models UID KEY URL _________________________________________________________________________________ cdsapirc ERA5, ERA5T uid key https://cds.climate.copernicus.eu/api/v2 -ecmwfapirc ERAI, HRES email key https://api.ecmwf.int/v1 -netrc GMAO, MERRA2 username password urs.earthdata.nasa.gov - HRRR [public access] +ecmwfapirc HRES email key https://api.ecmwf.int/v1 +netrc GMAO, MERRA2 username password urs.earthdata.nasa.gov + HRRR [public access] ''' import os @@ -17,7 +17,6 @@ # Filename for the hidden file per model API_FILENAME = {'ERA5' : 'cdsapirc', 'ERA5T' : 'cdsapirc', - 'ERAI' : 'ecmwfapirc', 'HRES' : 'ecmwfapirc', 'GMAO' : 'netrc', 'HRRR' : None @@ -69,14 +68,14 @@ def _check_envs(model): key = os.getenv('RAIDER_ECMWF_ERA5_API_KEY') host = API_URLS['cdsapirc'] - elif model in ('HRES'): - uid = os.getenv('RAIDER_HRES_EMAIL') + elif model in ('HRES',): + uid = os.getenv('RAIDER_HRES_EMAIL') key = os.getenv('RAIDER_HRES_API_KEY') host = os.getenv('RAIDER_HRES_URL') if host is None: host = API_URLS['ecmwfapirc'] - elif model in ('GMAO'): + elif model in ('GMAO',): uid = os.getenv('EARTHDATA_USERNAME') # same as in DockerizedTopsApp key = os.getenv('EARTHDATA_PASSWORD') host = API_URLS['netrc'] @@ -107,19 +106,19 @@ def check_api(model: str, hidden_ext = '_' if system()=="Windows" else '.' # skip below if model is HRRR as it does not need API - if api_filename: + if api_filename: # Check if the credential api file exists api_filename_path = Path(output_dir) / (hidden_ext + api_filename) api_filename_path = api_filename_path.expanduser() - # if update flag is on, overwrite existing file + # if update flag is on, overwrite existing file if update_flag is True: api_filename_path.unlink(missing_ok=True) - + # Check if API_RC file already exists if api_filename_path.exists(): return None - + # if it does not exist, put UID/KEY inserted, create it elif not api_filename_path.exists() and UID and KEY: # Create file with inputs, do it only once @@ -145,3 +144,8 @@ def check_api(model: str, f'{api_filename_path}, API ENVIRONMENTALS' f' and API UID and KEY, do not exist !!' f'\nGet API info from ' + '\033[1m' f'{help_url}' + '\033[0m, and add it!') + + +def setup_from_env(): + for model in API_FILENAME.keys(): + check_api(model) diff --git a/tools/RAiDER/models/customExceptions.py b/tools/RAiDER/models/customExceptions.py new file mode 100644 index 000000000..3bc9b69a9 --- /dev/null +++ b/tools/RAiDER/models/customExceptions.py @@ -0,0 +1,66 @@ +class DatetimeFailed(Exception): + def __init__(self, model, time): + msg = f"Weather model {model} failed to download for datetime {time}" + super().__init__(msg) + + +class DatetimeNotAvailable(Exception): + def __init__(self, model, time): + msg = f"Weather model {model} was not found for datetime {time}" + super().__init__(msg) + + +class DatetimeOutsideRange(Exception): + def __init__(self, model, time): + msg = f"Time {time} is outside the available date range for weather model {model}" + super().__init__(msg) + + +class ExistingWeatherModelTooSmall(Exception): + def __init__(self): + msg = 'The weather model passed does not cover all of the input ' \ + 'points; you may need to download a larger area.' + super().__init__(msg) + + +class TryToKeepGoingError(Exception): + def __init__(self, date=None): + if date is not None: + msg = 'The weather model does not exist for date {date}, so I will try to use the closest available date.' + else: + msg = 'I will try to keep going' + super().__init__(msg) + +class CriticalError(Exception): + def __init__(self): + msg = 'I have experienced a critical error, please take a look at the log files' + super().__init__(msg) + +class WrongNumberOfFiles(Exception): + def __init__(self, Nexp, Navail): + msg = 'The number of files downloaded does not match the requested, ' + 'I expected {} and got {}, aborting'.format(Nexp, Navail) + super().__init__(msg) + + +class NoWeatherModelData(Exception): + def __init__(self, custom_msg=None): + if custom_msg is None: + msg = 'No weather model files were available to download, aborting' + else: + msg = custom_msg + super().__init__(msg) + + +class NoStationDataFoundError(Exception): + def __init__(self, station_list=None, years=None): + if (station_list is None) and (years is None): + msg = 'No GNSS station data was found' + elif (years is None): + msg = 'No data was found for GNSS stations {}'.format(station_list) + elif station_list is None: + msg = 'No data was found for years {}'.format(years) + else: + msg = 'No data was found for GNSS stations {} and years {}'.format(station_list, years) + + super().__init__(msg) diff --git a/tools/RAiDER/models/weatherModel.py b/tools/RAiDER/models/weatherModel.py index 267eacb65..37165f07f 100755 --- a/tools/RAiDER/models/weatherModel.py +++ b/tools/RAiDER/models/weatherModel.py @@ -18,6 +18,7 @@ from RAiDER.interpolator import fillna3D from RAiDER.logger import logger from RAiDER.models import plotWeather as plots, weatherModel +from RAiDER.models.customExceptions import * from RAiDER.utilFcns import ( robmax, robmin, write2NETCDF4core, calcgeoh, transform_coords, clip_bbox ) @@ -31,7 +32,6 @@ 'HRRR-AK': 3, } - class WeatherModel(ABC): ''' Implement a generic weather model for getting estimated SAR delays @@ -57,7 +57,7 @@ def __init__(self): self._classname = None self._dataset = None - self._name = None + self._Name = None self._wmLoc = None self._model_level_type = 'ml' @@ -160,12 +160,9 @@ def fetch(self, out, time): # write the error raised by the weather model API to the log try: self._fetch(out) - err = False - except Exception as E: - err = E - - return err + logger.exception(E) + raise @abstractmethod @@ -308,22 +305,17 @@ def checkTime(self, time): self.Model(), self._valid_range[0].date(), end_time ) - msg = f"Weather model {self.Model()} is not available at: {time}" - if time < self._valid_range[0]: - logger.error(msg) - raise RuntimeError(msg) + raise DatetimeOutsideRange(self.Model(), time) if self._valid_range[1] is not None: if self._valid_range[1] == 'Present': pass elif self._valid_range[1] < time: - logger.error(msg) - raise RuntimeError(msg) + raise DatetimeOutsideRange(self.Model(), time) if time > datetime.datetime.utcnow() - self._lag_time: - logger.error(msg) - raise RuntimeError(msg) + raise DatetimeOutsideRange(self.Model(), time) def setLevelType(self, levelType): @@ -748,6 +740,7 @@ def write(self): "datetime": datetime.datetime.strftime(self._time, "%Y_%m_%dT%H_%M_%S"), 'date_created': datetime.datetime.now().strftime("%Y_%m_%dT%H_%M_%S"), 'title': 'Weather model data and delay calculations', + 'model_name': self._Name } diff --git a/tools/RAiDER/processWM.py b/tools/RAiDER/processWM.py index b1b034061..2b1bdc742 100755 --- a/tools/RAiDER/processWM.py +++ b/tools/RAiDER/processWM.py @@ -15,7 +15,7 @@ from RAiDER.logger import logger from RAiDER.utilFcns import getTimeFromFile from RAiDER.models.weatherModel import make_raw_weather_data_filename, checkContainment_raw - +from RAiDER.models.customExceptions import * def prepareWeatherModel( weather_model, @@ -67,10 +67,11 @@ def prepareWeatherModel( # if no weather model files supplied, check the standard location else: - E = weather_model.fetch(path_wm_raw, time) - if E: - logger.warning (E) - raise RuntimeError + os.makedirs(os.path.dirname(path_wm_raw), exist_ok=True) + try: + weather_model.fetch(path_wm_raw, time) + except DatetimeOutsideRange: + raise TryToKeepGoingError # If only downloading, exit now if download_only: @@ -89,11 +90,9 @@ def prepareWeatherModel( ) containment = weather_model.checkContainment(ll_bounds) - if not containment and weather_model.Model() in 'GMAO ERA5 ERA5T HRES'.split(): - msg = 'The weather model passed does not cover all of the input ' \ - 'points; you may need to download a larger area.' - logger.error(msg) - raise RuntimeError(msg) + if not containment and weather_model.Model() not in 'HRRR'.split(): + raise ExistingWeatherModelTooSmall + return f # Logging some basic info @@ -131,17 +130,14 @@ def prepareWeatherModel( except Exception as e: logger.exception("Unable to save weathermodel to file") logger.exception(e) - raise RuntimeError("Unable to save weathermodel to file") + raise CriticalError finally: wm = weather_model.Model() del weather_model - if not containment and wm in 'GMAO ERA5 ERA5T HRES'.split(): - msg = 'The weather model passed does not cover all of the input ' \ - 'points; you may need to download a larger area.' - logger.error(msg) - raise RuntimeError(msg) + if not containment and wm not in 'HRRR'.split(): + raise ExistingWeatherModelTooSmall else: return f diff --git a/tools/RAiDER/s1_azimuth_timing.py b/tools/RAiDER/s1_azimuth_timing.py index 70a44ee24..3fc20d353 100644 --- a/tools/RAiDER/s1_azimuth_timing.py +++ b/tools/RAiDER/s1_azimuth_timing.py @@ -2,13 +2,17 @@ import warnings import asf_search as asf -import hyp3lib.get_orb -import isce3.ext.isce3 as isce import numpy as np import pandas as pd from shapely.geometry import Point -from .losreader import get_orbit as get_isce_orbit +try: + import isce3.ext.isce3 as isce +except ImportError: + isce = None + +from RAiDER.losreader import get_orbit as get_isce_orbit +from RAiDER.s1_orbits import get_orbits_from_slc_ids def _asf_query(point: Point, @@ -76,7 +80,7 @@ def get_slc_id_from_point_and_time(lon: float, def get_azimuth_time_grid(lon_mesh: np.ndarray, lat_mesh: np.ndarray, hgt_mesh: np.ndarray, - orb: isce.core.Orbit) -> np.ndarray: + orb: 'isce.core.Orbit') -> np.ndarray: ''' Source: https://github.com/dbekaert/RAiDER/blob/dev/tools/RAiDER/losreader.py#L601C1-L674C22 @@ -84,6 +88,8 @@ def get_azimuth_time_grid(lon_mesh: np.ndarray, Technically, this is "sensor neutral" since it uses an orb object. ''' + if isce is None: + raise ImportError(f'isce3 is required for this function. Use conda to install isce3`') num_iteration = 100 residual_threshold = 1.0e-7 @@ -176,10 +182,13 @@ def get_s1_azimuth_time_grid(lon: np.ndarray, np.datetime64('NaT'), dtype='datetime64[ms]') return az_arr - orb_files = list(map(lambda slc_id: hyp3lib.get_orb.downloadSentinelOrbitFile(slc_id)[0], slc_ids)) - orb = get_isce_orbit(orb_files, dt, pad=600) + orb_files = get_orbits_from_slc_ids(slc_ids) + orb_files = [str(of) for of in orb_files] + + orb = get_isce_orbit(orb_files, dt, pad=600) az_arr = get_azimuth_time_grid(lon_mesh, lat_mesh, hgt_mesh, orb) + return az_arr @@ -330,6 +339,8 @@ def get_inverse_weights_for_dates(azimuth_time_array: np.ndarray, n_dates = len(dates) if n_unique_dates != n_dates: raise ValueError('Dates provided must be unique') + if n_dates == 0: + raise ValueError('No dates provided') if not all([isinstance(date, datetime.datetime) for date in dates]): raise TypeError('dates must be all datetimes') diff --git a/tools/RAiDER/s1_orbits.py b/tools/RAiDER/s1_orbits.py new file mode 100644 index 000000000..3aaa278a0 --- /dev/null +++ b/tools/RAiDER/s1_orbits.py @@ -0,0 +1,68 @@ +import netrc +import os +import re +from pathlib import Path +from platform import system +from typing import List, Optional, Tuple + +import eof.download + + +ESA_CDSE_HOST = 'dataspace.copernicus.eu' + + +def _netrc_path() -> Path: + netrc_name = '_netrc' if system().lower() == 'windows' else '.netrc' + return Path.home() / netrc_name + + +def ensure_orbit_credentials() -> Optional[int]: + """Ensure credentials exist for ESA's CDSE to download orbits + + This method will prefer to use CDSE credentials from your `~/.netrc` file if they exist, + otherwise will look for ESA_USERNAME and ESA_PASSWORD environment variables and + update or create your `~/.netrc` file. + + Returns `None` if the `~/.netrc` file did not need to be updated and the number of characters written if it did. + """ + netrc_file = _netrc_path() + + # netrc needs a netrc file; if missing create an empty one. + if not netrc_file.exists(): + netrc_file.touch() + + netrc_credentials = netrc.netrc(netrc_file) + if ESA_CDSE_HOST in netrc_credentials.hosts: + return + + username = os.environ.get('ESA_USERNAME') + password = os.environ.get('ESA_PASSWORD') + if username is None and password is None: + raise ValueError('Credentials are required for fetching orbit data from dataspace.copernicus.eu!\n' + 'Either add your credentials to ~/.netrc or set the ESA_USERNAME and ESA_PASSWORD ' + 'environment variables.') + + netrc_credentials.hosts[ESA_CDSE_HOST] = (username, None, password) + return netrc_file.write_text(str(netrc_credentials)) + + +def get_orbits_from_slc_ids(slc_ids: List[str], directory=Path.cwd()) -> List[Path]: + """Download all orbit files for a set of SLCs + + This method will ensure that the downloaded orbit files cover the entire acquisition start->stop time + + Returns a list of orbit file paths + """ + _ = ensure_orbit_credentials() + + missions = {slc_id[0:3] for slc_id in slc_ids} + start_times = {re.split(r'_+', slc_id)[4] for slc_id in slc_ids} + stop_times = {re.split(r'_+', slc_id)[5] for slc_id in slc_ids} + + orb_files = eof.download.download_eofs( + sorted(list(start_times | stop_times)), + sorted(list(missions)), + save_dir=str(directory) + ) + + return orb_files diff --git a/tools/RAiDER/utilFcns.py b/tools/RAiDER/utilFcns.py index d9f2c23cb..dce21673b 100755 --- a/tools/RAiDER/utilFcns.py +++ b/tools/RAiDER/utilFcns.py @@ -134,14 +134,7 @@ def rio_profile(fname): with rasterio.open(fname) as src: profile = src.profile - # if 'S1-GUNW' in fname: - # profile['length'] = profile['width'] - # profile['width'] = profile['height'] - if profile["crs"] is None: - raise AttributeError( - f"{fname} does not contain geotransform information" - ) return profile @@ -150,9 +143,6 @@ def rio_extents(profile): gt = profile["transform"].to_gdal() xSize = profile["width"] ySize = profile["height"] - - if profile["crs"] is None or not gt: - raise AttributeError('Profile does not contain geotransform information') W, E = gt[0], gt[0] + (xSize - 1) * gt[1] + (ySize - 1) * gt[2] N, S = gt[3], gt[3] + (xSize - 1) * gt[4] + (ySize - 1) * gt[5] return S, N, W, E @@ -280,12 +270,15 @@ def writeArrayToRaster(array, filename, noDataValue=0., fmt='ENVI', proj=None, g # Geotransform trans = None if gt is not None: - trans = rasterio.Affine.from_gdal(*gt) + try: + trans = rasterio.Affine.from_gdal(*gt) + except TypeError: + trans = gt ## cant write netcdfs with rasterio in a simple way if fmt == 'nc': fmt = 'GTiff' - filename = filename.replace('.nc', '.GTiff') + filename = filename.replace('.nc', '.tif') with rasterio.open(filename, mode="w", count=1, width=array_shp[1], height=array_shp[0], @@ -296,17 +289,6 @@ def writeArrayToRaster(array, filename, noDataValue=0., fmt='ENVI', proj=None, g return -def writeArrayToFile(lats, lons, array, filename, noDataValue=-9999): - ''' - Write a single-dim array of values to a file - ''' - array[np.isnan(array)] = noDataValue - with open(filename, 'w') as f: - f.write('Lat,Lon,Hgt_m\n') - for lat, lon, height in zip(lats, lons, array): - f.write('{},{},{}\n'.format(lat, lon, height)) - - def round_date(date, precision): # First try rounding up # Timedelta since the beginning of time @@ -342,20 +324,14 @@ def robmin(a): ''' Get the minimum of an array, accounting for empty lists ''' - try: - return np.nanmin(a) - except ValueError: - return 'N/A' + return np.nanmin(a) def robmax(a): ''' Get the minimum of an array, accounting for empty lists ''' - try: - return np.nanmax(a) - except ValueError: - return 'N/A' + return np.nanmax(a) def _get_g_ll(lats): @@ -427,25 +403,6 @@ def padLower(invar): return np.concatenate((new_var[:, :, np.newaxis], invar), axis=2) -def checkShapes(los, lats, lons, hts): - ''' - Make sure that by the time the code reaches here, we have a - consistent set of line-of-sight and position data. - ''' - from RAiDER.losreader import Zenith - test1 = hts.shape == lats.shape == lons.shape - try: - test2 = los.shape[:-1] == hts.shape - except AttributeError: - test2 = los is Zenith - - if not test1 and test2: - raise ValueError( - 'I need lats, lons, heights, and los to all be the same shape. ' + - 'lats had shape {}, lons had shape {}, '.format(lats.shape, lons.shape) + - 'heights had shape {}, and los was not Zenith'.format(hts.shape)) - - def round_time(dt, roundTo=60): ''' Round a datetime object to any time lapse in seconds @@ -471,11 +428,7 @@ def writeDelays(aoi, wetDelay, hydroDelay, # Do different things, depending on the type of input if aoi.type() == 'station_file': - #TODO: why is this a try/except? - try: - df = pd.read_csv(aoi._filename).drop_duplicates(subset=["Lat", "Lon"]) - except ValueError: - df = pd.read_csv(aoi._filename).drop_duplicates(subset=["Lat", "Lon"]) + df = pd.read_csv(aoi._filename).drop_duplicates(subset=["Lat", "Lon"]) df['wetDelay'] = wetDelay df['hydroDelay'] = hydroDelay @@ -510,11 +463,9 @@ def getTimeFromFile(filename): ''' fmt = '%Y_%m_%d_T%H_%M_%S' p = re.compile(r'\d{4}_\d{2}_\d{2}_T\d{2}_\d{2}_\d{2}') - try: - out = p.search(filename).group() - return datetime.strptime(out, fmt) - except BaseException: # TODO: Which error(s)? - raise RuntimeError('The filename for {} does not include a datetime in the correct format'.format(filename)) + out = p.search(filename).group() + return datetime.strptime(out, fmt) + # Part of the following UTM and WGS84 converter is borrowed from https://gist.github.com/twpayne/4409500