diff --git a/.circleci/config.yml b/.circleci/config.yml index 40a247c06..21ee62ba7 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -43,7 +43,7 @@ jobs: eval "$($HOME/bin/micromamba shell hook -s posix)" micromamba activate RAiDER python -m pip install . - python -c "import RAiDER; from RAiDER.delay import main" + python -c "import RAiDER; from RAiDER.delay import tropo_delay" python -c "import RAiDER; from RAiDER.interpolator import interp_along_axis" - run: name: Run unit tests diff --git a/CHANGELOG.md b/CHANGELOG.md index 185a25b9b..3919476ce 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,25 @@ 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.3.0] +RAiDER package was refactored to expose the main functionality as a Python library, including the `prepareWeatherModel` +and `tropo_delay` functions, as well as anciliarry functions needed for defining AOIs, look vectors, etc. + +### New/Updated Features ++ Python library access to main functions for accessing weather model data and calculating delays ++ Slant delay calculation through projection is supported for cubes with orbit files ++ Upgrade dem-stitcher to [`>=2.3.1`](https://github.com/ACCESS-Cloud-Based-InSAR/dem-stitcher/blob/dev/CHANGELOG.md#231) so that the updated urls for the GLO-30 DEM are used. ++ `raider.py ++calcDelaysGUNW GUNWFILE` is enabled as a placeholder only. ++ Upgraded ISCE3 to `>=v0.9.0` to fix a conda build issue as described in [#425](https://github.com/dbekaert/RAiDER/issues/425) ++ Allow user to specify --download_only or download_only=True in the configure file ++ Added documentation for the Python library interface. ++ Added some unit tests. ++ Fixed some bugs and tweaked the CLI. ++ Added unit tests, docstrings, initial API reference ++ __main__ file to allow calls to different functionality. `raider.py ++process downloadGNSS ...` can now perform the functionality of `raiderDownloadGNSS.py ... + + + ## [0.2.0] RAiDER package was refactored to use a configure file (yaml) to parse parameters. In addition, ocker container images diff --git a/README.md b/README.md old mode 100755 new mode 100644 index a401a581c..1a00261b8 --- a/README.md +++ b/README.md @@ -1,4 +1,5 @@ # RAiDER + Raytracing Atmospheric Delay Estimation for RADAR [![Language](https://img.shields.io/badge/python-3.7%2B-blue.svg)](https://www.python.org/) @@ -13,15 +14,18 @@ Copyright (c) 2019-2022, California Institute of Technology ("Caltech"). All rig THIS IS RESEARCH CODE PROVIDED TO YOU "AS IS" WITH NO WARRANTIES OF CORRECTNESS. USE AT YOUR OWN RISK. ## Contents -- [1. Getting Started](#1-getting-started) - - [Installing With Conda](#installing-with-conda) - - [Using the Docker Image](#using-the-docker-image) - - [Installing from Source](#installing-from-source) -- [2. Setup of third party weather model access](#2-setup-of-third-party-weather-model-access) -- [3. Running RAiDER and Documentation](#3-running-raider-and-documentation) -- [4. Citing](#4-citation) -- [5. Contributors](#5-contributors) + +1. [Getting Started](#1-getting-started) + - [Installing With Conda](#installing-with-conda) + - [Using the Docker Image](#using-the-docker-image) + - [Installing from Source](#installing-from-source) +2. [Setup of third party weather model access](#2-setup-of-third-party-weather-model-access) +3. [Running RAiDER and Documentation](#3-running-raider-and-documentation) +4. [Citing](#4-citation) +5. [Development](#5-development) + - [Contributors](#contributors) ------ + ## 1. Getting Started RAiDER has been tested on the following systems: @@ -31,6 +35,7 @@ RAiDER has been tested on the following systems: RAiDER does **not** currently run on arm64 processors on Mac. We will update this note once the build becomes available. ### Installing With Conda + RAiDER is available on [conda-forge](https://anaconda.org/conda-forge/raider). __[Conda](https://docs.conda.io/en/latest/index.html)__ is a cross-platform way to use Python that allows you to setup and use "virtual environments." These can help to keep dependencies for different sets of code separate. We recommend using [Miniforge](https://github.com/conda-forge/miniforge), a conda environment manager that uses conda-forge as its default code repo. Alternatively,see __[here](https://docs.anaconda.com/anaconda/install/)__ for help installing Anaconda and __[here](https://docs.conda.io/en/latest/miniconda.html)__ for installing Miniconda. Installing RAiDER: @@ -40,6 +45,7 @@ conda activate RAiDER ``` ### Using the Docker image + RAiDER provides a [docker container image](https://docs.docker.com/get-started/) with all the necessary dependencies pre-installed. To get the latest released version: ``` docker pull ghcr.io/dbekaert/raider:latest @@ -64,32 +70,43 @@ cd work ``` For more docker run options, see: . -### Installing from source -You can also install RAiDER directly from source. Doing so is recommended for those who would like to [contribute to the source code](https://github.com/dbekaert/RAiDER/blob/dev/CONTRIBUTING.md), which we heartily encourage! For more details on installing from source see [here](https://github.com/dbekaert/RAiDER/blob/dev/Installing_from_source.md). -``` -git clone https://github.com/dbekaert/RAiDER.git -cd RAiDER -conda create -f environment.yml -conda activate RAiDER -python -m pip install -e . -``` + ------ ## 2. Setup of third party weather model access -RAiDER has the ability to download weather models from third-parties; some of which require license agreements. See [here](WeatherModels.md) for details. + +RAiDER has the ability to download weather models from third-parties; some of which require license agreements. See [here](https://github.com/dbekaert/RAiDER/blob/dev/docs/WeatherModels.md) for details. ------ ## 3. Running RAiDER and Documentation + For detailed documentation, examples, and Jupyter notebooks see the [RAiDER-docs repository](https://github.com/dbekaert/RAiDER-docs). -We welcome contributions of other examples on how to leverage the RAiDER (see [here](https://github.com/dbekaert/RAiDER/blob/master/CONTRIBUTING.md) for instructions). +We welcome contributions of other examples on how to leverage the RAiDER (see [here](https://github.com/dbekaert/RAiDER/blob/dev/CONTRIBUTING.md) for instructions). ``` raiderDelay.py -h ``` provides a help menu and list of example commands to get started. -The RAiDER scripts are highly modulized in Python and allows for building your own processing workflow. +The RAiDER scripts are highly modularized in Python and allows for building your own processing workflow. ------ ## 4. Citation TODO ------ -## 5. Contributors +## 5. Development + +Contributions are welcome and heartily encourage! See our [contributing guide](https://github.com/dbekaert/RAiDER/blob/dev/CONTRIBUTING.md). + +### Development install +For development, we recommend installing directly from source. +``` +git clone https://github.com/dbekaert/RAiDER.git +cd RAiDER +conda create -f environment.yml +conda activate RAiDER +python -m pip install -e . +``` +For more details on installing from source see [here](https://github.com/dbekaert/RAiDER/blob/dev/docs/Installing_from_source.md). + +------ +### Contributors + * David Bekaert * Jeremy Maurer * Raymond Hogenson diff --git a/Installing_from_source.md b/docs/Installing_from_source.md similarity index 100% rename from Installing_from_source.md rename to docs/Installing_from_source.md diff --git a/WeatherModels.md b/docs/WeatherModels.md similarity index 62% rename from WeatherModels.md rename to docs/WeatherModels.md index eae7cbd5b..ddbc87eea 100644 --- a/WeatherModels.md +++ b/docs/WeatherModels.md @@ -1,30 +1,28 @@ # Accessing weather model data -RAiDER has built-in support for a number of different weather models. RAiDER provides all of the interfacing to data servers required to access data for the different weather models, although some weather models require a license agreement and accounts to be set-up. Instructions for accessing data, including license-limited data, are provided below. It is the user's responsibility to accept license agreements for whatever model is desired. +RAiDER has built-in support for a number of different weather models. RAiDER provides all the interfacing to data servers required to access data for the different weather models, although some weather models require a license agreement and accounts to be set-up. Instructions for accessing data, including license-limited data, are provided below. It is the user's responsibility to accept license agreements for whatever model is desired. In addition, RAiDER provides functionality for adding additional weather models. See the [RAiDER-docs repository](https://github.com/dbekaert/RAiDER-docs) page on how to do this. We would love to expand the suite of supported models, and welcome any contributions. Please see the contributing guidelines or reach out through an issue ticket for help. ------ -## Contents +## 1. Usage -1. [NOAA weather models (HRRR)](#noaa-weather-models-(hrrr)) -2. [ECMWF weather models (ERA5, ERA5T, ERAI, HRES)](#ecmwf-weather-models-(era5,-era5t,-erai,-hres)) -3. [NASA weather models (GMAO, MERRA2)](#nasa-weather-models-(gmao,-merra2)) +::: RAiDER.processWM.prepareWeatherModel + options: + show_root_heading: true + heading_level: 3 -#Potential download failure:# +### Potential download failure ERA-5/ERA-I products require access to the ESA Copernicus servers. GMAO and MERRA-2 products require access to the NASA Earthdata servers. If you are unable to download products, ensure that you have registered and have downloaded the public API key, and accepted/added the license/application for type of product you wish to download as detailed below. ------ -## 1. NOAA weather models (HRRR) -High-resolution rapid refresh (HRRR) weather model data products are generated by __[NOAA](https://rapidrefresh.noaa.gov/hrrr/)__ for the coninental US (CONUS) but not archived beyond three days. However a public __[archive](home.chpc.utah.edu/~u0553130/Brian_Blaylock/hrrr_FAQ.html)__ is available at the University of Utah. This archive does not require a license agreement. This model has the highest spatial resolution available in RAiDER, with a horizontal grid spacing of about 3 km, and is provided in a Lambert conformal conic projection. - -[return to table of content](#contents) - +## 2. NOAA weather models (HRRR) +High-resolution rapid refresh (HRRR) weather model data products are generated by __[NOAA](https://rapidrefresh.noaa.gov/hrrr/)__ for the coninental US (CONUS) but not archived beyond three days. However, a public __[archive](home.chpc.utah.edu/~u0553130/Brian_Blaylock/hrrr_FAQ.html)__ is available at the University of Utah. This archive does not require a license agreement. This model has the highest spatial resolution available in RAiDER, with a horizontal grid spacing of about 3 km, and is provided in a Lambert conformal conic projection. ------ -## 2. ECMWF weather models (ERA5, ERA5T, ERAI, HRES) +## 3. ECMWF weather models (ERA5, ERA5T, ERAI, HRES) The Copernicus Climate Data Store (CDS) provides access to the European Centre for Medium-Range Weather Forecasts (__[ECMWF](https://www.ecmwf.int/)__) provides a number of different weather models, including ERA5 and ERA5T reanalysis models. The ECMWF provides access to both reanalysis and real-time prediction models. You can read more information about their reanalysis models __[here](https://www.ecmwf.int/en/research/climate-reanalysis)__ and real-time model __[here](https://www.ecmwf.int/en/forecasts/datasets/catalogue-ecmwf-real-time-products)__. ECMWF models are global, with horizontal resolution of about 30 km for ERA-I, ERA-5, and ERA-5T, and 6 km for Hi-RES. All of these models come in a global projection (EPSG 4326, WGS-84). @@ -34,14 +32,19 @@ The ECMWF provides access to both reanalysis and real-time prediction models. Yo 2. Confirm your email, etc. 3. Install the public API key and client as instructed __[here](https://cds.climate.copernicus.eu/api-how-to)__: - a. Copy the URL and API key from the webpage into a file in your home directory name ~/.cdsapirc - url: https://cds.climate.copernicus.eu/api/v2 - key: your_key_here - __**Note**: the key represents the API key obtained upon the registration of CDS API, and should be replaced with the user's own information.__ - - b. Install the CDS API using pip: - pip install cdsapi - ___**Note**: this step has been included in the conda install of RAiDER, thus can be omitted if one uses the recommended conda install of RAiDER___ + a. Copy the URL and API key from the webpage into a file in your home directory name `~/.cdsapirc` + + url: https://cds.climate.copernicus.eu/api/v2 + key: your_key_here + + **Note**: the key represents the API key obtained upon the registration of CDS API, and should be replaced with the user's own information. + + b. Install the CDS API using pip + + pip install cdsapi + + **Note**: this step has been included in the conda install of RAiDER, thus can be omitted if one uses the recommended conda install of RAiDER + 4. You must accept the [license](https://cds.climate.copernicus.eu/cdsapp/#!/terms/licence-to-use-copernicus-products) for each product you wish to download. @@ -49,26 +52,27 @@ The ECMWF provides access to both reanalysis and real-time prediction models. Yo ECMWF requires a license agreement to be able to access, download, and use their products. Instructions for completing this process is below. -1. Create an account on the ECMWF servers __[here](https://accounts.ecmwf.int/auth/realms/ecmwf/protocol/openid-connect/auth?response_type=code&scope=openid%20email&client_id=apache-www&state=sBYlpcTRPhat8d6uuM9swLCxuP8&redirect_uri=https%3A%2F%2Fwww.ecmwf.int%2Foidc.cgi&nonce=RyEzBUy4m6oo_HxRQEmJxbc5jrKY4KFZd1Usgi8cpnM)__. The ERA-I model is open-access, while HRES requires a special liscence agreement. +1. Create an account on the ECMWF servers __[here](https://accounts.ecmwf.int/auth/realms/ecmwf/protocol/openid-connect/auth?response_type=code&scope=openid%20email&client_id=apache-www&state=sBYlpcTRPhat8d6uuM9swLCxuP8&redirect_uri=https%3A%2F%2Fwww.ecmwf.int%2Foidc.cgi&nonce=RyEzBUy4m6oo_HxRQEmJxbc5jrKY4KFZd1Usgi8cpnM)__. The ERA-I model is open-access, while HRES requires a special licence agreement. 2. Confirm your email, etc. 3. Install the public API key and client as instructed __[here](https://confluence.ecmwf.int/display/WEBAPI/Access+ECMWF+Public+Datasets#AccessECMWFPublicDatasets-key)__: - a. Copy the URL and API key from the webpage into a file in your home directory name ~/.ecmwfapirc - { - "url" : "https://api.ecmwf.int/v1", - "key" : your key here, - "email" : your email here - } + a. Copy the URL and API key from the webpage into a file in your home directory name `~/.ecmwfapirc` + + { + "url" : "https://api.ecmwf.int/v1", + "key" : your key here, + "email" : your email here + } - __**Note**: the email that is used to register the user account, and the key represents the API key obtained upon the registration of ECMWF API, and should be replaced with the user's own information.__ + **Note**: the email that is used to register the user account, and the key represents the API key obtained upon the registration of ECMWF API, and should be replaced with the user's own information. - b. Install the ECMWF API using pip: - ```pip install ecmwf-api-client``` - ___**Note**: this step has been included in the conda install of RAiDER, thus can be omitted if one uses the recommended conda install of RAiDER___ + b. Install the ECMWF API using pip: -[return to table of content](#contents) + pip install ecmwf-api-client` -### 3. NASA weather models (GMAO, MERRA2) + **Note**: this step has been included in the conda install of RAiDER, thus can be omitted if one uses the recommended conda install of RAiDER + +## 4. NASA weather models (GMAO, MERRA2) 1. The Global Modeling and Assimilation Office (__[GMAO](https://www.nccs.nasa.gov/services/data-collections/coupled-products/geos5-forecast#:~:text=The%20Global%20Modeling%20and%20Assimilation,near%2Dreal%2Dtime%20production.)__) at NASA generates reanalysis weather models. GMAO products can also be accessed without a license agreement through the pyDAP interface implemented in RAiDER. GMAO has a horizontal grid spacing of approximately 33 km, and its projection is EPSG code 4326 (WGS-84). @@ -85,7 +89,7 @@ Reference: __[The Modern-Era Retrospective Analysis for Research and Application login password - __**Note**: the username and password represent the user's username and password.__ + **Note**: the username and password represent the user's username and password. 4. Add the application `NASA GESDISC DATA ARCHIVE` by clicking on the `Applications->Authorized Apps` on the menu after logging into your Earthdata profile, and then scrolling down to the application `NASA GESDISC DATA ARCHIVE` to approve it. _This seems not required for GMAO for now, but recommended to do so for all OpenDAP-based weather models._ 5. Install the OpenDAP using pip: @@ -93,9 +97,6 @@ Reference: __[The Modern-Era Retrospective Analysis for Research and Application pip install pydap==3.2.1 - ___**Note**: this step has been included in the conda install of RAiDER, thus can be omitted if one uses the recommended conda install of RAiDER___ + **Note**: this step has been included in the conda install of RAiDER, thus can be omitted if one uses the recommended conda install of RAiDER - ___**Note**: PyDAP v3.2.1 is required for now (thus specified in the above pip install command) because the latest v3.2.2 (as of now) has a known [bug](https://colab.research.google.com/drive/1f_ss1Oa3VzgAOd_p8sgekdnLVE5NW6s5) in accessing and slicing the GMAO data. This bug is expected to be fixed in newer versions of PyDAP.___ - -[return to table of content](#contents) - + **Note**: PyDAP v3.2.1 is required for now (thus specified in the above pip install command) because the latest v3.2.2 (as of now) has a known [bug](https://colab.research.google.com/drive/1f_ss1Oa3VzgAOd_p8sgekdnLVE5NW6s5) in accessing and slicing the GMAO data. This bug is expected to be fixed in newer versions of PyDAP. diff --git a/docs/index.md b/docs/index.md new file mode 120000 index 000000000..32d46ee88 --- /dev/null +++ b/docs/index.md @@ -0,0 +1 @@ +../README.md \ No newline at end of file diff --git a/docs/macros.py b/docs/macros.py new file mode 100644 index 000000000..dc8fd6674 --- /dev/null +++ b/docs/macros.py @@ -0,0 +1,17 @@ +import requests + +import RAiDER + + +def define_env(env): + """Macros Hook""" + + @env.macro + def raider_version(): + return RAiDER.__version__ + + @env.macro + def get_content(url): + response = requests.get(url) + response.raise_for_status() + return response.content.decode() diff --git a/docs/reference.md b/docs/reference.md new file mode 100644 index 000000000..da94b41f2 --- /dev/null +++ b/docs/reference.md @@ -0,0 +1,5 @@ +# RAiDER *v{{ raider_version() }}* API Reference + +::: RAiDER + options: + show_submodules: true diff --git a/docs/reference_td.md b/docs/reference_td.md new file mode 100644 index 000000000..8832496ca --- /dev/null +++ b/docs/reference_td.md @@ -0,0 +1,3 @@ +***`tropo_delay`*** + +::: RAiDER.delay \ No newline at end of file diff --git a/docs/tutorials.md b/docs/tutorials.md new file mode 100644 index 000000000..a56121378 --- /dev/null +++ b/docs/tutorials.md @@ -0,0 +1,3 @@ +### Tutorials + +{{ get_content('https://raw.githubusercontent.com/dbekaert/RAiDER-docs/main/README.md') }} diff --git a/environment.yml b/environment.yml index 8f5e10837..b018451be 100644 --- a/environment.yml +++ b/environment.yml @@ -17,12 +17,12 @@ dependencies: - cxx-compiler - cython - dask - - dem_stitcher>=2.3.0 + - dem_stitcher>=2.3.1 - ecmwf-api-client - h5py - herbie-data - h5netcdf - - isce3==0.8.0=*_2 + - isce3>=0.9.0 - lxml - matplotlib - netcdf4 @@ -31,11 +31,14 @@ dependencies: - progressbar - pydap>3.2.2 - pyproj>=2.2.0 + - pyyaml - rasterio>=1.3.0 + - requests - s3fs - scipy - shapely - sysroot_linux-64 + - tqdm - xarray # For packaging and testing - autopep8 @@ -44,7 +47,13 @@ dependencies: - pytest-timeout - pytest-console-scripts - setuptools_scm >=6.2 - + # For docs website + - mkdocs + - mkdocstrings + - mkdocstrings-python + - mkdocs-macros-plugin + - mkdocs-material + - mkdocs-material-extensions # For RAiDER-docs - jupyterlab - jupyter_contrib_nbextensions diff --git a/mkdocs.yml b/mkdocs.yml new file mode 100644 index 000000000..1d799865b --- /dev/null +++ b/mkdocs.yml @@ -0,0 +1,29 @@ +site_name: RAiDER API Documentation +site_description: +repo_url: https://github.com/dbekaert/RAiDER +repo_name: RAiDER + +nav: + - Home page: index.md + - Tutorials: + - Overview: tutorials.md + - RAiDER tutorial: https://nbviewer.org/github/dbekaert/RAiDER-docs/blob/main/notebooks/RAiDER_tutorial/RAiDER_tutorial.ipynb" target="_blank + - raiderStats tutorial: https://nbviewer.org/github/dbekaert/RAiDER-docs/blob/main/notebooks/raiderStats/raiderStats_tutorial.ipynb" target="_blank + - Downloading GNSS tropospheric delays: https://nbviewer.org/github/dbekaert/RAiDER-docs/blob/main/notebooks/raiderDownloadGNSS/raiderDownloadGNSS_tutorial.ipynb" target="_blank + - GNSS delay manipulation with Pandas: https://nbviewer.org/github/dbekaert/RAiDER-docs/blob/main/notebooks/Pandas_tutorial/Pandas_tutorial.ipynb" target="_blank + - Weather Model Access: WeatherModels.md + - Delay Calculation: reference_td.md + - API Reference: reference.md + + +theme: + name: material + +plugins: + - search + - macros: + module_name: docs/macros + - mkdocstrings: + handlers: + python: + inherited_members: true diff --git a/pyproject.toml b/pyproject.toml index 170f9830e..19b81705d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -15,7 +15,8 @@ authors = [ {name="David Bekaert", email="david.bekaert@jpl.nasa.gov"}, {name="Jeremy Maurer", email="maurer.jeremy@gmail.com"}, {name="Piyush Agram", email="piyush@descarteslabs.com"}, - {name="Simran Sangha", email="simran.s.sangha@jpl.nasa.gov"} + {name="Simran Sangha", email="simran.s.sangha@jpl.nasa.gov"}, + {name="Brett Buzzanga", email="buzzanga@jpl.nasa.gov"}, ] description = "Raytracing Atmospheric Delay Estimation for RADAR" readme = "README.md" @@ -40,11 +41,13 @@ repository = "https://github.com/dbekaert/RAiDER" "Bug Tracker" = "https://github.com/dbekaert/RAiDER/issues" [project.scripts] -"generateGACOSVRT.py" = "RAiDER.models.generateGACOSVRT:main" -"raiderDownloadGNSS.py" = "RAiDER.gnss.downloadGNSSDelays:main" -"raider.py" = "RAiDER.cli.raider:main" -"prepFromAria.py" = "RAiDER.cli.prepFromAria:main" +"raider.py" = "RAiDER.cli.__main__:main" +"calcDelays.py" = "RAiDER.cli.__main__:calcDelays" +"calcDelaysGUNW.py" = "RAiDER.cli.__main__:calcDelaysGUNW" +"raiderDownloadGNSS.py" = "RAiDER.cli.__main__:downloadGNSS" +"downloadGNSS.py" = "RAiDER.cli.__main__:downloadGNSS" "raiderStats.py" = "RAiDER.cli.statsPlot:main" +"generateGACOSVRT.py" = "RAiDER.models.generateGACOSVRT:main" [tool.setuptools] zip-safe = false diff --git a/test/_checkArgs.py b/test/_checkArgs.py deleted file mode 100644 index 8be5b17e0..000000000 --- a/test/_checkArgs.py +++ /dev/null @@ -1,442 +0,0 @@ -import datetime -import os -import pytest - -import multiprocessing as mp -import numpy as np -import pandas as pd - -from test import TEST_DIR, pushd - -# import RAiDER.runProgram -from RAiDER.cli.raider import ( - parseCMD, read_template_file, create_parser, read_template_file, -) -from RAiDER.checkArgs import checkArgs, makeDelayFileNames -from RAiDER.constants import _ZREF -from RAiDER.losreader import Zenith, Conventional, Raytracing - - -SCENARIO_1 = os.path.join(TEST_DIR, "scenario_1") -SCENARIO_2 = os.path.join(TEST_DIR, "scenario_2") - - -def isWriteable(dirpath): - '''Test whether a directory is writeable''' - try: - filehandle = open(os.path.join(dirpath, 'tmp.txt'), 'w') - filehandle.close() - return True - except IOError: - return False - - -@pytest.fixture -def parsed_args(tmp_path): - parser = create_parser() - args = parser.parse_args([ - '--date', '20200103', - '--time', '23:00:00', - # '--latlon', 'latfile.dat', 'lonfile.dat', - '--bbox', '-1', '1', '-1', '1', - '--model', 'ERA5', - '--outformat', 'hdf5' - ]) - return args, parser - - -def test_checkArgs_outfmt_1(parsed_args): - '''Test that passing height levels with hdf5 outformat works''' - args, p = parsed_args - args.outformat = 'hdf5' - args.heightlvs = [10, 100, 1000] - checkArgs(args, p) - assert True - - -def test_checkArgs_outfmt_2(parsed_args): - '''Test that passing a raster format with height levels throws an error''' - args, p = parsed_args - args.heightlvs = [10, 100, 1000] - args.outformat = 'envi' - with pytest.raises(ValueError): - checkArgs(args, p) - - -def test_checkArgs_outfmt_3(parsed_args): - '''Test that passing a raster format with height levels throws an error''' - args, p = parsed_args - args.query_area = os.path.join(SCENARIO_2, 'stations.csv') - argDict = checkArgs(args, p) - assert argDict['flag'] == 'station_file' - - -def test_checkArgs_outfmt_4(parsed_args): - '''Test that passing a raster format with height levels throws an error''' - args, p = parsed_args - args.query_area = [os.path.join(SCENARIO_1, 'geom', 'lat.dat'), os.path.join(SCENARIO_1, 'geom', 'lat.dat')] - argDict = checkArgs(args, p) - assert argDict['flag'] == 'files' - - -def test_checkArgs_outfmt_5(parsed_args): - '''Test that passing a raster format with height levels throws an error''' - args, p = parsed_args - args.query_area = os.path.join(SCENARIO_2, 'stations.csv') - argDict = checkArgs(args, p) - assert pd.read_csv(argDict['wetFilenames'][0]).shape == (8, 4) - - -def test_checkArgs_outloc_1(parsed_args): - '''Test that the default output and weather model directories are correct''' - args, p = parsed_args - argDict = checkArgs(args, p) - out = argDict['out'] - wmLoc = argDict['wmLoc'] - assert os.path.abspath(out) == os.getcwd() - assert os.path.abspath(wmLoc) == os.path.join(os.getcwd(), 'weather_files') - - -def test_checkArgs_outloc_2(parsed_args, tmp_path): - '''Tests that the correct output location gets assigned when provided''' - with pushd(tmp_path): - args, p = parsed_args - args.out = tmp_path - argDict = checkArgs(args, p) - out = argDict['out'] - assert out == tmp_path - - -def test_checkArgs_outloc_2b(parsed_args, tmp_path): - ''' Tests that the weather model directory gets passed through by itself''' - with pushd(tmp_path): - args, p = parsed_args - args.out = tmp_path - args.wmLoc = 'weather_dir' - argDict = checkArgs(args, p) - assert argDict['wmLoc'] == 'weather_dir' - - -def test_checkArgs_outloc_3(parsed_args): - '''Tests that the weather model directory gets created when needed''' - args, p = parsed_args - argDict = checkArgs(args, p) - assert os.path.isdir(argDict['wmLoc']) - - -def test_checkArgs_outloc_4(parsed_args): - '''Tests for creating writeable weather model directory''' - args, p = parsed_args - argDict = checkArgs(args, p) - - assert isWriteable(argDict['wmLoc']) - - -def test_ll_bounds_1(parsed_args): - '''Tests that lats out of bounds raises error''' - args, p = parsed_args - args.query_area[0] = -91 - with pytest.raises(ValueError): - checkArgs(args, p) - - -def test_ll_bounds_2(parsed_args): - '''Tests that lats out of bounds raises error''' - args, p = parsed_args - args.query_area[1] = 91 - with pytest.raises(ValueError): - checkArgs(args, p) - - -def test_los_1(parsed_args): - '''Tests that lats out of bounds raises error''' - args, p = parsed_args - args.lineofsight = 'los.rdr' - argDict = checkArgs(args, p) - assert isinstance(argDict['los'], Conventional) - assert argDict['los']._file == 'los.rdr' - - -def test_los_2(parsed_args): - '''Tests that lats out of bounds raises error''' - args, p = parsed_args - args.statevectors = 'sv.txt' - argDict = checkArgs(args, p) - assert isinstance(argDict['los'], Conventional) - assert argDict['los']._file == 'sv.txt' - - -def test_los_3(parsed_args): - '''Tests that lats out of bounds raises error''' - args, p = parsed_args - argDict = checkArgs(args, p) - assert isinstance(argDict['los'], Zenith) - - -def test_models_1a(parsed_args): - '''Tests that the weather model gets passed through correctly''' - args, p = parsed_args - argDict = checkArgs(args, p) - assert argDict['weather_model']['type'].Model() == 'ERA-5' - assert argDict['weather_model']['name'] == 'era5' - - -def test_models_1b(parsed_args): - '''Tests that the weather model gets passed through correctly''' - args, p = parsed_args - args.model = 'HRRR' - argDict = checkArgs(args, p) - assert argDict['weather_model']['type'].Model() == 'HRRR' - assert argDict['weather_model']['name'] == 'hrrr' - - -def test_models_1c(parsed_args): - '''Tests that the weather model gets passed through correctly''' - args, p = parsed_args - args.model = 'NCMR' - argDict = checkArgs(args, p) - assert argDict['weather_model']['type'].Model() == 'NCMR' - assert argDict['weather_model']['name'] == 'ncmr' - - -def test_models_1d(parsed_args): - '''Tests that the weather model gets passed through correctly''' - args, p = parsed_args - args.model = 'era-5' - argDict = checkArgs(args, p) - assert argDict['weather_model']['type'].Model() == 'ERA-5' - assert argDict['weather_model']['name'] == 'era5' - - -def test_models_1e(parsed_args): - '''Tests that the weather model gets passed through correctly''' - args, p = parsed_args - args.model = 'ERA-5' - argDict = checkArgs(args, p) - assert argDict['weather_model']['type'].Model() == 'ERA-5' - assert argDict['weather_model']['name'] == 'era5' - - -def test_models_1f(parsed_args): - '''Tests that the weather model gets passed through correctly''' - args, p = parsed_args - args.model = 'Era-5' - argDict = checkArgs(args, p) - assert argDict['weather_model']['type'].Model() == 'ERA-5' - assert argDict['weather_model']['name'] == 'era5' - - -def test_models_2(parsed_args): - '''Tests that unknown weather models get rejected''' - args, p = parsed_args - args.model = 'unknown' - with pytest.raises(NotImplementedError): - checkArgs(args, p) - - -def test_models_3a(parsed_args): - '''Tests that WRF weather models requires files''' - args, p = parsed_args - args.model = 'WRF' - with pytest.raises(RuntimeError): - checkArgs(args, p) - - -def test_models_3b(parsed_args): - '''Tests that HDF5 weather models requires files''' - args, p = parsed_args - args.model = 'HDF5' - with pytest.raises(RuntimeError): - checkArgs(args, p) - - -def test_models_3c(parsed_args): - '''Tests that WRF weather models requires files''' - args, p = parsed_args - args.model = 'WRF' - args.files = ['file1.wrf', 'file2.wrf'] - # argDict = checkArgs(args, p) - # TODO - assert True - - -def test_zref_1(parsed_args): - '''tests that default zref gets generated''' - args, p = parsed_args - argDict = checkArgs(args, p) - assert argDict['zref'] == _ZREF - - -def test_zref_2(parsed_args): - '''tests that default zref gets generated''' - ztest = 20000 - args, p = parsed_args - args.zref = ztest - argDict = checkArgs(args, p) - assert argDict['zref'] == ztest - - -def test_parallel_1(parsed_args): - '''tests that parallel options are handled correctly''' - args, p = parsed_args - argDict = checkArgs(args, p) - assert argDict['parallel'] == 1 - - -def test_parallel_2(parsed_args): - '''tests that parallel options are handled correctly''' - args, p = parsed_args - args.parallel = 'all' - argDict = checkArgs(args, p) - assert argDict['parallel'] == mp.cpu_count() - - -def test_parallel_3(parsed_args): - '''tests that parallel options are handled correctly''' - args, p = parsed_args - args.parallel = 2 - argDict = checkArgs(args, p) - assert argDict['parallel'] == 2 - - -def test_parallel_4(parsed_args): - '''tests that parallel options are handled correctly''' - args, p = parsed_args - args.parallel = 2000 - argDict = checkArgs(args, p) - assert argDict['parallel'] == mp.cpu_count() - - -def test_verbose_1(parsed_args): - '''tests that verbose option is handled correctly''' - args, p = parsed_args - argDict = checkArgs(args, p) - assert not argDict['verbose'] - - -def test_verbose_2(parsed_args): - '''tests that verbose option is handled correctly''' - args, p = parsed_args - args.verbose = True - argDict = checkArgs(args, p) - assert argDict['verbose'] - - -def test_download_only_1(parsed_args): - '''tests that the download-only option is handled correctly''' - args, p = parsed_args - argDict = checkArgs(args, p) - assert not argDict['download_only'] - - -def test_download_only_2(parsed_args): - '''tests that the download-only option is handled correctly''' - args, p = parsed_args - args.download_only = True - argDict = checkArgs(args, p) - assert argDict['download_only'] - - -def test_useWeatherNodes_1(parsed_args): - '''tests that the correct flag gets passed''' - args, p = parsed_args - argDict = checkArgs(args, p) - assert argDict['flag'] == 'bounding_box' # default arguments use a bounding box - - -def test_filenames_1(parsed_args): - '''tests that the correct filenames are generated''' - args, p = parsed_args - argDict = checkArgs(args, p) - assert 'Delay' not in argDict['wetFilenames'][0] - assert 'wet' in argDict['wetFilenames'][0] - assert 'hydro' in argDict['hydroFilenames'][0] - assert '20200103' in argDict['wetFilenames'][0] - assert '20200103' in argDict['hydroFilenames'][0] - assert len(argDict['hydroFilenames']) == 1 - - -def test_filenames_2(parsed_args): - '''tests that the correct filenames are generated''' - args, p = parsed_args - args.query_area = os.path.join(SCENARIO_2, 'stations.csv') - argDict = checkArgs(args, p) - assert 'Delay' in argDict['wetFilenames'][0] - assert '20200103' in argDict['wetFilenames'][0] - assert len(argDict['wetFilenames']) == 1 - - -def test_makeDelayFileNames_1(): - assert makeDelayFileNames(None, None, "h5", "name", "dir") == \ - ("dir/name_wet_ztd.h5", "dir/name_hydro_ztd.h5") - - -def test_makeDelayFileNames_2(): - assert makeDelayFileNames(None, (), "h5", "name", "dir") == \ - ("dir/name_wet_std.h5", "dir/name_hydro_std.h5") - - -def test_makeDelayFileNames_3(): - assert makeDelayFileNames(datetime.datetime(2020, 1, 1, 1, 2, 3), None, "h5", "model_name", "dir") == \ - ( - "dir/model_name_wet_20200101T010203_ztd.h5", - "dir/model_name_hydro_20200101T010203_ztd.h5" - ) - - -def test_makeDelayFileNames_4(): - assert makeDelayFileNames(datetime.datetime(1900, 12, 31, 1, 2, 3), "los", "h5", "model_name", "dir") == \ - ( - "dir/model_name_wet_19001231T010203_std.h5", - "dir/model_name_hydro_19001231T010203_std.h5" - ) - - -def test_model2module(): - model_module_name, model_obj = modelName2Module('ERA5') - assert model_obj().Model() == 'ERA-5' - - -def test_dem_1(parsed_args): - '''Test that passing a raster format with height levels throws an error''' - args, p = parsed_args - argDict = checkArgs(args, p) - assert argDict['heights'][0] == 'skip' - assert argDict['heights'][1] is None - - -def test_dem_2(parsed_args): - '''Test that passing a raster format with height levels throws an error''' - args, p = parsed_args - args.heightlvs = [10, 100, 1000] - argDict = checkArgs(args, p) - assert argDict['heights'][0] == 'lvs' - assert np.allclose(argDict['heights'][1], [10, 100, 1000]) - - -def test_dem_3(parsed_args): - '''Test that passing a raster format with height levels throws an error''' - args, p = parsed_args - args.heightlvs = [10, 100, 1000] - args.query_area = os.path.join(SCENARIO_2, 'stations.csv') - argDict = checkArgs(args, p) - assert argDict['heights'][0] == 'lvs' - assert np.allclose(argDict['heights'][1], [10, 100, 1000]) - - -def test_dem_4(parsed_args): - '''Test that passing a raster format with height levels throws an error''' - args, p = parsed_args - args.query_area = os.path.join(SCENARIO_2, 'stations.csv') - argDict = checkArgs(args, p) - assert argDict['heights'][0] == 'pandas' - assert argDict['heights'][1][0] == argDict['wetFilenames'][0] - - -def test_dem_5(parsed_args): - '''Test that passing a raster format with height levels throws an error''' - args, p = parsed_args - args.query_area = [os.path.join(SCENARIO_1, 'geom', 'lat.dat'), os.path.join(SCENARIO_1, 'geom', 'lat.dat')] - argDict = checkArgs(args, p) - assert argDict['heights'][0] == 'download' - assert argDict['heights'][1] == os.path.join(argDict['out'], 'geom', 'warpedDEM.dem') diff --git a/test/scenario_3/HRRR_tropo_20181113T230000_ray.nc b/test/scenario_3/HRRR_tropo_20181113T230000_ray.nc new file mode 100644 index 000000000..8ba4eed64 Binary files /dev/null and b/test/scenario_3/HRRR_tropo_20181113T230000_ray.nc differ diff --git a/test/scenario_3/HRRR_tropo_20181113T230000_std.nc b/test/scenario_3/HRRR_tropo_20181113T230000_std.nc index 8ba4eed64..063de4f1a 100644 Binary files a/test/scenario_3/HRRR_tropo_20181113T230000_std.nc and b/test/scenario_3/HRRR_tropo_20181113T230000_std.nc differ diff --git a/test/scenario_3/raider_example_3.yaml b/test/scenario_3/raider_example_3.yaml index b336ab452..94d298c8c 100644 --- a/test/scenario_3/raider_example_3.yaml +++ b/test/scenario_3/raider_example_3.yaml @@ -11,7 +11,7 @@ aoi_group: height_group: height_levels: 0 100 500 1000 los_group: - ray_trace: True # Use projected slant delay by default + ray_trace: True orbit_file: test/scenario_3/S1A_OPER_AUX_POEORB_OPOD_20181203T120749_V20181112T225942_20181114T005942.EOF runtime_group: output_projection: 4326 diff --git a/test/scenario_3/raider_example_3_proj.yaml b/test/scenario_3/raider_example_3_proj.yaml new file mode 100644 index 000000000..5f26ffb2b --- /dev/null +++ b/test/scenario_3/raider_example_3_proj.yaml @@ -0,0 +1,18 @@ +# vim: set filetype=yaml: +look_dir: right +date_group: + date_list: [20181113] +time_group: + time: "23:00:00" + end_time: +weather_model: HRRR +aoi_group: + bounding_box: 36.8 36.85 -76.15 -76.05 +height_group: + height_levels: 0 100 500 1000 +los_group: + ray_trace: False # Use projected slant delay by default + orbit_file: test/scenario_3/S1A_OPER_AUX_POEORB_OPOD_20181203T120749_V20181112T225942_20181114T005942.EOF +runtime_group: + output_projection: 4326 + cube_spacing_in_m: 5000.0 diff --git a/test/scenario_3/raider_example_4.yaml b/test/scenario_3/raider_example_4.yaml deleted file mode 100644 index ce6dc46bf..000000000 --- a/test/scenario_3/raider_example_4.yaml +++ /dev/null @@ -1,17 +0,0 @@ -# vim: set filetype=yaml: - look_dir: right - date_group: - date_list: [20181113] - time_group: - time: "23:00:00" - end_time: - weather_model: HRRR - aoi_group: - bounding_box: 36.8 36.85 -76.15 -76.05 - los_group: - ray_trace: False # Use projected slant delay by default - orbit_file: test/scenario_3/S1A_OPER_AUX_POEORB_OPOD_20181203T120749_V20181112T225942_20181114T005942.EOF - runtime_group: - output_projection: 4326 - cube_spacing_in_m: 250.0 - \ No newline at end of file diff --git a/test/test_checkArgs.py b/test/test_checkArgs.py new file mode 100644 index 000000000..e0699e733 --- /dev/null +++ b/test/test_checkArgs.py @@ -0,0 +1,180 @@ +import datetime +import os +import pytest + +import multiprocessing as mp +import numpy as np +import pandas as pd + +from test import TEST_DIR, pushd + +from RAiDER.cli import DEFAULT_DICT +from RAiDER.checkArgs import checkArgs, makeDelayFileNames +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") + +@pytest.fixture +def args(): + d = DEFAULT_DICT + d['date_list'] = [datetime.datetime(2018, 1, 1)] + d['time'] = datetime.time(12,0,0) + d['aoi'] = BoundingBox([38, 39, -92, -91]) + d['los'] = Zenith() + d['weather_model'] = GMAO() + + return d + + +def isWriteable(dirpath): + '''Test whether a directory is writeable''' + try: + filehandle = open(os.path.join(dirpath, 'tmp.txt'), 'w') + filehandle.close() + return True + except IOError: + return False + +def test_checkArgs_outfmt_1(args): + '''Test that passing height levels with hdf5 outformat works''' + args = args + args.file_format = 'h5' + args.heightlvls = [10, 100, 1000] + checkArgs(args) + assert os.path.splitext(args.wetFilenames[0])[-1] == '.h5' + + +def test_checkArgs_outfmt_2(args): + '''Test that passing a raster format with height levels throws an error''' + args = args + args.heightlvs = [10, 100, 1000] + args.file_format = 'GTiff' + args = checkArgs(args) + assert os.path.splitext(args.wetFilenames[0])[-1] == '.nc' + + +def test_checkArgs_outfmt_3(args): + '''Test that passing a raster format with height levels throws an error''' + args = args + with pytest.raises(FileNotFoundError): + args.aoi = StationFile(os.path.join('fake_dir', 'stations.csv')) + + +def test_checkArgs_outfmt_4(args): + '''Test that passing a raster format with height levels throws an error''' + args = args + args.aoi = RasterRDR( + lat_file = os.path.join(SCENARIO_1, 'geom', 'lat.dat'), + lon_file = os.path.join(SCENARIO_1, 'geom', 'lon.dat'), + ) + argDict = checkArgs(args) + assert argDict.aoi.type()=='radar_rasters' + + +def test_checkArgs_outfmt_5(args): + '''Test that passing a raster format with height levels throws an error''' + args = args + args.aoi = StationFile(os.path.join(SCENARIO_2, 'stations.csv')) + argDict = checkArgs(args) + assert pd.read_csv(argDict['wetFilenames'][0]).shape == (8, 4) + + +def test_checkArgs_outloc_1(args): + '''Test that the default output and weather model directories are correct''' + args = args + argDict = checkArgs(args) + out = argDict['output_directory'] + wmLoc = argDict['weather_model_directory'] + assert os.path.abspath(out) == os.getcwd() + assert os.path.abspath(wmLoc) == os.path.join(os.getcwd(), 'weather_files') + + +def test_checkArgs_outloc_2(args, tmp_path): + '''Tests that the correct output location gets assigned when provided''' + with pushd(tmp_path): + args = args + args.output_directory = tmp_path + argDict = checkArgs(args) + out = argDict['output_directory'] + assert out == tmp_path + + +def test_checkArgs_outloc_2b(args, tmp_path): + ''' Tests that the weather model directory gets passed through by itself''' + with pushd(tmp_path): + args = args + args.output_directory = tmp_path + args.weather_model_directory = 'weather_dir' + argDict = checkArgs(args) + assert argDict['weather_model_directory'] == 'weather_dir' + + +def test_checkArgs_outloc_3(args, tmp_path): + '''Tests that the weather model directory gets created when needed''' + with pushd(tmp_path): + args = args + args.output_directory = tmp_path + argDict = checkArgs(args) + assert os.path.isdir(argDict['weather_model_directory']) + + +def test_checkArgs_outloc_4(args): + '''Tests for creating writeable weather model directory''' + args = args + argDict = checkArgs(args) + + assert isWriteable(argDict['weather_model_directory']) + + +def test_filenames_1(args): + '''tests that the correct filenames are generated''' + args = args + argDict = checkArgs(args) + assert 'Delay' not in argDict['wetFilenames'][0] + assert 'wet' in argDict['wetFilenames'][0] + assert 'hydro' in argDict['hydroFilenames'][0] + assert '20180101' in argDict['wetFilenames'][0] + assert '20180101' in argDict['hydroFilenames'][0] + assert len(argDict['hydroFilenames']) == 1 + + +def test_filenames_2(args): + '''tests that the correct filenames are generated''' + args = args + args.aoi = StationFile(os.path.join(SCENARIO_2, 'stations.csv')) + argDict = checkArgs(args) + assert '20180101' in argDict['wetFilenames'][0] + assert len(argDict['wetFilenames']) == 1 + + +def test_makeDelayFileNames_1(): + assert makeDelayFileNames(None, None, "h5", "name", "dir") == \ + ("dir/name_wet_ztd.h5", "dir/name_hydro_ztd.h5") + + +def test_makeDelayFileNames_2(): + assert makeDelayFileNames(None, (), "h5", "name", "dir") == \ + ("dir/name_wet_std.h5", "dir/name_hydro_std.h5") + + +def test_makeDelayFileNames_3(): + assert makeDelayFileNames(datetime.datetime(2020, 1, 1, 1, 2, 3), None, "h5", "model_name", "dir") == \ + ( + "dir/model_name_wet_20200101T010203_ztd.h5", + "dir/model_name_hydro_20200101T010203_ztd.h5" + ) + + +def test_makeDelayFileNames_4(): + assert makeDelayFileNames(datetime.datetime(1900, 12, 31, 1, 2, 3), "los", "h5", "model_name", "dir") == \ + ( + "dir/model_name_wet_19001231T010203_std.h5", + "dir/model_name_hydro_19001231T010203_std.h5" + ) + + + diff --git a/test/test_llreader.py b/test/test_llreader.py index dff8681a6..71da892dd 100644 --- a/test/test_llreader.py +++ b/test/test_llreader.py @@ -6,7 +6,7 @@ from test import GEOM_DIR, TEST_DIR -from RAiDER.cli.raider import create_parser +from RAiDER.cli.raider import calcDelays from RAiDER.utilFcns import rio_open from RAiDER.llreader import ( @@ -20,7 +20,7 @@ @pytest.fixture def parser(): - return create_parser() + return calcDelays() @pytest.fixture @@ -33,6 +33,14 @@ def llfiles(): return os.path.join(SCENARIO1_DIR, 'lat.dat'), os.path.join(SCENARIO1_DIR, 'lon.dat') +def test_latlon_reader_2(): + with pytest.raises(ValueError): + RasterRDR(lat_file=None, lon_file=None) + + with pytest.raises(ValueError): + RasterRDR(lat_file='doesnotexist.rdr', lon_file='doesnotexist.rdr') + + def test_latlon_reader(): latfile = os.path.join(GEOM_DIR, 'lat.rdr') lonfile = os.path.join(GEOM_DIR, 'lon.rdr') @@ -88,3 +96,9 @@ def test_bounds_from_csv(station_file): bounds_true = [33.746, 36.795, -118.312, -114.892] snwe = bounds_from_csv(station_file) assert all([np.allclose(b, t) for b, t in zip(snwe, bounds_true)]) + + +def test_readZ_sf(station_file): + aoi = StationFile(station_file) + assert np.allclose(aoi.readZ(), .1) + diff --git a/test/test_scenario_3.py b/test/test_scenario_3.py new file mode 100644 index 000000000..b168508dd --- /dev/null +++ b/test/test_scenario_3.py @@ -0,0 +1,27 @@ +import os +import pytest +import subprocess + +from test import TEST_DIR + +import numpy as np +import xarray as xr + + + +def test_scenario_3(): + SCENARIO_DIR = os.path.join(TEST_DIR, "scenario_3") + + test_path = os.path.join(SCENARIO_DIR, 'raider_example_3.yaml') + process = subprocess.run(['raider.py', test_path],stdout=subprocess.PIPE, universal_newlines=True) + assert process.returncode == 0 + + new_data = xr.load_dataset('HRRR_tropo_20181113T230000_ray.nc') + golden_data = xr.load_dataset(os.path.join(SCENARIO_DIR, 'HRRR_tropo_20181113T230000_ray.nc')) + + assert np.allclose(golden_data['wet'], new_data['wet']) + assert np.allclose(golden_data['hydro'], new_data['hydro']) + + # Clean up files + subprocess.run(['rm', '-f', './HRRR*']) + subprocess.run(['rm', '-rf', './weather_files']) diff --git a/test/test_scenario_3_proj.py b/test/test_scenario_3_proj.py new file mode 100644 index 000000000..5e3c07af1 --- /dev/null +++ b/test/test_scenario_3_proj.py @@ -0,0 +1,27 @@ +import os +import pytest +import subprocess + +from test import TEST_DIR + +import numpy as np +import xarray as xr + + +# @pytest.mark.skip +def test_scenario_3_proj(): + SCENARIO_DIR = os.path.join(TEST_DIR, "scenario_3") + + test_path = os.path.join(SCENARIO_DIR, 'raider_example_3_proj.yaml') + process = subprocess.run(['raider.py', test_path],stdout=subprocess.PIPE, universal_newlines=True) + assert process.returncode == 0 + + new_data = xr.load_dataset('HRRR_tropo_20181113T230000_std.nc') + golden_data = xr.load_dataset(os.path.join(SCENARIO_DIR, 'HRRR_tropo_20181113T230000_std.nc')) + + assert np.allclose(golden_data['wet'], new_data['wet']) + assert np.allclose(golden_data['hydro'], new_data['hydro']) + + # Clean up files + subprocess.run(['rm', '-f', './HRRR*']) + subprocess.run(['rm', '-rf', './weather_files']) \ No newline at end of file diff --git a/test/test_scenarios.py b/test/test_scenarios.py index 42adde142..26ce04511 100644 --- a/test/test_scenarios.py +++ b/test/test_scenarios.py @@ -17,6 +17,7 @@ def test_scenario_1(): new_data = xr.load_dataset('HRRR_tropo_20200101T120000_ztd.nc') golden_data = xr.load_dataset(os.path.join(SCENARIO_DIR, 'HRRR_tropo_20200101T120000_ztd.nc')) + assert np.allclose(golden_data['wet'], new_data['wet']) assert np.allclose(golden_data['hydro'], new_data['hydro']) @@ -25,32 +26,3 @@ def test_scenario_1(): subprocess.run(['rm', '-f', './HRRR*']) subprocess.run(['rm', '-rf', './weather_files']) - -#TODO: Next release include GNSS station test -# @pytest.mark.long -# def test_scenario_2(): -# test_path = os.path.join(TEST_DIR, "scenario_2", 'raider_example_2.yaml') -# process = subprocess.run(['raider.py', test_path],stdout=subprocess.PIPE, universal_newlines=True) -# assert process.returncode == 0 - -# # Clean up files -# subprocess.run(['rm', '-f', './HRRR*']) -# subprocess.run(['rm', '-rf', './weather_files']) - - -def test_scenario_3(): - SCENARIO_DIR = os.path.join(TEST_DIR, "scenario_3") - - test_path = os.path.join(SCENARIO_DIR, 'raider_example_3.yaml') - process = subprocess.run(['raider.py', test_path],stdout=subprocess.PIPE, universal_newlines=True) - assert process.returncode == 0 - - new_data = xr.load_dataset('HRRR_tropo_20181113T230000_std.nc') - golden_data = xr.load_dataset(os.path.join(SCENARIO_DIR, 'HRRR_tropo_20181113T230000_std.nc')) - - assert np.allclose(golden_data['wet'], new_data['wet']) - assert np.allclose(golden_data['hydro'], new_data['hydro']) - - # Clean up files - subprocess.run(['rm', '-f', './HRRR*']) - subprocess.run(['rm', '-rf', './weather_files']) diff --git a/test/weather_model/test_downloaders.py b/test/weather_model/test_downloaders.py index 6f86dc118..bd1237197 100644 --- a/test/weather_model/test_downloaders.py +++ b/test/weather_model/test_downloaders.py @@ -10,7 +10,8 @@ from RAiDER.models.era5t import ERA5T from RAiDER.models.erai import ERAI -@pytest.mark.skip + +@pytest.mark.long def test_era5(): wm = ERA5() wm.fetch( @@ -19,7 +20,8 @@ def test_era5(): datetime(2020, 1, 1, 0, 0, 0) ) -@pytest.mark.skip + +@pytest.mark.long def test_era5t(): wm = ERA5T() wm.fetch( @@ -29,6 +31,7 @@ def test_era5t(): ) +@pytest.mark.long def test_erai(): wm = ERAI() wm.fetch( diff --git a/test/weather_model/test_weather_model.py b/test/weather_model/test_weather_model.py index 0cf5cc633..74ddd4f28 100644 --- a/test/weather_model/test_weather_model.py +++ b/test/weather_model/test_weather_model.py @@ -13,7 +13,7 @@ from pyproj import CRS from RAiDER.constants import _ZMIN, _ZREF -from RAiDER.delay import build_cube_ray +from RAiDER.delay import _build_cube_ray from RAiDER.losreader import state_to_los from RAiDER.models.weatherModel import ( WeatherModel, @@ -240,7 +240,7 @@ def test_build_cube_ray(setup_fake_raytracing, model): #TODO: Check that the look vectors are not nans lv, xyz = state_to_los(svs, np.stack([_Y.ravel(), _X.ravel(), _Z.ravel()], axis=-1),out="ecef") - out = build_cube_ray(xs, ys, zs, orb, look_dir, CRS(4326), CRS(4326), [m.interpWet(), m.interpHydro()], elp=elp) + out = _build_cube_ray(xs, ys, zs, orb, look_dir, CRS(4326), CRS(4326), [m.interpWet(), m.interpHydro()], elp=elp) assert out.shape == out_true.shape diff --git a/tools/RAiDER/aria/calcGUNW.py b/tools/RAiDER/aria/calcGUNW.py new file mode 100644 index 000000000..79ea59ca1 --- /dev/null +++ b/tools/RAiDER/aria/calcGUNW.py @@ -0,0 +1,8 @@ +def main(args): + """ Calculate interferometric phase delay at dates in ARIA standard product (GUNW) """ + pass + + + +if __name__ == '__main__': + main() diff --git a/tools/RAiDER/aria/prepFromGUNW.py b/tools/RAiDER/aria/prepFromGUNW.py new file mode 100644 index 000000000..6076d5b02 --- /dev/null +++ b/tools/RAiDER/aria/prepFromGUNW.py @@ -0,0 +1,11 @@ +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Author: Jeremy Maurer, Brett Buzzanga +# Copyright 2022, by the California Institute of Technology. ALL RIGHTS +# RESERVED. United States Government Sponsorship acknowledged. +# +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +def main(args): + """ Get needed RAiDER parameters from an ARIA standard product (GUNW) """ + pass diff --git a/tools/RAiDER/checkArgs.py b/tools/RAiDER/checkArgs.py index 774b548f5..04443f7c8 100644 --- a/tools/RAiDER/checkArgs.py +++ b/tools/RAiDER/checkArgs.py @@ -8,11 +8,14 @@ import os import pandas as pd +import rasterio.drivers as rd from datetime import datetime + from RAiDER.losreader import Zenith from RAiDER.llreader import BoundingBox +from RAiDER.logger import logger def checkArgs(args): @@ -34,36 +37,54 @@ def checkArgs(args): # filenames wetNames, hydroNames = [], [] for d in args.date_list: - if not args.aoi is not BoundingBox: - if args.station_file is not None: + if (args.aoi.type() != 'bounding_box'): + + # Handle the GNSS station file + if (args.aoi.type()=='station_file'): wetFilename = os.path.join( args.output_directory, '{}_Delay_{}.csv' .format( - args.weather_model, - args.time.strftime('%Y%m%dT%H%M%S'), + args.weather_model.Model(), + d.strftime('%Y%m%dT%H%M%S'), ) ) - hydroFilename = wetFilename + hydroFilename = None # only the 'wetFilename' is used for the station_file - # copy the input file to the output location for editing - indf = pd.read_csv(args.query_area).drop_duplicates(subset=["Lat", "Lon"]) + # copy the input station file to the output location for editing + indf = pd.read_csv(args.aoi._filename).drop_duplicates(subset=["Lat", "Lon"]) indf.to_csv(wetFilename, index=False) else: - wetNames.append(None) - hydroNames.append(None) + # This implies rasters + fmt = get_raster_ext(args.file_format) + wetFilename, hydroFilename = makeDelayFileNames( + d, + args.los, + fmt, + args.weather_model._dataset.upper(), + args.output_directory, + ) + + else: + # In this case a cube file format is needed + if args.file_format not in '.nc .h5 h5 hdf5 .hdf5 nc'.split(): + fmt = 'nc' + logger.debug('Invalid extension %s for cube. Defaulting to .nc', args.file_format) + else: + fmt = args.file_format.strip('.').replace('df', '') + wetFilename, hydroFilename = makeDelayFileNames( d, args.los, - args.raster_format, + fmt, args.weather_model._dataset.upper(), args.output_directory, ) - - wetNames.append(wetFilename) - hydroNames.append(hydroFilename) + + wetNames.append(wetFilename) + hydroNames.append(hydroFilename) args.wetFilenames = wetNames args.hydroFilenames = hydroNames @@ -71,6 +92,20 @@ def checkArgs(args): return args +def get_raster_ext(fmt): + drivers = rd.raster_driver_extensions() + extensions = {value.upper():key for key, value in drivers.items()} + + # add in ENVI/ISCE formats with generic extension + extensions['ENVI'] = '.dat' + extensions['ISCE'] = '.dat' + + try: + return extensions[fmt.upper()] + except KeyError: + raise ValueError('{} is not a valid gdal/rasterio file format for rasters'.format(fmt)) + + def makeDelayFileNames(time, los, outformat, weather_model_name, out): ''' return names for the wet and hydrostatic delays. diff --git a/tools/RAiDER/cli/__init__.py b/tools/RAiDER/cli/__init__.py index e69de29bb..79a9cc19e 100644 --- a/tools/RAiDER/cli/__init__.py +++ b/tools/RAiDER/cli/__init__.py @@ -0,0 +1,45 @@ +import os +from RAiDER.constants import _ZREF, _CUBE_SPACING_IN_M + +class AttributeDict(dict): + __getattr__ = dict.__getitem__ + __setattr__ = dict.__setitem__ + __delattr__ = dict.__delitem__ + +DEFAULT_DICT = AttributeDict( + dict( + look_dir='right', + date_start=None, + date_end=None, + date_step=None, + date_list=None, + time=None, + end_time=None, + weather_model=None, + lat_file=None, + lon_file=None, + station_file=None, + bounding_box=None, + geocoded_file=None, + dem=None, + use_dem_latlon=False, + height_levels=None, + height_file_rdr=None, + ray_trace=False, + zref=_ZREF, + cube_spacing_in_m=_CUBE_SPACING_IN_M, + los_file=None, + los_convention='isce', + los_cube=None, + orbit_file=None, + verbose=True, + raster_format='GTiff', + file_format='GTiff', + output_directory=os.getcwd(), + weather_model_directory=os.path.join( + os.getcwd(), + 'weather_files' + ), + output_projection='EPSG:4236', + ) + ) diff --git a/tools/RAiDER/cli/__main__.py b/tools/RAiDER/cli/__main__.py new file mode 100644 index 000000000..b1cc91a95 --- /dev/null +++ b/tools/RAiDER/cli/__main__.py @@ -0,0 +1,25 @@ +import argparse +import os +from importlib.metadata import entry_points + +from RAiDER.cli.raider import calcDelays, downloadGNSS, calcDelaysGUNW + +def main(): + parser = argparse.ArgumentParser( + prefix_chars='+', + formatter_class=argparse.ArgumentDefaultsHelpFormatter + ) + parser.add_argument( + '++process', choices=['calcDelays', 'downloadGNSS', 'calcDelaysGUNW'], + default='calcDelays', + help='Select the entrypoint to use' + ) + args, unknowns = parser.parse_known_args() + os.sys.argv = [args.process, *unknowns] + + process_entry_point = entry_points(group='console_scripts', name=f'{args.process}.py')[0] + process_entry_point.load()() + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/tools/RAiDER/cli/raider.py b/tools/RAiDER/cli/raider.py index b417732b7..8613a4f82 100644 --- a/tools/RAiDER/cli/raider.py +++ b/tools/RAiDER/cli/raider.py @@ -8,11 +8,10 @@ import RAiDER from RAiDER.constants import _ZREF, _CUBE_SPACING_IN_M from RAiDER.logger import logger, logging -from RAiDER.cli.validators import (enforce_time, enforce_bbox, parse_dates, - get_query_region, get_heights, get_los, enforce_wm) - -from RAiDER.checkArgs import checkArgs -from RAiDER.delay import main as main_delay +from RAiDER.cli import DEFAULT_DICT, AttributeDict +from RAiDER.cli.validators import ( + enforce_time, enforce_bbox, parse_dates, get_query_region, get_heights, get_los, enforce_wm +) HELP_MESSAGE = """ @@ -32,47 +31,6 @@ raider.py customTemplatefile.cfg """ -class AttributeDict(dict): - __getattr__ = dict.__getitem__ - __setattr__ = dict.__setitem__ - __delattr__ = dict.__delitem__ - -DEFAULT_DICT = dict( - look_dir='right', - date_start=None, - date_end=None, - date_step=None, - date_list=None, - time=None, - end_time=None, - weather_model=None, - lat_file=None, - lon_file=None, - station_file=None, - bounding_box=None, - geocoded_file=None, - dem=None, - use_dem_latlon=False, - height_levels=None, - height_file_rdr=None, - ray_trace=False, - zref=_ZREF, - cube_spacing_in_m=_CUBE_SPACING_IN_M, # TODO - Where are these parsed? - los_file=None, - los_convention='isce', - los_cube=None, - orbit_file=None, - verbose=True, - raster_format='GTiff', - output_directory=os.getcwd(), - weather_model_directory=os.path.join( - os.getcwd(), - 'weather_files' - ), - output_projection='EPSG:4236', - ) - - def create_parser(): """Parse command line arguments using argparse.""" p = argparse.ArgumentParser( @@ -156,14 +114,17 @@ def parseCMD(iargs=None): def read_template_file(fname): """ Read the template file into a dictionary structure. - Parameters: fname - str, full path to the template file - delimiter - str, string to separate the key and value - skip_chars - list of str, skip certain charaters in values - Returns: template - dict, file content - Examples: template = read_template('raider.yaml') + Args: + fname (str): full path to the template file + Returns: + dict: arguments to pass to RAiDER functions + + Examples: + >>> template = read_template_file('raider.yaml') - Modified from MintPy's 'read_template' """ + from RAiDER.cli.validators import (enforce_time, enforce_bbox, parse_dates, + get_query_region, get_heights, get_los, enforce_wm) with open(fname, 'r') as f: try: params = yaml.safe_load(f) @@ -231,30 +192,303 @@ def drop_nans(d): return d -########################################################################## -def main(iargs=None): - # parse - inps = parseCMD(iargs) +def calcDelays(iargs=None): + """ Parse command line arguments using argparse. """ + import RAiDER + from RAiDER.delay import tropo_delay + from RAiDER.checkArgs import checkArgs + from RAiDER.processWM import prepareWeatherModel + from RAiDER.utilFcns import writeDelays + examples = 'Examples of use:' \ + '\n\t raider.py customTemplatefile.cfg' \ + '\n\t raider.py -g' + + p = argparse.ArgumentParser( + description = + 'Command line interface for RAiDER processing with a configure file.' + 'Default options can be found by running: raider.py --generate_config', + epilog=examples, formatter_class=argparse.RawDescriptionHelpFormatter) + + p.add_argument( + 'customTemplateFile', nargs='?', + help='custom template with option settings.\n' + + "ignored if the default smallbaselineApp.cfg is input." + ) + + p.add_argument( + '-g', '--generate_template', action='store_true', + help='generate default template (if it does not exist) and exit.' + ) + + p.add_argument( + '--download_only', action='store_true', + help='only download a weather model.' + ) + + ## if not None, will replace first argument (customTemplateFile) + args = p.parse_args(args=iargs) + + # default input file + template_file = os.path.join(os.path.dirname(RAiDER.__file__), + 'cli', 'raider.yaml') + + if args.generate_template: + dst = os.path.join(os.getcwd(), 'raider.yaml') + shutil.copyfile(template_file, dst) + logger.info('Wrote %s', dst) + os.sys.exit() + + + # check: existence of input template files + if (not args.customTemplateFile + and not os.path.isfile(os.path.basename(template_file)) + and not args.generate_template): + msg = "No template file found! It requires that either:" + msg += "\n a custom template file, OR the default template " + msg += "\n file 'raider.yaml' exists in current directory." + + p.print_usage() + print(examples) + raise SystemExit(f'ERROR: {msg}') + + if args.customTemplateFile: + # check the existence + if not os.path.isfile(args.customTemplateFile): + raise FileNotFoundError(args.customTemplateFile) + + args.customTemplateFile = os.path.abspath(args.customTemplateFile) + else: + args.customTemplateFile = template_file # Read the template file - params = read_template_file(inps.customTemplateFile) + params = read_template_file(args.customTemplateFile) # Argument checking - params = checkArgs(params) - - params['download_only'] = inps.download_only + params = checkArgs(params) if not params.verbose: logger.setLevel(logging.INFO) - + delay_dct = {} for t, w, f in zip( params['date_list'], params['wetFilenames'], params['hydroFilenames'] ): + + los = params['los'] + aoi = params['aoi'] + model = params['weather_model'] + + if los.ray_trace(): + ll_bounds = aoi.add_buffer(buffer=1) # add a buffer for raytracing + else: + ll_bounds = aoi.bounds() + + ########################################################### + # weather model calculation + logger.debug('Starting to run the weather model calculation') + logger.debug('Time: {}'.format(t.strftime('%Y%m%d'))) + logger.debug('Beginning weather model pre-processing') try: - (_, _) = main_delay(t, w, f, params) + weather_model_file = prepareWeatherModel( + model, t, + ll_bounds=ll_bounds, # SNWE + wmLoc=params['weather_model_directory'], + makePlots=params['verbose'], + ) except RuntimeError: logger.exception("Date %s failed", t) continue + + # Now process the delays + try: + wet_delay, hydro_delay = tropo_delay( + t, weather_model_file, aoi, los, + height_levels = params['height_levels'], + out_proj = params['output_projection'], + look_dir = params['look_dir'], + cube_spacing_m = params['cube_spacing_in_m'], + ) + except RuntimeError: + logger.exception("Date %s failed", t) + continue + + ########################################################### + # Write the delays to file + # Different options depending on the inputs + + if los.is_Projected(): + out_filename = w.replace("_ztd", "_std") + f = f.replace("_ztd", "_std") + elif los.ray_trace(): + out_filename = w.replace("_std", "_ray") + f = f.replace("_std", "_ray") + else: + out_filename = w + + if hydro_delay is None: + # means that a dataset was returned + ds = wet_delay + ext = os.path.splitext(out_filename)[1] + if ext not in ['.nc', '.h5']: + out_filename = f'{os.path.splitext(out_filename)[0]}.nc' + + out_filename = out_filename.replace("wet", "tropo") + + if out_filename.endswith(".nc"): + ds.to_netcdf(out_filename, mode="w") + elif out_filename.endswith(".h5"): + ds.to_netcdf(out_filename, engine="h5netcdf", invalid_netcdf=True) + + else: + if aoi.type() == 'station_file': + out_filename = f'{os.path.splitext(out_filename)[0]}.csv' + + if aoi.type() in ['station_file', 'radar_rasters', 'geocoded_file']: + writeDelays(aoi, wet_delay, hydro_delay, out_filename, f, outformat=params['raster_format']) + + logger.info('Wrote hydro delays to: %s', f) + + logger.info('Wrote wet delays to: %s', out_filename) + + # delay_dct[t] = wet_delay, hydro_delay + delay_dct[t] = out_filename, f + + return delay_dct + + +## ------------------------------------------------------ downloadGNSSDelays.py +def downloadGNSS(): + """Parse command line arguments using argparse.""" + from RAiDER.gnss.downloadGNSSDelays import main as dlGNSS + p = argparse.ArgumentParser( + formatter_class=argparse.RawDescriptionHelpFormatter, + description=""" \ + Check for and download tropospheric zenith delays for a set of GNSS stations from UNR + + Example call to virtually access and append zenith delay information to a CSV table in specified output + directory, across specified range of time (in YYMMDD YYMMDD) and all available times of day, and confined to specified + geographic bounding box : + downloadGNSSdelay.py --out products -y 20100101 20141231 -b '39 40 -79 -78' + + Example call to virtually access and append zenith delay information to a CSV table in specified output + directory, across specified range of time (in YYMMDD YYMMDD) and specified time of day, and distributed globally : + downloadGNSSdelay.py --out products -y 20100101 20141231 --returntime '00:00:00' + + + Example call to virtually access and append zenith delay information to a CSV table in specified output + directory, across specified range of time in 12 day steps (in YYMMDD YYMMDD days) and specified time of day, and distributed globally : + downloadGNSSdelay.py --out products -y 20100101 20141231 12 --returntime '00:00:00' + + Example call to virtually access and append zenith delay information to a CSV table in specified output + directory, across specified range of time (in YYMMDD YYMMDD) and specified time of day, and distributed globally but restricted + to list of stations specified in input textfile : + downloadGNSSdelay.py --out products -y 20100101 20141231 --returntime '00:00:00' -f station_list.txt + + NOTE, following example call to physically download zenith delay information not recommended as it is not + necessary for most applications. + Example call to physically download and append zenith delay information to a CSV table in specified output + directory, across specified range of time (in YYMMDD YYMMDD) and specified time of day, and confined to specified + geographic bounding box : + downloadGNSSdelay.py --download --out products -y 20100101 20141231 --returntime '00:00:00' -b '39 40 -79 -78' + """) + + # Stations to check/download + area = p.add_argument_group( + 'Stations to check/download. Can be a lat/lon bounding box or file, or will run the whole world if not specified') + area.add_argument( + '--station_file', '-f', default=None, dest='station_file', + help=('Text file containing a list of 4-char station IDs separated by newlines')) + area.add_argument( + '-b', '--bounding_box', dest='bounding_box', type=str, default=None, + help="Provide either valid shapefile or Lat/Lon Bounding SNWE. -- Example : '19 20 -99.5 -98.5'") + area.add_argument( + '--gpsrepo', '-gr', default='UNR', dest='gps_repo', + help=('Specify GPS repository you wish to query. Currently supported archives: UNR.')) + + misc = p.add_argument_group("Run parameters") + add_out(misc) + + misc.add_argument( + '--date', dest='dateList', + help=dedent("""\ + Date to calculate delay. + Can be a single date, a list of two dates (earlier, later) with 1-day interval, or a list of two dates and interval in days (earlier, later, interval). + Example accepted formats: + YYYYMMDD or + YYYYMMDD YYYYMMDD + YYYYMMDD YYYYMMDD N + """), + nargs="+", + action=DateListAction, + type=date_type, + required=True + ) + + misc.add_argument( + '--returntime', dest='returnTime', + help="Return delays closest to this specified time. If not specified, the GPS delays for all times will be returned. Input in 'HH:MM:SS', e.g. '16:00:00'", + default=None) + + misc.add_argument( + '--download', + help='Physically download data. Note this option is not necessary to proceed with statistical analyses, as data can be handled virtually in the program.', + action='store_true', dest='download', default=False) + + add_cpus(misc) + add_verbose(misc) + + args = p.parse_args() + + dlGNSS(args) + return + + +## ------------------------------------------------------------ prepFromGUNW.py +def calcDelaysGUNW(iargs=None): + from RAiDER.aria.prepFromGUNW import main as GUNW_prep + from RAiDER.aria.calcGUNW import main as GUNW_calc + + p = argparse.ArgumentParser( + description='Calculate a cube of interferometic delays for GUNW files') + + p.add_argument( + 'file', type=str, + help='1 ARIA GUNW netcdf file' + ) + + p.add_argument( + '-m', '--model', default='HRRR', type=str, + help='Weather model (Default=HRRR).' + ) + + + p.add_argument( + '-o', '--output_directory', default=os.getcwd(), type=str, + help='Directory to store results (Default=./).' + ) + + p.add_argument( + '-w', '--write', default=True, + help='Optionally write the delays into the given GUNW product (Default=True).' + ) + + + args = p.parse_args(args=iargs) + args.argv = iargs if iargs else os.sys.argv[1:] + # args.files = glob.glob(args.files) # eventually support multiple files + + ## below are placeholders and not yet implemented + ## prep the config needed for delay calcs + # path_cfg, wavelength = GUNW_prep(args) + + ## write the delays to disk using config and return dictionary of: + # date: wet/hydro filename + # dct_delays = calcDelays([path_cfg]) + + ## calculate the interferometric phase and write it out + # GUNW_calc(tropoDelayFile, args.file, wavelength, args.output_directory, args.write) + + return + diff --git a/tools/RAiDER/cli/raider.yaml b/tools/RAiDER/cli/raider.yaml index 1ebe8ee24..8e135917e 100644 --- a/tools/RAiDER/cli/raider.yaml +++ b/tools/RAiDER/cli/raider.yaml @@ -121,7 +121,7 @@ los_group: ## runtime_group: verbose: True - raster_format: GTiff # Can be any rasterio-compatible format + file_format: GTiff # Can be any rasterio-compatible format output_directory: . # uses the runtime directory by default weather_model_directory: # Defaults to /weather_files/ output_projection: # Specify the PROJ-compatible projection of the output delays as an EPSG code diff --git a/tools/RAiDER/cli/statsPlot.py b/tools/RAiDER/cli/statsPlot.py index 9b0aa45ff..ff1b2866f 100755 --- a/tools/RAiDER/cli/statsPlot.py +++ b/tools/RAiDER/cli/statsPlot.py @@ -181,7 +181,7 @@ def convert_SI(val, unit_in, unit_out): # adjust if input isn't datetime, and assume it to be part of workflow # e.g. sigZTD filter, already extracted datetime object try: - return eval('val.apply(pd.to_datetime).dt.{}.astype(np.float).astype("Int32")'.format(unit_out)) + return eval('val.apply(pd.to_datetime).dt.{}.astype(float).astype("Int32")'.format(unit_out)) except BaseException: # TODO: Which error(s)? return val @@ -987,8 +987,7 @@ def create_DF(self): # estimate central longitude lines if '--time_lines' specified if self.time_lines and 'Datetime' in self.df.keys(): - self.df['Date_hr'] = self.df['Datetime'].dt.hour.astype( - np.float).astype("Int32") + self.df['Date_hr'] = self.df['Datetime'].dt.hour.astype(float).astype("Int32") # get list of unique times all_hrs = sorted(set(self.df['Date_hr'])) diff --git a/tools/RAiDER/cli/validators.py b/tools/RAiDER/cli/validators.py index 02793f122..21352d312 100755 --- a/tools/RAiDER/cli/validators.py +++ b/tools/RAiDER/cli/validators.py @@ -15,6 +15,7 @@ from RAiDER.llreader import BoundingBox, Geocube, RasterRDR, StationFile, GeocodedFile, Geocube from RAiDER.losreader import Zenith, Conventional, Raytracing from RAiDER.utilFcns import rio_extents, rio_profile +from RAiDER.logger import logger _BUFFER_SIZE = 0.2 # default buffer size in lat/lon degrees @@ -33,17 +34,18 @@ def enforce_wm(value): def get_los(args): - if ('orbit_file' in args.keys()) and (args['orbit_file'] is not None): - if args.ray_trace: + if args.get('orbit_file'): + if args.get('ray_trace'): los = Raytracing(args.orbit_file) else: los = Conventional(args.orbit_file) - elif ('los_file' in args.keys()) and (args['los_file'] is not None): + elif args.get('los_file'): if args.ray_trace: los = Raytracing(args.los_file, args.los_convention) else: los = Conventional(args.los_file, args.los_convention) - elif ('los_cube' in args.keys()) and (args['los_cube'] is not None): + + elif args.get('los_cube'): raise NotImplementedError('LOS_cube is not yet implemented') # if args.ray_trace: # los = Raytracing(args.los_cube) @@ -114,24 +116,32 @@ def get_query_region(args): ''' # Get bounds from the inputs # make sure this is first - if ('use_dem_latlon' in args.keys()) and args['use_dem_latlon']: + if args.get('use_dem_latlon'): query = GeocodedFile(args.dem, is_dem=True) - elif 'lat_file' in args.keys(): - hgt_file = args.get('hgt_file_rdr', None) # only get it if exists - query = RasterRDR(args.lat_file, args.lon_file, hgt_file) + elif args.get('lat_file'): + hgt_file = args.get('height_file_rdr') # only get it if exists + dem_file = args.get('dem') + query = RasterRDR(args.lat_file, args.lon_file, hgt_file, dem_file) - elif 'station_file' in args.keys(): + elif args.get('station_file'): query = StationFile(args.station_file) - elif 'bounding_box' in args.keys(): + elif args.get('bounding_box'): bbox = enforce_bbox(args.bounding_box) if (np.min(bbox[0]) < -90) | (np.max(bbox[1]) > 90): raise ValueError('Lats are out of N/S bounds; are your lat/lon coordinates switched? Should be SNWE') query = BoundingBox(bbox) - elif 'geocoded_file' in args.keys(): - query = GeocodedFile(args.geocoded_file, is_dem=False) + elif args.get('geocoded_file'): + gfile = os.path.basename(args.geocoded_file).upper() + if (gfile.startswith('SRTM') or gfile.startswith('GLO')): + logger.debug('Using user DEM: %s', gfile) + is_dem = True + else: + is_dem = False + + query = GeocodedFile(args.geocoded_file, is_dem=is_dem) ## untested elif 'los_cube' in args.keys(): @@ -149,7 +159,10 @@ def enforce_bbox(bbox): """ Enforce a valid bounding box """ - bbox = [float(d) for d in bbox.strip().split()] + if isinstance(bbox, str): + bbox = [float(d) for d in bbox.strip().split()] + else: + bbox = [float(d) for d in bbox] # Check the bbox if len(bbox) != 4: diff --git a/tools/RAiDER/constants.py b/tools/RAiDER/constants.py index 6f3b61575..7d7c5cbe7 100644 --- a/tools/RAiDER/constants.py +++ b/tools/RAiDER/constants.py @@ -11,7 +11,7 @@ _ZREF = np.float64(15000) # maximum requierd height _STEP = np.float64(15.0) # integration step size in meters -_CUBE_SPACING_IN_M = np.float(2000) # Horizontal spacing of cube +_CUBE_SPACING_IN_M = float(2000) # Horizontal spacing of cube _g0 = np.float64(9.80665) _RE = np.float64(6371008.7714) diff --git a/tools/RAiDER/delay.py b/tools/RAiDER/delay.py index 10a957163..bd593dcec 100755 --- a/tools/RAiDER/delay.py +++ b/tools/RAiDER/delay.py @@ -5,82 +5,118 @@ # RESERVED. United States Government Sponsorship acknowledged. # # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +"""RAiDER tropospheric delay calculation + +This module provides the main RAiDER functionality for calculating +tropospheric wet and hydrostatic delays from a weather model. Weather +models are accessed as NETCDF files and should have "wet" "hydro" +"wet_total" and "hydro_total" fields specified. +""" import os import datetime -import h5py -import numpy as np +import pyproj import xarray -from netCDF4 import Dataset + from pyproj import CRS, Transformer -from pyproj.exceptions import CRSError +from typing import List import isce3.ext.isce3 as isce -from RAiDER.constants import _STEP +import numpy as np + from RAiDER.delayFcns import ( - getInterpolators, - calculate_start_points, - get_delays, + getInterpolators ) -from RAiDER.dem import getHeights -from RAiDER.logger import logger -from RAiDER.llreader import BoundingBox -from RAiDER.losreader import Zenith, Conventional, Raytracing, get_sv, getTopOfAtmosphere -from RAiDER.processWM import prepareWeatherModel +from RAiDER.logger import logger, logging +from RAiDER.losreader import get_sv, getTopOfAtmosphere from RAiDER.utilFcns import ( - writeDelays, projectDelays, writePnts2HDF5, lla2ecef, transform_bbox, clip_bbox, rio_profile, ) -def tropo_delay_cube(dt, wf, args, model_file=None): +############################################################################### +def tropo_delay( + dt, + weather_model_file: str, + aoi, + los, + height_levels: List[float]=None, + out_proj: int | str=4326, + cube_spacing_m: int=None, + look_dir: str='right', + ): """ - raider cube generation function. - - Same as tropo_delay() above. + Calculate integrated delays on query points. + + Args: + dt: Datetime - Datetime object for determining when to calculate delays + weather_model_File: string - Name of the NETCDF file containing a pre-processed weather model + aoi: AOI object - AOI object + los: LOS object - LOS object + height_levels: list - (optional) list of height levels on which to calculate delays. Only needed for cube generation. + out_proj: int,str - (optional) EPSG code for output projection + look_dir: str - (optional) Satellite look direction. Only needed for slant delay calculation + cube_spacing_m: int - (optional) Horizontal spacing in meters when generating cubes + + Returns: + xarray Dataset *or* ndarrays: - wet and hydrostatic delays at the grid nodes / query points. """ - los = args['los'] - weather_model = args['weather_model'] - wmLoc = args['weather_model_directory'] - zref = args['zref'] - download_only = args['download_only'] - verbose = args['verbose'] - aoi = args['aoi'] - cube_spacing = args["cube_spacing_in_m"] - ll_bounds = aoi.bounds() - - try: - crs = CRS(args['output_projection']) - except CRSError: - raise ValueError('output_projection argument is not a valid CRS specifier') + # get heights + if height_levels is None: + with xarray.load_dataset(weather_model_file) as ds: + height_levels = ds.z.values + + #TODO: expose this as library function + ds = _get_delays_on_cube(dt, weather_model_file, aoi.bounds(), height_levels, + los, out_proj=out_proj, cube_spacing_m=cube_spacing_m, look_dir=look_dir) + + if (aoi.type() == 'bounding_box') or (aoi.type() == 'Geocube'): + return ds, None + + else: + # CRS can be an int, str, or CRS object + try: + out_proj = CRS.from_epsg(out_proj) + except pyproj.exceptions.CRSError: + out_proj = out_proj + pnt_proj = CRS.from_epsg(4326) + lats, lons = aoi.readLL() + hgts = aoi.readZ() + pnts = transformPoints(lats, lons, hgts, pnt_proj, out_proj) + if pnts.ndim == 3: + pnts = pnts.transpose(1,2,0) + elif pnts.ndim == 2: + pnts = pnts.T + ifWet, ifHydro = getInterpolators(ds, 'ztd') # the cube from get_delays_on_cube calls the total delays 'wet' and 'hydro' + wetDelay = ifWet(pnts) + hydroDelay = ifHydro(pnts) + + # return the delays (ZTD or STD) + if los.is_Projected(): + los.setTime(dt) + los.setPoints(lats, lons, hgts) + wetDelay = los(wetDelay) + hydroDelay = los(hydroDelay) + + return wetDelay, hydroDelay + + +def _get_delays_on_cube(dt, weather_model_file, ll_bounds, heights, los, out_proj=4326, cube_spacing_m=None, look_dir='right', nproc=1): + """ + raider cube generation function. + """ # For testing multiprocessing # TODO - move this to configuration - nproc = 1 - - # logging - logger.debug('Starting to run the weather model cube calculation') - logger.debug(f'Time: {dt}') - logger.debug(f'Max integration height is {zref:1.1f} m') - logger.debug(f'Output cube projection is {crs.to_wkt()}') - logger.debug(f'Output cube spacing is {cube_spacing} m') - - # We are using zenith model only for now - logger.debug('Beginning weather model pre-processing') - - # If weather model file is not provided - if model_file is None: - weather_model_file = prepareWeatherModel( - weather_model, - dt, - wmLoc=wmLoc, - ll_bounds=ll_bounds, - zref=zref, - download_only=download_only, - makePlots=verbose - ) + crs = CRS(out_proj) + + # Load CRS from weather model file + wm_proj = rio_profile(f"netcdf:{weather_model_file}:t")["crs"] + if wm_proj is None: + print("WARNING: I can't find a CRS in the weather model file, so I will assume you are using WGS84") + wm_proj = CRS.from_epsg(4326) else: - weather_model_file = model_file + wm_proj = CRS.from_wkt(wm_proj.to_wkt()) # Determine the output grid extent here wesn = ll_bounds[2:] + ll_bounds[:2] @@ -89,66 +125,45 @@ def tropo_delay_cube(dt, wf, args, model_file=None): ) # Clip output grid to multiples of spacing - # If output is desired in degrees + if cube_spacing_m is None: + with xarray.load_dataset(weather_model_file) as ds: + xpts = ds.x.values + ypts = ds.y.values + cube_spacing_m = np.nanmean([np.nanmean(np.diff(xpts)), np.nanmean(np.diff(ypts))]) + if wm_proj.axis_info[0].unit_name == "degree": + cube_spacing_m = cube_spacing_m * 1.0e5 # Scale by 100km + if crs.axis_info[0].unit_name == "degree": - out_spacing = cube_spacing / 1.0e5 # Scale by 100km - out_snwe = clip_bbox(out_snwe, out_spacing) + out_spacing = cube_spacing_m / 1.0e5 # Scale by 100km else: - out_spacing = cube_spacing - out_snwe = clip_bbox(out_snwe, out_spacing) + out_spacing = cube_spacing_m + out_snwe = clip_bbox(out_snwe, out_spacing) logger.debug(f"Output SNWE: {out_snwe}") logger.debug(f"Output cube spacing: {out_spacing}") - # Load downloaded weather model file to get projection info - with xarray.load_dataset(weather_model_file) as ds: - # Output grid points - North up grid - if args['height_levels'] is not None: - heights = args['height_levels'] - else: - heights = ds.z.values - - logger.debug(f'Output height range is {min(heights)} to {max(heights)}') - - # Load CRS from weather model file - wm_proj = rio_profile(f"netcdf:{weather_model_file}:t")["crs"] - if wm_proj is None: - print("WARNING: I can't find a CRS in the weather model file, so I will assume you are using WGS84") - wm_proj = CRS.from_epsg(4326) - else: - wm_proj = CRS.from_wkt(wm_proj.to_wkt()) - # Build the output grid zpts = np.array(heights) xpts = np.arange(out_snwe[2], out_snwe[3] + out_spacing, out_spacing) ypts = np.arange(out_snwe[1], out_snwe[0] - out_spacing, -out_spacing) + # If no orbit is provided # Build zenith delay cube - if los.is_Zenith(): - out_type = "zenith" - out_filename = wf.replace("wet", "tropo") + if los.is_Zenith() or los.is_Projected(): + out_type = ["zenith" if los.is_Zenith() else 'slant - projected'][0] # Get ZTD interpolators ifWet, ifHydro = getInterpolators(weather_model_file, "total") # Build cube - wetDelay, hydroDelay = build_cube( + wetDelay, hydroDelay = _build_cube( xpts, ypts, zpts, wm_proj, crs, [ifWet, ifHydro]) else: - out_type = "slant range" - if not los.ray_trace(): - out_filename = wf.replace("_ztd", "_std").replace("wet", "tropo") - else: - out_filename = wf.replace("_ztd", "_ray").replace("wet", "tropo") - - if args["look_dir"].lower() not in ["right", "left"]: - raise ValueError( - f'Unknown look direction: {args["look_dir"]}' - ) + out_type = "slant - raytracing" # Get pointwise interpolators ifWet, ifHydro = getInterpolators( @@ -157,25 +172,22 @@ def tropo_delay_cube(dt, wf, args, model_file=None): shared=(nproc > 1), ) - if los.ray_trace(): - # Build cube - if nproc == 1: - wetDelay, hydroDelay = build_cube_ray( - xpts, ypts, zpts, - dt, args["los"]._file, args["look_dir"], - wm_proj, crs, - [ifWet, ifHydro]) - - ### Use multi-processing here - else: - # Pre-build output arrays + # Build cube + if nproc == 1: + wetDelay, hydroDelay = _build_cube_ray( + xpts, ypts, zpts, + dt, los._file, look_dir, + wm_proj, crs, + [ifWet, ifHydro]) + + ### Use multi-processing here + else: + # Pre-build output arrays - # Create worker pool + # Create worker pool - # Loop over heights - raise NotImplementedError - else: - raise NotImplementedError('Conventional STD is not yet implemented on the cube') + # Loop over heights + raise NotImplementedError # Write output file # Modify this as needed for NISAR / other projects @@ -205,7 +217,7 @@ def tropo_delay_cube(dt, wf, args, model_file=None): source=os.path.basename(weather_model_file), history=str(datetime.datetime.utcnow()) + " RAiDER", description=f"RAiDER geo cube - {out_type}", - reference_time=str(args["time"]), + reference_time=str(dt), ), ) @@ -240,73 +252,49 @@ def tropo_delay_cube(dt, wf, args, model_file=None): ds.x.attrs["long_name"] = "x-coordinate in projected coordinate system" ds.x.attrs["units"] = "m" + return ds - ext = os.path.splitext(out_filename) - if ext not in '.nc .h5'.split(): - out_filename = f'{os.path.splitext(out_filename)[0]}.nc' - logger.debug('Invalid extension %s for cube. Defaulting to .nc', ext) - - if out_filename.endswith(".nc"): - ds.to_netcdf(out_filename, mode="w") - elif out_filename.endswith(".h5"): - ds.to_netcdf(out_filename, engine="h5netcdf", invalid_netcdf=True) - logger.info('Finished writing data to: %s', out_filename) - return - - -def checkQueryPntsFile(pnts_file, query_shape): - ''' - Check whether the query points file exists, and if it - does, check that the shapes are all consistent - ''' - write_flag = True - if os.path.exists(pnts_file): - # Check whether the number of points is consistent with the new inputs - with h5py.File(pnts_file, 'r') as f: - if query_shape == tuple(f['lon'].attrs['Shape']): - write_flag = False - - return write_flag - - -def transformPoints(lats, lons, hgts, old_proj, new_proj): +def transformPoints(lats: np.ndarray, lons: np.ndarray, hgts: np.ndarray, old_proj: CRS, new_proj: CRS) -> np.ndarray: ''' Transform lat/lon/hgt data to an array of points in a new projection - Parameters - ---------- - lats - WGS-84 latitude (EPSG: 4326) - lons - ditto for longitude - hgts - Ellipsoidal height in meters - old_proj - the original projection of the points - new_proj - the new projection in which to return the points - - Returns - ------- - the array of query points in the weather model coordinate system (YX) + Args: + lats: ndarray - WGS-84 latitude (EPSG: 4326) + lons: ndarray - ditto for longitude + hgts: ndarray - Ellipsoidal height in meters + old_proj: CRS - the original projection of the points + new_proj: CRS - the new projection in which to return the points + + Returns: + ndarray: the array of query points in the weather model coordinate system (YX) ''' t = Transformer.from_crs(old_proj, new_proj) # Flags for flipping inputs or outputs - in_flip = old_proj.axis_info[0].direction == "east" - out_flip = new_proj.axis_info[0].direction == "east" + if not isinstance(new_proj, pyproj.CRS): + new_proj = CRS.from_epsg(new_proj.lstrip('EPSG:')) + if not isinstance(old_proj, pyproj.CRS): + old_proj = CRS.from_epsg(old_proj.lstrip('EPSG:')) + + in_flip = old_proj.axis_info[0].direction + out_flip = new_proj.axis_info[0].direction - if in_flip: + if in_flip == 'east': res = t.transform(lons, lats, hgts) else: res = t.transform(lats, lons, hgts) - if out_flip: + if out_flip == 'east': return np.stack((res[1], res[0], res[2]), axis=-1).T else: return np.stack(res, axis=-1).T -def build_cube(xpts, ypts, zpts, model_crs, pts_crs, interpolators): +def _build_cube(xpts, ypts, zpts, model_crs, pts_crs, interpolators): """ - Iterate over interpolators and build a cube + Iterate over interpolators and build a cube using Zenith """ # Create a regular 2D grid xx, yy = np.meshgrid(xpts, ypts) @@ -352,10 +340,12 @@ def build_cube(xpts, ypts, zpts, model_crs, pts_crs, interpolators): return outputArrs -def build_cube_ray(xpts, ypts, zpts, ref_time, orbit_file, look_dir, model_crs, - pts_crs, interpolators, outputArrs=None): +def _build_cube_ray( + xpts, ypts, zpts, ref_time, orbit_file, look_dir, model_crs, + pts_crs, interpolators, outputArrs=None +): """ - Iterate over interpolators and build a cube + Iterate over interpolators and build a cube using raytracing """ # Some constants for this module @@ -552,158 +542,3 @@ def build_cube_ray(xpts, ypts, zpts, ref_time, orbit_file, look_dir, model_crs, if output_created_here: return outputArrs - - -############################################################################### -def main(dt, wetFilename, hydroFilename, args): - """ - raider main function for calculating delays. - - Parameterss - ---------- - args: dict Parameters and inputs needed for processing, - containing the following key-value pairs: - - los - tuple, Zenith class object, ('los', 2-band los file), or ('sv', orbit_file) - lats - ndarray or str - lons - ndarray or str - heights - see checkArgs for format - weather_model - type of weather model to use - wmLoc - Directory containing weather model files - zref - max integration height - outformat - File format to use for raster outputs - time - list of datetimes to calculate delays - download_only - Only download the raw weather model data and exit - wetFilename - - hydroFilename - - pnts_file - Input a points file from previous run - verbose - verbose printing - """ - # unpacking the dictionairy - los = args['los'] - heights = args['dem'] - weather_model = args['weather_model'] - wmLoc = args['weather_model_directory'] - zref = args['zref'] - outformat = args['raster_format'] - verbose = args['verbose'] - aoi = args['aoi'] - - download_only = args['download_only'] - - if los.ray_trace(): - ll_bounds = aoi.add_buffer(buffer=1) # add a buffer for raytracing - else: - ll_bounds = aoi.bounds() - - # logging - logger.debug('Starting to run the weather model calculation') - logger.debug('Time type: {}'.format(type(dt))) - logger.debug('Time: {}'.format(dt.strftime('%Y%m%d'))) - logger.debug('Max integration height is {:1.1f} m'.format(zref)) - - ########################################################### - # weather model calculation - logger.debug('Beginning weather model pre-processing') - - weather_model_file = prepareWeatherModel( - weather_model, - dt, - wmLoc=wmLoc, - ll_bounds=ll_bounds, # SNWE - zref=zref, - download_only=download_only, - makePlots=verbose, - ) - - if download_only: - logger.debug('Weather model has downloaded. Finished.') - return None, None - - - if aoi.type() == 'bounding_box' or \ - (args['height_levels'] and aoi.type() != 'station_file'): - # This branch is specifically for cube generation - try: - tropo_delay_cube( - dt, wetFilename, args, - model_file=weather_model_file, - ) - except Exception as e: - logger.error(e) - raise RuntimeError('Something went wrong in calculating delays on the cube') - return None, None - - ########################################################### - # Load the downloaded model file for CRS information - wm_proj = rio_profile(f"netcdf:{weather_model_file}:t")["crs"] - if wm_proj is None: - print("WARNING: I can't find a CRS in the weather model file, so I will assume you are using WGS84") - wm_proj = CRS.from_epsg(4326) - else: - wm_proj = CRS.from_wkt(wm_proj.to_wkt()) - #################################################################### - - #################################################################### - # Calculate delays - if isinstance(los, (Zenith, Conventional)): - # Start actual processing here - logger.debug("Beginning DEM calculation") - # Lats, Lons will be translated from file to array here if needed - - # read lats/lons - lats, lons = aoi.readLL() - hgts = aoi.readZ() - - los.setPoints(lats, lons, hgts) - - # Transform query points if needed - pnt_proj = CRS.from_epsg(4326) - if wm_proj != pnt_proj: - pnts = transformPoints( - lats, - lons, - hgts, - pnt_proj, - wm_proj, - ).T - else: - # interpolators require y, x, z - pnts = np.stack([lats, lons, hgts], axis=-1) - - # either way I'll need the ZTD - ifWet, ifHydro = getInterpolators(weather_model_file, 'total') - wetDelay = ifWet(pnts) - hydroDelay = ifHydro(pnts) - - # return the delays (ZTD or STD) - wetDelay = los(wetDelay) - hydroDelay = los(hydroDelay) - - elif isinstance(los, Raytracing): - raise NotImplementedError - else: - raise ValueError("Unknown operation type") - - ########################################################### - # Write the delays to file - # Different options depending on the inputs - - if not isinstance(wetFilename, str): - wetFilename = wetFilename[0] - hydroFilename = hydroFilename[0] - - if aoi.type() == 'station_file': - wetFilename = f'{os.path.splitext(wetFilename)[0]}.csv' - - writeDelays( - aoi, - wetDelay, - hydroDelay, - wetFilename, - hydroFilename, - outformat=outformat, - ) - logger.info('Finished writing data to file') - - return wetDelay, hydroDelay diff --git a/tools/RAiDER/delayFcns.py b/tools/RAiDER/delayFcns.py index ca1ebc347..60e5f1e88 100755 --- a/tools/RAiDER/delayFcns.py +++ b/tools/RAiDER/delayFcns.py @@ -12,19 +12,18 @@ from netCDF4 import Dataset import numpy as np from pyproj import CRS, Transformer -from scipy.interpolate import RegularGridInterpolator +from scipy.interpolate import RegularGridInterpolator as Interpolator from RAiDER.constants import _STEP -from RAiDER.interpolator import RegularGridInterpolator as Interpolator from RAiDER.makePoints import makePoints1D def calculate_start_points(x, y, z, ds): ''' - Parameters + Args: ---------- wm_file: str - A file containing a regularized weather model. - Returns + Returns: ------- SP: ndarray - a * x 3 array containing the XYZ locations of the pixels in ECEF coordinates. Note the ordering of the array is [Y X Z] @@ -92,18 +91,19 @@ def getInterpolators(wm_file, kind='pointwise', shared=False): an interpolator ''' # Get the weather model data - with Dataset(wm_file, mode='r') as f: - xs_wm = np.array(f.variables['x'][:]) - ys_wm = np.array(f.variables['y'][:]) - zs_wm = np.array(f.variables['z'][:]) - - # Can get the point-wise or total delays, depending on what is requested - if kind == 'pointwise': - wet = np.array(f.variables['wet'][:]).transpose(1, 2, 0) - hydro = np.array(f.variables['hydro'][:]).transpose(1, 2, 0) - elif kind == 'total': - wet = np.array(f.variables['wet_total'][:]).transpose(1, 2, 0) - hydro = np.array(f.variables['hydro_total'][:]).transpose(1, 2, 0) + try: + ds = xarray.load_dataset(wm_file) + except: + ds = wm_file + + xs_wm = np.array(ds.variables['x'][:]) + ys_wm = np.array(ds.variables['y'][:]) + zs_wm = np.array(ds.variables['z'][:]) + wet = ds.variables['wet_total' if kind=='total' else 'wet'][:] + hydro = ds.variables['hydro_total' if kind=='total' else 'hydro'][:] + + wet = np.array(wet).transpose(1, 2, 0) + hydro = np.array(hydro).transpose(1, 2, 0) # If shared interpolators are requested # The arrays are not modified - so turning off lock for performance @@ -114,9 +114,8 @@ def getInterpolators(wm_file, kind='pointwise', shared=False): wet = make_shared_raw(wet) hydro = make_shared_raw(hydro) - - ifWet = Interpolator((ys_wm, xs_wm, zs_wm), wet, fill_value=np.nan) - ifHydro = Interpolator((ys_wm, xs_wm, zs_wm), hydro, fill_value=np.nan) + ifWet = Interpolator((ys_wm, xs_wm, zs_wm), wet, fill_value=np.nan, bounds_error = False) + ifHydro = Interpolator((ys_wm, xs_wm, zs_wm), hydro, fill_value=np.nan, bounds_error = False) return ifWet, ifHydro diff --git a/tools/RAiDER/dem.py b/tools/RAiDER/dem.py index cd15939d8..0327f4bab 100644 --- a/tools/RAiDER/dem.py +++ b/tools/RAiDER/dem.py @@ -30,7 +30,7 @@ def getHeights(ll_bounds, dem_type, dem_file, lats=None, lons=None): 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 @@ -41,16 +41,13 @@ def getHeights(ll_bounds, dem_type, dem_file, lats=None, lons=None): elif (dem_type == 'download') or (dem_type == 'dem'): if ~os.path.exists(dem_file): - download_dem( - ll_bounds, - writeDEM = True, - outName=dem_file, - ) + download_dem(ll_bounds, writeDEM=True, outName=dem_file) + #TODO: interpolate heights to query lats/lons # Interpolate to the query points hts = interpolateDEM( dem_file, - lats, + lats, lons, ) @@ -59,31 +56,34 @@ def getHeights(ll_bounds, dem_type, dem_file, lats=None, lons=None): def download_dem( ll_bounds, - save_flag='new', writeDEM=False, outName='warpedDEM', - buf=0.02 + buf=0.02, + overwrite=False, ): - ''' Download a DEM if one is not already present. ''' - folder = os.path.dirname(outName) - # inExtent is SNWE - # dem-stitcher wants WSEN - bounds = [ - np.floor(ll_bounds[2]) - buf, np.floor(ll_bounds[0]) - buf, - np.ceil(ll_bounds[3]) + buf, np.ceil(ll_bounds[1]) + buf - ] + """ Download a DEM if one is not already present. """ + if os.path.exists(outName) and not overwrite: + logger.info('Using existing DEM: %s', outName) + zvals, metadata = rio_open(outName, returnProj=True) - zvals, metadata = stitch_dem( - bounds, - dem_name='glo_30', - dst_ellipsoidal_height=True, - dst_area_or_point='Area', - ) - if writeDEM: - dem_file = os.path.join(folder, 'GLO30_fullres_dem.tif') - with rasterio.open(dem_file, 'w', **metadata) as ds: - ds.write(zvals, 1) - ds.update_tags(AREA_OR_POINT='Point') + else: + # inExtent is SNWE + # dem-stitcher wants WSEN + bounds = [ + np.floor(ll_bounds[2]) - buf, np.floor(ll_bounds[0]) - buf, + np.ceil(ll_bounds[3]) + buf, np.ceil(ll_bounds[1]) + buf + ] - return zvals, metadata + zvals, metadata = stitch_dem( + bounds, + dem_name='glo_30', + dst_ellipsoidal_height=True, + dst_area_or_point='Area', + ) + if writeDEM: + with rasterio.open(outName, 'w', **metadata) as ds: + ds.write(zvals, 1) + ds.update_tags(AREA_OR_POINT='Point') + logger.info('Wrote DEM: %s', outName) + return zvals, metadata diff --git a/tools/RAiDER/gnss/downloadGNSSDelays.py b/tools/RAiDER/gnss/downloadGNSSDelays.py index 971dc778b..04f1ed283 100755 --- a/tools/RAiDER/gnss/downloadGNSSDelays.py +++ b/tools/RAiDER/gnss/downloadGNSSDelays.py @@ -5,15 +5,11 @@ # RESERVED. United States Government Sponsorship acknowledged. # # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -import argparse import itertools import multiprocessing import os import pandas as pd -from textwrap import dedent -from RAiDER.cli.parser import add_cpus, add_out, add_verbose -from RAiDER.cli.validators import DateListAction, date_type from RAiDER.logger import logger, logging from RAiDER.getStationDelays import get_station_data from RAiDER.utilFcns import requests_retry_session @@ -22,92 +18,6 @@ _UNR_URL = "http://geodesy.unr.edu/" -def create_parser(): - """Parse command line arguments using argparse.""" - p = argparse.ArgumentParser( - formatter_class=argparse.RawDescriptionHelpFormatter, - description=""" -Check for and download tropospheric zenith delays for a set of GNSS stations from UNR - -Example call to virtually access and append zenith delay information to a CSV table in specified output -directory, across specified range of time (in YYMMDD YYMMDD) and all available times of day, and confined to specified -geographic bounding box : -downloadGNSSdelay.py --out products -y 20100101 20141231 -b '39 40 -79 -78' - -Example call to virtually access and append zenith delay information to a CSV table in specified output -directory, across specified range of time (in YYMMDD YYMMDD) and specified time of day, and distributed globally : -downloadGNSSdelay.py --out products -y 20100101 20141231 --returntime '00:00:00' - - -Example call to virtually access and append zenith delay information to a CSV table in specified output -directory, across specified range of time in 12 day steps (in YYMMDD YYMMDD days) and specified time of day, and distributed globally : -downloadGNSSdelay.py --out products -y 20100101 20141231 12 --returntime '00:00:00' - -Example call to virtually access and append zenith delay information to a CSV table in specified output -directory, across specified range of time (in YYMMDD YYMMDD) and specified time of day, and distributed globally but restricted -to list of stations specified in input textfile : -downloadGNSSdelay.py --out products -y 20100101 20141231 --returntime '00:00:00' -f station_list.txt - -NOTE, following example call to physically download zenith delay information not recommended as it is not -necessary for most applications. -Example call to physically download and append zenith delay information to a CSV table in specified output -directory, across specified range of time (in YYMMDD YYMMDD) and specified time of day, and confined to specified -geographic bounding box : -downloadGNSSdelay.py --download --out products -y 20100101 20141231 --returntime '00:00:00' -b '39 40 -79 -78' -""") - - # Stations to check/download - area = p.add_argument_group( - 'Stations to check/download. Can be a lat/lon bounding box or file, or will run the whole world if not specified') - area.add_argument( - '--station_file', '-f', default=None, dest='station_file', - help=('Text file containing a list of 4-char station IDs separated by newlines')) - area.add_argument( - '-b', '--bounding_box', dest='bounding_box', type=str, default=None, - help="Provide either valid shapefile or Lat/Lon Bounding SNWE. -- Example : '19 20 -99.5 -98.5'") - area.add_argument( - '--gpsrepo', '-gr', default='UNR', dest='gps_repo', - help=('Specify GPS repository you wish to query. Currently supported archives: UNR.')) - - misc = p.add_argument_group("Run parameters") - add_out(misc) - - misc.add_argument( - '--date', dest='dateList', - help=dedent("""\ - Date to calculate delay. - Can be a single date, a list of two dates (earlier, later) with 1-day interval, or a list of two dates and interval in days (earlier, later, interval). - Example accepted formats: - YYYYMMDD or - YYYYMMDD YYYYMMDD - YYYYMMDD YYYYMMDD N - """), - nargs="+", - action=DateListAction, - type=date_type, - required=True - ) - - misc.add_argument( - '--returntime', dest='returnTime', - help="Return delays closest to this specified time. If not specified, the GPS delays for all times will be returned. Input in 'HH:MM:SS', e.g. '16:00:00'", - default=None) - - misc.add_argument( - '--download', - help='Physically download data. Note this option is not necessary to proceed with statistical analyses, as data can be handled virtually in the program.', - action='store_true', dest='download', default=False) - - add_cpus(misc) - add_verbose(misc) - - return p - - -def cmd_line_parse(iargs=None): - parser = create_parser() - return parser.parse_args(args=iargs) - def get_station_list(bbox=None, writeLoc=None, userstatList=None, name_appendix=''): ''' @@ -286,16 +196,14 @@ def get_ID(line): return stat_id, float(lat), float(lon), float(height) -def main(params=None): +def main(inps=None): """ Main workflow for querying supported GPS repositories for zenith delay information. """ - if params is not None: - inps = params + try: dateList = inps.date_list returnTime = inps.time - else: - inps = cmd_line_parse() + except: dateList = inps.dateList returnTime = inps.returnTime diff --git a/tools/RAiDER/llreader.py b/tools/RAiDER/llreader.py index 25b720b09..4b00816be 100644 --- a/tools/RAiDER/llreader.py +++ b/tools/RAiDER/llreader.py @@ -10,12 +10,14 @@ import numpy as np import pandas as pd +import rasterio from pyproj import CRS from RAiDER.dem import download_dem from RAiDER.interpolator import interpolateDEM from RAiDER.utilFcns import rio_extents, rio_open, rio_profile, rio_stats, get_file_and_band +from RAiDER.logger import logger class AOI(object): @@ -27,6 +29,10 @@ def __init__(self): self._proj = CRS.from_epsg(4326) + def type(self): + return self._type + + def bounds(self): return self._bounding_box @@ -40,7 +46,10 @@ def add_buffer(self, buffer): Check whether an extra lat/lon buffer is needed for raytracing ''' # if raytracing, add a 1-degree buffer all around - ll_bounds = self._bounding_box.copy() + try: + ll_bounds = self._bounding_box.copy() + except AttributeError: + ll_bounds = list(self._bounding_box) ll_bounds[0] = np.max([ll_bounds[0] - buffer, -90]) ll_bounds[1] = np.min([ll_bounds[1] + buffer, 90]) ll_bounds[2] = np.max([ll_bounds[2] - buffer,-180]) @@ -51,13 +60,10 @@ def add_buffer(self, buffer): class StationFile(AOI): '''Use a .csv file containing at least Lat, Lon, and optionally Hgt_m columns''' def __init__(self, station_file): - AOI.__init__(self) + super().__init__() self._filename = station_file self._bounding_box = bounds_from_csv(station_file) - - - def type(self): - return 'station_file' + self._type = 'station_file' def readLL(self): @@ -80,42 +86,43 @@ def readZ(self): class RasterRDR(AOI): - def __init__(self, lat_file, lon_file=None, hgt_file=None, convention='isce'): - AOI.__init__(self) - # allow for 2-band lat/lon raster - if (lon_file is None): - self._file = lat_file - else: - self._latfile = lat_file - self._lonfile = lon_file + def __init__(self, lat_file, lon_file=None, hgt_file=None, dem_file=None, convention='isce'): + super().__init__() + self._type = 'radar_rasters' + self._latfile = lat_file + self._lonfile = lon_file + if (self._latfile is None) and (self._lonfile is None): + raise ValueError('You need to specify a 2-band file or two single-band files') + + try: self._proj, self._bounding_box, _ = bounds_from_latlon_rasters(lat_file, lon_file) + except rasterio.errors.RasterioIOError: + raise ValueError('Could not open {}, does it exist?'.format(self._latfile)) - # keep track of the height file + # keep track of the height file, dem and convention self._hgtfile = hgt_file + self._demfile = dem_file self._convention = convention - - def type(self): - return 'radar_rasters' - - def readLL(self): - if self._latfile is not None: - return rio_open(self._latfile), rio_open(self._lonfile) - elif self._file is not None: - return rio_open(self._file) + # allow for 2-band lat/lon raster + if self._lonfile is None: + return rio_open(self._latfile) else: - raise ValueError('lat/lon files are not defined') + return rio_open(self._latfile), rio_open(self._lonfile) def readZ(self): - if self._hgtfile is not None: + if self._hgtfile is not None and os.path.exists(self._hgtfile): + logger.info('Using existing heights at: %s', self._hgtfile) return rio_open(self._hgtfile) + else: + demFile = 'GLO30_fullres_dem.tif' if self._demfile is None else self._demfile zvals, metadata = download_dem( self._bounding_box, - writeDEM = True, - outName = os.path.join('GLO30_fullres_dem.tif'), + writeDEM=True, + outName=os.path.join(demFile), ) z_bounds = get_bbox(metadata) z_out = interpolateDEM(zvals, z_bounds, self.readLL(), method='nearest') @@ -127,24 +134,19 @@ class BoundingBox(AOI): def __init__(self, bbox): AOI.__init__(self) self._bounding_box = bbox - - def type(self): - return 'bounding_box' + self._type = 'bounding_box' class GeocodedFile(AOI): '''Parse a Geocoded file for coordinates''' def __init__(self, filename, is_dem=False): - AOI.__init__(self) + super().__init__() self._filename = filename self.p = rio_profile(filename) self._bounding_box = rio_extents(self.p) self._is_dem = is_dem _, self._proj, self._gt = rio_stats(filename) - - - def type(self): - return 'geocoded_file' + self._type = 'geocoded_file' def readLL(self): @@ -160,31 +162,24 @@ def readLL(self): def readZ(self): - if self._is_dem: - return rio_open(self._filename) - - else: - zvals, metadata = download_dem( - self._bounding_box, - writeDEM = True, - outName = os.path.join('GLO30_fullres_dem.tif'), - ) - z_bounds = get_bbox(metadata) - z_out = interpolateDEM(zvals, z_bounds, self.readLL(), method='nearest') - return z_out + demFile = self._filename if self._is_dem else 'GLO30_fullres_dem.tif' + bbox = self._bounding_box + zvals, metadata = download_dem(bbox, writeDEM=True, outName=demFile) + z_bounds = get_bbox(metadata) + z_out = interpolateDEM(zvals, z_bounds, self.readLL(), method='nearest') + return z_out class Geocube(AOI): '''Parse a georeferenced data cube''' def __init__(self): - AOI.__init__(self) + super().__init__() + self._type = 'geocube' raise NotImplementedError - def type(self): - return 'geocube' def readLL(self): - raise NotImplementedError + return None, None def bounds_from_latlon_rasters(latfile, lonfile): diff --git a/tools/RAiDER/losreader.py b/tools/RAiDER/losreader.py index 98bd78b84..c8235eae1 100644 --- a/tools/RAiDER/losreader.py +++ b/tools/RAiDER/losreader.py @@ -23,13 +23,16 @@ class LOS(ABC): - '''LOS Class definition for handling look vectors''' - + ''' + LOS Class definition for handling look vectors + ''' def __init__(self): self._lats, self._lons, self._heights = None, None, None self._look_vecs = None self._ray_trace = False self._is_zenith = False + self._is_projected = False + def setPoints(self, lats, lons=None, heights=None): '''Set the pixel locations''' @@ -50,23 +53,33 @@ def setPoints(self, lats, lons=None, heights=None): self._lats = lats self._lons = lons self._heights = heights - + + def setTime(self, dt): self._time = dt - + + def is_Zenith(self): return self._is_zenith - + + + def is_Projected(self): + return self._is_projected + + def ray_trace(self): return self._ray_trace class Zenith(LOS): - """Special value indicating a look vector of "zenith".""" + """ + Class definition for a "Zenith" object. + """ def __init__(self): - LOS.__init__(self) + super().__init__() self._is_zenith = True + def setLookVectors(self): '''Set point locations and calculate Zenith look vectors''' if self._lats is None: @@ -74,6 +87,7 @@ def setLookVectors(self): if self._look_vecs is None: self._look_vecs = getZenithLookVecs(self._lats, self._lons, self._heights) + def __call__(self, delays): '''Placeholder method for consistency with the other classes''' return delays @@ -84,18 +98,21 @@ class Conventional(LOS): Special value indicating that the zenith delay will be projected using the standard cos(inc) scaling. """ - - def __init__(self, filename=None, los_convention='isce', time=None, pad=None): - LOS.__init__(self) + def __init__(self, filename=None, los_convention='isce', time=None, pad=600): + super().__init__() self._file = filename self._time = time - self._pad = pad - self._convention = los_convention - if self._convention == 'hyp3': + self._pad = pad + self._is_projected = True + self._convention = los_convention + if self._convention.lower() != 'isce': raise NotImplementedError() + def __call__(self, delays): - '''Read the LOS file and convert it to look vectors''' + ''' + Read the LOS file and convert it to look vectors + ''' if self._lats is None: raise ValueError('Target points not set') if self._file is None: @@ -104,6 +121,7 @@ def __call__(self, delays): try: # if an ISCE-style los file is passed open it with GDAL LOS_enu = inc_hd_to_enu(*rio_open(self._file)) + except OSError: # Otherwise, treat it as an orbit / statevector file svs = np.stack( @@ -133,7 +151,7 @@ class Raytracing(LOS): because they are in an ECEF reference frame instead of a local ENU. This is done because the construction of rays is done in ECEF rather than the local ENU. - Parameters + Args: ---------- time: python datetime - user-requested query time. Must be compatible with the orbit file passed. @@ -142,36 +160,34 @@ class Raytracing(LOS): the user-specified time; default 3 hours Only required for a statevector file. - Returns + Returns: ------- - look_vecs: ndarray - an x 3 array of unit look vectors, defined in - an Earth-centered, earth-fixed reference frame (ECEF). - Convention is vectors point from the target pixel to the - sensor. - lengths: ndarray - array of of the distnce from the surface to - the top of the troposphere (denoted by zref) + ndarray - an x 3 array of unit look vectors, defined in + an Earth-centered, earth-fixed reference frame (ECEF). + Convention is vectors point from the target pixel to the + sensor. + ndarray - array of of the distnce from the surface to + the top of the troposphere (denoted by zref) Example: -------- >>> from RAiDER.losreader import Raytracing >>> import numpy as np - ->>> #TODO - >>> # + >>> TODO """ - def __init__(self, filename=None, los_convention='isce', time=None, pad=None): + def __init__(self, filename=None, los_convention='isce', time=None, pad=600): '''read in and parse a statevector file''' - LOS.__init__(self) + super().__init__() self._ray_trace = True self._file = filename self._time = time - self._pad = pad + self._pad = pad self._convention = los_convention - self._pad = 600 - if self._convention == 'hyp3': + if self._convention.lower() != 'isce': raise NotImplementedError() + def getLookVectors(self, time, pad=3 * 60): ''' Calculate look vectors for raytracing @@ -210,6 +226,7 @@ def getLookVectors(self, time, pad=3 * 60): out="ecef" ) + def getIntersectionWithHeight(self, height): """ This function computes the intersection point of a ray at a height @@ -218,6 +235,7 @@ def getIntersectionWithHeight(self, height): # We just leverage the same code as finding top of atmosphere here return getTopOfAtmosphere(self._xyz, self._look_vecs, height) + def getIntersectionWithLevels(self, levels): """ This function returns the points at which rays intersect the @@ -243,6 +261,7 @@ def getIntersectionWithLevels(self, levels): return rays + def calculateDelays(self, delays): ''' Here "delays" is point-wise delays (i.e. refractivities), not @@ -259,13 +278,11 @@ def getZenithLookVecs(lats, lons, heights): ''' Returns look vectors when Zenith is used. - Parameters - ---------- - lats/lons/heights: ndarray - Numpy arrays containing WGS-84 target locations + Args: + lats/lons/heights (ndarray): - Numpy arrays containing WGS-84 target locations - Returns - ------- - zenLookVecs: ndarray - (in_shape) x 3 unit look vectors in an ECEF reference frame + Returns: + zenLookVecs (ndarray): - (in_shape) x 3 unit look vectors in an ECEF reference frame ''' x = np.cos(np.radians(lats)) * np.cos(np.radians(lons)) y = np.cos(np.radians(lats)) * np.sin(np.radians(lons)) @@ -278,20 +295,17 @@ def get_sv(los_file, ref_time, pad=3 * 60): """ Read an LOS file and return orbital state vectors - Parameters - ---------- - los_file: str - user-passed file containing either look - vectors or statevectors for the sensor - ref_time: python datetime - User-requested datetime; if not encompassed - by the orbit times will raise a ValueError - pad: int - number of seconds to keep around the - requested time - - Returns - ------- - svs: 7 x 1 list of Nt x 1 ndarrays - the times, x/y/z positions and - velocities of the sensor for the given - window around the reference time + Args: + los_file (str): - user-passed file containing either look + vectors or statevectors for the sensor + ref_time (datetime):- User-requested datetime; if not encompassed + by the orbit times will raise a ValueError + pad (int): - number of seconds to keep around the + requested time + + Returns: + svs (list of ndarrays): - the times, x/y/z positions and velocities + of the sensor for the given window around the reference time """ try: svs = read_txt_file(los_file) @@ -318,15 +332,13 @@ def inc_hd_to_enu(incidence, heading): Convert incidence and heading to line-of-sight vectors from the ground to the top of the troposphere. - Parameters - ---------- - incidence: ndarray - incidence angle in deg from vertical - heading: ndarray - heading angle in deg clockwise from north - lats/lons/heights: ndarray - WGS84 ellipsoidal target (ground pixel) locations + Args: + incidence: ndarray - incidence angle in deg from vertical + heading: ndarray - heading angle in deg clockwise from north + lats/lons/heights: ndarray - WGS84 ellipsoidal target (ground pixel) locations - Returns - ------- - LOS: ndarray - (input_shape) x 3 array of unit look vectors in local ENU + Returns: + LOS: ndarray - (input_shape) x 3 array of unit look vectors in local ENU Algorithm referenced from http://earthdef.caltech.edu/boards/4/topics/327 ''' @@ -376,16 +388,15 @@ def read_txt_file(filename): should be denoted as integer time in seconds since the reference epoch (user-requested time). - Parameters - ---------- - filename: str - user-supplied space-delimited text file with no header - containing orbital statevectors as 7 columns: - - time in seconds since the user-supplied epoch - - x / y / z locations in ECEF cartesian coordinates - - vx / vy / vz velocities in m/s in ECEF coordinates - Returns - svs: list - a length-7 list of numpy vectors containing the above - variables + Args: + filename (str): - user-supplied space-delimited text file with no header + containing orbital statevectors as 7 columns: + - time in seconds since the user-supplied epoch + - x / y / z locations in ECEF cartesian coordinates + - vx / vy / vz velocities in m/s in ECEF coordinates + Returns: + svs (list): - a length-7 list of numpy vectors containing the above + variables ''' t = list() x = list() @@ -423,11 +434,11 @@ def read_ESA_Orbit_file(filename): ''' Read orbit data from an orbit file supplied by ESA - Parameters + Args: ---------- filename: str - string of the orbit filename - Returns + Returns: ------- t: Nt x 1 ndarray - a numpy vector with Nt elements containing time in python datetime @@ -471,12 +482,12 @@ def state_to_los(svs, llh_targets, out="lookangle"): Converts information from a state vector for a satellite orbit, given in terms of position and velocity, to line-of-sight information at each (lon,lat, height) coordinate requested by the user. - Parameters + Args: ---------- svs - t, x, y, z, vx, vy, vz - time, position, and velocity in ECEF of the sensor llh_targets - lats, lons, heights - Ellipsoidal (WGS84) positions of target ground pixels - Returns + Returns: ------- LOS - * x 3 matrix of LOS unit vectors in ECEF (*not* ENU) Example: @@ -533,12 +544,12 @@ def cut_times(times, ref_time, pad=3600 * 3): Slice the orbit file around the reference aquisition time. This is done by default using a three-hour window, which for Sentinel-1 empirically works out to be roughly the largest window allowed by the orbit time. - Parameters + Args: ---------- times: Nt x 1 ndarray - Vector of orbit times as datetime ref_time: datetime - Reference time pad: int - integer time in seconds to use as padding - Returns + Returns: ------- idx: Nt x 1 logical ndarray - a mask of times within the padded request time. """ @@ -552,13 +563,13 @@ def get_radar_pos(llh, orb, out="lookangle"): ''' Calculate the coordinate of the sensor in ECEF at the time corresponding to ***. - Parameters + Args: ---------- orb: isce3.core.Orbit - Nt x 7 matrix of statevectors: [t x y z vx vy vz] llh: ndarray - position of the target in LLH out: str - either lookangle or ecef for vector - Returns + Returns: ------- if out == "lookangle" los: ndarray - Satellite incidence angle @@ -577,7 +588,7 @@ def get_radar_pos(llh, orb, out="lookangle"): # Get xyz positions of targets here targ_xyz = np.stack( - lla2ecef(llh[:, 0], llh[:, 1], llh[:, 2]), axis=-1 + lla2ecef(llh[:, 1], llh[:, 0], llh[:, 2]), axis=-1 ) # Get some isce3 constants for this inversion diff --git a/tools/RAiDER/models/hrrr.py b/tools/RAiDER/models/hrrr.py index f9e9c252a..88599bf93 100644 --- a/tools/RAiDER/models/hrrr.py +++ b/tools/RAiDER/models/hrrr.py @@ -22,7 +22,7 @@ class HRRR(WeatherModel): def __init__(self): # initialize a weather model - WeatherModel.__init__(self) + super().__init__() self._humidityType = 'q' self._model_level_type = 'pl' # Default, pressure levels are 'pl' @@ -66,9 +66,12 @@ def __init__(self): x0 = 0 y0 = 0 earth_radius = 6371229 - p1 = CRS('+proj=lcc +lat_1={lat1} +lat_2={lat2} +lat_0={lat0} +lon_0={lon0} +x_0={x0} +y_0={y0} +a={a} +b={a} +units=m +no_defs'.format(lat1=lat1, lat2=lat2, lat0=lat0, lon0=lon0, x0=x0, y0=y0, a=earth_radius)) + p1 = CRS(f'+proj=lcc +lat_1={lat1} +lat_2={lat2} +lat_0={lat0} '\ + f'+lon_0={lon0} +x_0={x0} +y_0={y0} +a={earth_radius} '\ + f'+b={earth_radius} +units=m +no_defs') self._proj = p1 + def _fetch(self, out): ''' Fetch weather model data from HRRR @@ -88,7 +91,11 @@ def load_weather(self, *args, filename=None, **kwargs): filename = self.files # read data from grib file - ds = xarray.open_dataset(filename) + try: + ds = xarray.open_dataset(filename, engine='cfgrib') + except EOFError: + ds = xarray.open_dataset(filename, engine='netcdf4') + pl = np.array([self._convertmb2Pa(p) for p in ds.levels.values]) xArr = ds['x'].values yArr = ds['y'].values @@ -135,7 +142,7 @@ def _makeDataCubes(self, filename, out=None): # Get profile information from gdal prof = rio_profile(str(filename)) - + # Now get bounds S, N, W, E = self._ll_bounds @@ -175,15 +182,6 @@ def _makeDataCubes(self, filename, out=None): lats = lats[y_min:y_max, x_min:x_max] lons = lons[y_min:y_max, x_min:x_max] - # TODO: Remove this block once we know we handle - # different flavors of HRRR correctly - # self._proj currently used but these are available - # as attrs of ds["t"] for example - # llhtolcc = Transformer.from_crs(4326, self._proj) - # res = llhtolcc.transform(lats, lons) - # print("ERROR in X: ", np.abs(res[0] - xArr[None, :]).max()) - # print("ERROR in Y: ", np.abs(res[1] - yArr[:, None]).max()) - # Data variables t = ds['t'][:, y_min:y_max, x_min:x_max].to_numpy() z = ds['gh'][:, y_min:y_max, x_min:x_max].to_numpy() @@ -257,7 +255,8 @@ def _makeDataCubes(self, filename, out=None): for k, v in self._proj.to_cf().items(): ds_new.proj.attrs[k] = v - ds_new.to_netcdf(out) + ds_new.to_netcdf(out, engine='netcdf4') + def _download_hrrr_file(self, DATE, out, model='hrrr', product='prs', fxx=0, verbose=False): ''' diff --git a/tools/RAiDER/models/weatherModel.py b/tools/RAiDER/models/weatherModel.py index 3f0ec1cf7..ef9e31e6f 100755 --- a/tools/RAiDER/models/weatherModel.py +++ b/tools/RAiDER/models/weatherModel.py @@ -127,7 +127,7 @@ def fetch(self, out, ll_bounds, time): Checks the input datetime against the valid date range for the model and then calls the model _fetch routine - Parameters + Args: ---------- out - ll_bounds - 4 x 1 array, SNWE @@ -198,7 +198,7 @@ def load( Calls the load_weather method. Each model class should define a load_weather method appropriate for that class. 'args' should be one or more filenames. ''' - self.set_latlon_bounds(ll_bounds, Nextra=0) + self.set_latlon_bounds(ll_bounds) # If the weather file has already been processed, do nothing self._out_name = self.out_file(outLoc) @@ -381,7 +381,7 @@ def bbox(self) -> list: """ Obtains the bounding box of the weather model in lat/lon CRS. - Returns + Returns: ------- list xmin, ymin, xmax, ymax @@ -417,7 +417,7 @@ def checkContainment(self: weatherModel, Checks containment of weather model bbox of outLats and outLons provided. - Parameters + Args: ---------- weather_model : weatherModel outLats : np.ndarray @@ -429,7 +429,7 @@ def checkContainment(self: weatherModel, this ensures that translates have some overlap. The default is 1e-5 or ~11.1 meters. - Returns + Returns: ------- bool True if weather model contains bounding box of OutLats and outLons @@ -538,7 +538,7 @@ def _calculategeoh(self, z, lnsp): def getProjection(self): ''' - Returns the native weather projection, which should be a pyproj object + Returns: the native weather projection, which should be a pyproj object ''' return self._proj diff --git a/tools/RAiDER/prepFromAria.py b/tools/RAiDER/prepFromAria.py deleted file mode 100644 index a4a3f1712..000000000 --- a/tools/RAiDER/prepFromAria.py +++ /dev/null @@ -1,107 +0,0 @@ -# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -# -# Author: Jeremy Maurer -# Copyright 2020, by the California Institute of Technology. ALL RIGHTS -# RESERVED. United States Government Sponsorship acknowledged. -# -# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -import numpy as np -import rasterio -from RAiDER.utilFcns import rio_open, writeArrayToRaster - - -def parse_args(): - """Parse command line arguments using argparse.""" - import argparse - p = argparse.ArgumentParser( - description='Prepare files from ARIA-tools output to use with RAiDER') - - # Line of sight - p.add_argument( - '--incFile', '-i', type=str, - help='GDAL-readable raster image file of inclination', - metavar='INC', required=True) - p.add_argument( - '--azFile', '-a', type=str, - help='GDAL-readable raster image file of azimuth', - metavar='AZ', required=True) - - p.add_argument( - '--los_filename', '-f', default='los.rdr', type=str, dest='los_file', - help=('Output Line-of-sight filename')) - p.add_argument( - '--lat_filename', '-l', default='lat.rdr', type=str, dest='lat_file', - help=('Output latitude filename')) - p.add_argument( - '--lon_filename', '-L', default='lon.rdr', type=str, dest='lon_file', - help=('Output longitude filename')) - - p.add_argument( - '--format', '-t', default='ENVI', type=str, dest='fmt', - help=('Output file format')) - - return p.parse_args(), p - - -def makeLatLonGrid(inFile, lonFileName, latFileName, fmt='ENVI'): - ''' - Convert the geocoded grids to lat/lon files for input to RAiDER - ''' - ds = rasterio.open(inFile) - xSize = ds.width - ySize = ds.height - gt = ds.transform.to_gdal() - proj = ds.crs - - # Create the xy grid - xStart = gt[0] - yStart = gt[3] - xStep = gt[1] - yStep = gt[-1] - - xEnd = xStart + xStep * xSize - 0.5 * xStep - yEnd = yStart + yStep * ySize - 0.5 * yStep - - x = np.arange(xStart, xEnd, xStep) - y = np.arange(yStart, yEnd, yStep) - X, Y = np.meshgrid(x, y) - writeArrayToRaster(X, lonFileName, 0., fmt, proj, gt) - writeArrayToRaster(Y, latFileName, 0., fmt, proj, gt) - - -def makeLOSFile(incFile, azFile, fmt='ENVI', filename='los.rdr'): - ''' - Create a line-of-sight file from ARIA-derived azimuth and inclination files - ''' - az, az_prof = rio_open(azFile, returnProj=True) - az[az == 0] = np.nan - inc = rio_open(incFile) - - heading = 90 - az - heading[np.isnan(heading)] = 0. - - array_shp = np.shape(az)[:2] - - # Write the data to a file - with rasterio.open(filename, mode="w", count=2, - driver=fmt, width=array_shp[1], - height=array_shp[0], crs=az_prof.crs, - transform=az_prof.transform, - dtype=az.dtype, nodata=0.) as dst: - dst.write(inc, 1) - dst.write(heading, 2) - - return 0 - - -def prepFromAria(): - ''' - A command-line utility to convert ARIA standard product outputs from ARIA-tools to - RAiDER-compatible format - ''' - args, p = parse_args() - makeLOSFile(args.incFile, args.azFile, args.fmt, args.los_file) - makeLatLonGrid(args.incFile, args.lon_file, args.lat_file, args.fmt) - -def main(): - prepFromAria() diff --git a/tools/RAiDER/processWM.py b/tools/RAiDER/processWM.py index ac944fa06..11c6b9a25 100755 --- a/tools/RAiDER/processWM.py +++ b/tools/RAiDER/processWM.py @@ -6,25 +6,39 @@ # # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ import os + +import matplotlib.pyplot as plt import numpy as np + +from typing import List + from RAiDER.logger import logger from RAiDER.utilFcns import getTimeFromFile -import matplotlib.pyplot as plt def prepareWeatherModel( - weather_model, - time=None, - wmLoc=None, - ll_bounds=None, - zref=None, - download_only=False, - makePlots=False, - force_download=False, -): - ''' - Parse inputs to download and prepare a weather model grid for interpolation - ''' + weather_model, + time=None, + wmLoc: str=None, + ll_bounds: List[float]=None, + download_only: bool=False, + makePlots: bool=False, + force_download: bool=False, + ) -> str: + """Parse inputs to download and prepare a weather model grid for interpolation + + Args: + weather_model: WeatherModel - instantiated weather model object + time: datetime - Python datetime to request. Will be rounded to nearest available time + wmLoc: str - file path to which to write weather model file(s) + ll_bounds: list of float - bounding box to download in [S, N, W, E] format + download_only: bool - False if preprocessing weather model data + makePlots: bool - whether to write debug plots + force_download: bool - True if you want to download even when the weather model exists + + Returns: + str: filename of the netcdf file to which the weather model has been written + """ # Ensure the file output location exists if wmLoc is None: wmLoc = os.path.join(os.getcwd(), 'weather_files') @@ -73,7 +87,6 @@ def prepareWeatherModel( f = weather_model.load( wmLoc, ll_bounds = ll_bounds, - zref=zref, ) if f is not None: logger.warning( @@ -122,33 +135,17 @@ def prepareWeatherModel( del weather_model -def checkBounds(weather_model, outLats, outLons): - '''Check the bounds of a weather model''' - ds = xr.load_dataset(weather_model.files[0]) # TODO: xr is undefined - coords = ds.coords # coords is dict-like - keys = [k for k in coords.keys()] - xc = coords[keys[0]] - yc = coords[keys[1]] - lat_bounds = [yc.min(), yc.max()] - lon_bounds = [xc.min(), xc.max()] - self_extent = lat_bounds + lon_bounds - in_extent = weather_model._getExtent(outLats, outLons) - - return in_extent, self_extent - - -def weather_model_debug( - los, - lats, - lons, - ll_bounds, - weather_model, - wmLoc, - zref, - time, - out, - download_only -): +def _weather_model_debug( + los, + lats, + lons, + ll_bounds, + weather_model, + wmLoc, + time, + out, + download_only + ): """ raiderWeatherModelDebug main function. """ @@ -180,7 +177,6 @@ def weather_model_debug( lats=lats, lons=lons, ll_bounds=ll_bounds, - zref=zref, download_only=download_only, makePlots=True ) diff --git a/tools/RAiDER/utilFcns.py b/tools/RAiDER/utilFcns.py index e5cfb86f2..bf353884d 100755 --- a/tools/RAiDER/utilFcns.py +++ b/tools/RAiDER/utilFcns.py @@ -67,7 +67,7 @@ def enu2ecef( h0: ndarray, ): """ - Parameters + Args: ---------- e1 : float target east ENU coordinate (meters) @@ -262,11 +262,15 @@ def writeArrayToRaster(array, filename, noDataValue=0., fmt='ENVI', proj=None, g if gt is not None: trans = rasterio.Affine.from_gdal(*gt) + ## cant write netcdfs with rasterio in a simple way + driver = fmt if not fmt == 'nc' else 'GTiff' with rasterio.open(filename, mode="w", count=1, width=array_shp[1], height=array_shp[0], dtype=dtype, crs=proj, nodata=noDataValue, - driver=fmt, transform=trans) as dst: + driver=driver, transform=trans) as dst: dst.write(array, 1) + logger.info('Wrote: %s', filename) + return def writeArrayToFile(lats, lons, array, filename, noDataValue=-9999): @@ -341,7 +345,7 @@ def _get_g_ll(lats): def _get_Re(lats): ''' - Returns the ellipsoid as a fcn of latitude + Returns: the ellipsoid as a fcn of latitude ''' # TODO: verify constants, add to base class constants? return np.sqrt(1 / (((cosd(lats)**2) / Rmax**2) + ((sind(lats)**2) / Rmin**2))) @@ -671,7 +675,7 @@ def UTM_to_WGS84(z, l, x, y): def transform_bbox(wesn, dest_crs=4326, src_crs=4326, margin=100.): """ Transform bbox to lat/lon or another CRS for use with rest of workflow - Returns SNWE + Returns: SNWE """ # TODO - Handle dateline crossing if isinstance(src_crs, int): @@ -991,7 +995,7 @@ def calcgeoh(lnsp, t, q, z, a, b, R_d, num_levels): Calculate pressure, geopotential, and geopotential height from the surface pressure and model levels provided by a weather model. The model levels are numbered from the highest eleveation to the lowest. - Parameters + Args: ---------- lnsp: ndarray - [y, x] array of log surface pressure t: ndarray - [z, y, x] cube of temperatures @@ -1000,7 +1004,7 @@ def calcgeoh(lnsp, t, q, z, a, b, R_d, num_levels): a: ndarray - [z] vector of a values b: ndarray - [z] vector of b values num_levels: int - integer number of model levels - Returns + Returns: ------- geopotential - The geopotential in units of height times acceleration pressurelvs - The pressure at each of the model levels for each of