diff --git a/.bumpversion.cfg b/.bumpversion.cfg deleted file mode 100644 index af1c1d0f..00000000 --- a/.bumpversion.cfg +++ /dev/null @@ -1,3 +0,0 @@ -[bumpversion] -current_version = 0.5.2 -files = pyproject.toml aurora/__init__.py diff --git a/.github/workflows/publish.yaml b/.github/workflows/publish.yaml new file mode 100644 index 00000000..69623b05 --- /dev/null +++ b/.github/workflows/publish.yaml @@ -0,0 +1,139 @@ +name: Publish on PR merge + +on: + pull_request: + types: [opened, synchronize, reopened, labeled, closed] + branches: + - main + +permissions: + contents: write + id-token: write + pull-requests: write + +jobs: + update-version-and-changelog: + # Run only when a release label is added (not on every push) + if: | + github.event.action == 'labeled' && + github.event.pull_request.merged == false && ( + contains(github.event.pull_request.labels.*.name, 'release:major') || + contains(github.event.pull_request.labels.*.name, 'release:minor') || + contains(github.event.pull_request.labels.*.name, 'release:patch') + ) + runs-on: ubuntu-latest + + steps: + - name: Check out PR branch + uses: actions/checkout@v4 + with: + ref: ${{ github.head_ref }} + fetch-depth: 0 + token: ${{ secrets.GITHUB_TOKEN }} + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: "3.13" + + - name: Install tools + run: | + pip install bump-my-version build pandoc + sudo apt-get update + sudo apt-get install -y pandoc + cargo install git-cliff + + - name: Set up Git user + run: | + git config --local user.email "github-actions[bot]@users.noreply.github.com" + git config --local user.name "github-actions[bot]" + + - name: Determine version bump type + id: bump + run: | + if [[ "${{ contains(github.event.pull_request.labels.*.name, 'release:major') }}" == "true" ]]; then + echo "type=major" >> $GITHUB_OUTPUT + elif [[ "${{ contains(github.event.pull_request.labels.*.name, 'release:minor') }}" == "true" ]]; then + echo "type=minor" >> $GITHUB_OUTPUT + else + echo "type=patch" >> $GITHUB_OUTPUT + fi + + - name: Bump version + run: bump-my-version bump ${{ steps.bump.outputs.type }} + env: + BMV_ALLOW_DIRTY: "true" + + - name: Extract version + id: version + run: | + VERSION=$(bump-my-version show current) + echo "Extracted version: $VERSION" + echo "version=$VERSION" >> $GITHUB_OUTPUT + + - name: Generate changelog + run: | + git cliff -o CHANGELOG.md + echo "Generated changelog for version ${{ steps.version.outputs.version }}" + + + - name: Commit and push changes to PR branch + run: | + git add . + git commit -m "chore: bump version to ${{ steps.version.outputs.version }} and update changelog [skip ci]" || echo "No changes to commit" + git push origin HEAD:${{ github.head_ref }} + + publish-after-merge: + # Run only after PR is merged with release labels + if: | + github.event.pull_request.merged == true && ( + contains(github.event.pull_request.labels.*.name, 'release:major') || + contains(github.event.pull_request.labels.*.name, 'release:minor') || + contains(github.event.pull_request.labels.*.name, 'release:patch') + ) + runs-on: ubuntu-latest + + steps: + - name: Check out repository + uses: actions/checkout@v4 + with: + fetch-depth: 0 + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: "3.12" + + - name: Install tools + run: pip install bump-my-version build + + - name: Extract version + id: version + run: | + VERSION=$(bump-my-version show current) + echo "version=$VERSION" >> $GITHUB_OUTPUT + + - name: Create GitHub Release + uses: softprops/action-gh-release@v2 + with: + tag_name: v${{ steps.version.outputs.version }} + name: v${{ steps.version.outputs.version }} + body_path: CHANGELOG.md + generate_release_notes: false + + - name: Build package + run: python -m build + + - name: Publish to TestPyPI + id: test-pypi + uses: pypa/gh-action-pypi-publish@v1.13.0 + with: + repository-url: https://test.pypi.org/legacy/ + skip-existing: true + attestations: false + + - name: Publish to PyPI + if: steps.test-pypi.outcome == 'success' + uses: pypa/gh-action-pypi-publish@v1.13.0 + with: + skip-existing: true \ No newline at end of file diff --git a/.github/workflows/test_notebooks.yaml b/.github/workflows/test_notebooks.yaml new file mode 100644 index 00000000..ff3f0e90 --- /dev/null +++ b/.github/workflows/test_notebooks.yaml @@ -0,0 +1,136 @@ +# Notebook Testing Workflow +# +# This workflow tests Jupyter notebooks and runs by default on: +# - Pull requests targeting the 'main' branch +# +# To manually trigger this workflow on any branch: +# 1. Go to Actions tab in GitHub +# 2. Select "Test Notebooks" workflow +# 3. Click "Run workflow" dropdown +# 4. Select your branch and click "Run workflow" + +name: Test Notebooks + +on: + pull_request: + branches: + - main + # Allow manual triggering from the Actions tab + workflow_dispatch: + +jobs: + test-notebooks: + name: Notebooks (${{ matrix.python-version }}, ${{ matrix.os }}) + runs-on: ${{ matrix.os }} + defaults: + run: + shell: bash + strategy: + fail-fast: true + matrix: + os: ["ubuntu-latest"] + python-version: ["3.10", "3.11", "3.12"] + + steps: + - uses: actions/checkout@v4 + + - name: Install uv + uses: astral-sh/setup-uv@v3 + with: + version: "latest" + + - name: Set up Python ${{ matrix.python-version }} + run: uv python install ${{ matrix.python-version }} + + - name: Cache MTH5 test files + uses: actions/cache@v4 + with: + path: ~/.cache/aurora + key: mth5-test-files-${{ runner.os }}-${{ hashFiles('tests/conftest.py') }} + restore-keys: | + mth5-test-files-${{ runner.os }}- + + - name: Create virtual environment and install dependencies + run: | + uv venv --python ${{ matrix.python-version }} + source .venv/bin/activate + uv pip install --upgrade pip + uv pip install 'setuptools<68' + uv pip install -e ".[dev,test]" + uv pip install mt_metadata[obspy] + uv pip install "mt_metadata[obspy]" + uv pip install mth5 + uv pip install git+https://github.com/kujaku11/mth5_test_data.git + # Explicitly include nbconvert & ipykernel + uv pip install jupyter nbconvert nbformat ipykernel + python -m ipykernel install --user --name "python3" + + - name: Install system dependencies + run: | + sudo apt-get update + sudo apt-get install -y pandoc + + - name: Execute Jupyter Notebooks + run: | + source .venv/bin/activate + # debugging: print locations of key packages + python -c "import pkg_resources; print('pkg_resources location:', pkg_resources.__file__)" + python -c "import setuptools; print('setuptools location:', setuptools.__file__)" + + python << 'EOF' + import nbformat + import subprocess + import sys + + notebooks = [ + "docs/examples/dataset_definition.ipynb", + "docs/examples/operate_aurora.ipynb", + "docs/tutorials/pkd_units_check.ipynb", + "docs/tutorials/pole_zero_fitting/lemi_pole_zero_fitting_example.ipynb", + "docs/tutorials/processing_configuration.ipynb", + "docs/tutorials/process_cas04_multiple_station.ipynb", + "docs/tutorials/synthetic_data_processing.ipynb" + ] + + failures = [] + + for nb_path in notebooks: + # Update kernel spec + print(f"Updating kernel in {nb_path}") + try: + with open(nb_path, "r", encoding="utf-8") as f: + nb = nbformat.read(f, as_version=4) + + nb["metadata"]["kernelspec"]["name"] = "python3" + nb["metadata"]["kernelspec"]["display_name"] = "Python (python3)" + + with open(nb_path, "w", encoding="utf-8") as f: + nbformat.write(nb, f) + print(f"✓ Updated kernel in {nb_path}") + except Exception as e: + print(f"✗ Failed to update kernel in {nb_path}: {e}") + failures.append(nb_path) + continue + + # Execute notebook + print(f"Executing {nb_path}") + result = subprocess.run( + ["jupyter", "nbconvert", "--to", "notebook", "--execute", nb_path], + capture_output=True, + text=True + ) + + if result.returncode != 0: + print(f"✗ Failed to execute {nb_path}") + print(result.stderr) + failures.append(nb_path) + else: + print(f"✓ Successfully executed {nb_path}") + + if failures: + print("\n======= Summary =======") + print(f"Failed notebooks: {failures}") + sys.exit(1) + else: + print("\n✓ All notebooks executed successfully!") + EOF diff --git a/.github/workflows/tests.yaml b/.github/workflows/tests.yaml index 1255ed2e..da755086 100644 --- a/.github/workflows/tests.yaml +++ b/.github/workflows/tests.yaml @@ -17,7 +17,7 @@ jobs: run: shell: bash strategy: - fail-fast: false + fail-fast: true matrix: os: ["ubuntu-latest"] python-version: ["3.10", "3.11", "3.12"] @@ -46,10 +46,16 @@ jobs: uv venv --python ${{ matrix.python-version }} source .venv/bin/activate uv pip install --upgrade pip + uv pip install 'setuptools<68' uv pip install -e ".[dev,test]" - # uv pip install mt_metadata[obspy] - uv pip install "mt_metadata[obspy] @ git+https://github.com/kujaku11/mt_metadata.git@patches" - uv pip install git+https://github.com/kujaku11/mth5.git@patches + uv pip install mt_metadata[obspy] + # uv pip install "mt_metadata[obspy] @ git+https://github.com/kujaku11/mt_metadata.git@patches" + # uv pip install git+https://github.com/kujaku11/mth5.git@patches + # uv pip install "mt_metadata[obspy] @ git+https://github.com/kujaku11/mt_metadata.git" + # uv pip install git+https://github.com/kujaku11/mth5.git + uv pip install "mt_metadata[obspy]" + uv pip install mth5 + # uv pip install mth5 uv pip install git+https://github.com/kujaku11/mth5_test_data.git @@ -62,67 +68,6 @@ jobs: sudo apt-get update sudo apt-get install -y pandoc - - name: Set kernel and execute Jupyter Notebooks - run: | - source .venv/bin/activate - python << 'EOF' - import nbformat - import subprocess - import sys - - notebooks = [ - "docs/examples/dataset_definition.ipynb", - "docs/examples/operate_aurora.ipynb", - "docs/tutorials/pkd_units_check.ipynb", - "docs/tutorials/pole_zero_fitting/lemi_pole_zero_fitting_example.ipynb", - "docs/tutorials/processing_configuration.ipynb", - "docs/tutorials/process_cas04_multiple_station.ipynb", - "docs/tutorials/synthetic_data_processing.ipynb" - ] - - failures = [] - - for nb_path in notebooks: - # Update kernel spec - print(f"Updating kernel in {nb_path}") - try: - with open(nb_path, "r", encoding="utf-8") as f: - nb = nbformat.read(f, as_version=4) - - nb["metadata"]["kernelspec"]["name"] = "python3" - nb["metadata"]["kernelspec"]["display_name"] = "Python (python3)" - - with open(nb_path, "w", encoding="utf-8") as f: - nbformat.write(nb, f) - print(f"✓ Updated kernel in {nb_path}") - except Exception as e: - print(f"✗ Failed to update kernel in {nb_path}: {e}") - failures.append(nb_path) - continue - - # Execute notebook - print(f"Executing {nb_path}") - result = subprocess.run( - ["jupyter", "nbconvert", "--to", "notebook", "--execute", nb_path], - capture_output=True, - text=True - ) - - if result.returncode != 0: - print(f"✗ Failed to execute {nb_path}") - print(result.stderr) - failures.append(nb_path) - else: - print(f"✓ Successfully executed {nb_path}") - - if failures: - print("\n======= Summary =======") - print(f"Failed notebooks: {failures}") - sys.exit(1) - else: - print("\n✓ All notebooks executed successfully!") - EOF - - name: Run Tests run: | source .venv/bin/activate diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 00000000..7279e6f7 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,400 @@ +# Changelog + +All notable changes to this project will be documented in this file. + +## Bug Fixes + +- Fix filter additons to use new add_filter method ([56ab23a](56ab23a18023e01be92e5af5126351236c7dcdd2)) + +- Fixing bugs with feature weighting ([438ba49](438ba4971742e7a80d43c5a6d6143ee22744f699)) + +- Fix chained assignment warnings ([321f488](321f4884630f6f3ea7b86e2ed44c6b293c65fe89)) + +- Fix FutureWarnings ([f66ddf0](f66ddf0f1259cb8a5428a3661c163f030bf2141d)) + +- Fix future warning again ([937b659](937b6599127421483f872d993e4c08191a90a0cc)) + +- Fix invalid escape sequence ([b3dc583](b3dc58323716b3cd4380f055c65f8de519ff7638)) + +- Fix indentation in docstring ([c349ea3](c349ea339cc86a3ef40a68f0d98f9c14f41e489f)) + +- Fix docstrings ([f95ce81](f95ce816630eda4dabc720cb3fd1a09bbd3f0faf)) + +- Fix docstring format ([2ba97d7](2ba97d77e331bfddd95d0350a70d5348c6995253)) + +- Fix sphinx warnings ([de0d380](de0d380dc8bdf7a039571efa94cda934c4bd4c8e)) + +- Fix sphinx warnings ([620b6be](620b6be34fa7ab3ad17c395f8ec2dbb261905330)) + +- Fix sphinx warnings ([e0d2a31](e0d2a31bc4baf1765cdfb5d9452e261f449fc46e)) + +- Fix sphinx warnings ([ee86285](ee8628555cefed6b33ccc63b6eee4055de4cdf5a)) + +- Fix warning Title level inconsitent ([38fbec3](38fbec3048a528f7eaad20595b557f38677441a9)) + +- Fix escape sequence warning ([6cd2115](6cd2115647e9773ae5e6e0636586eafbe7668c8c)) + +- Fix warning Title level inconsitent ([cf908a9](cf908a9b4d93a394c07b26dd321dac905cce94a6)) + +- Fix \c warning ([4a7808f](4a7808f8c8f1666c7ec66caadb381fbebd976088)) + + +## Revert + +- Revert parkfield paths to .cache ([d610149](d610149a215ade7db17317082b0b4e4a9eec9664)) + + +## Testing + +- Tests.yaml points at patches versions of mt_metadata and mth5 ([f4cb06e](f4cb06e0c0c5f7f706832256fb53b9c8f94c3779)) + + +## Housekeeping + +- Move method to tests ([dfe214a](dfe214ae8a386cc08524a830f64459512218e155)) + + +## Bug Fixes + +- Fix doi ([5a8b4c1](5a8b4c198deb7a5c19fa88bcc20635f59cc4534a)) + +- Fix failing notebook + +- comment out work in progress ([6118645](611864514b94d97ded9b59a78848d1c987b2dcfe)) + +- Fix typehint ([49b1f33](49b1f33d03dc073bcefdb036aa9b3d6eec5f25fa)) + +- Fix typo ([5a630e2](5a630e2f0225e6f2e5ce6cad6bd686e9928b0818)) + +- Fixing mergeing errors ([c1f96e8](c1f96e8da642917fa86a3c572490cc21563ab5ea)) + +- Fix issue #381 ([818b959](818b9591e7af3eb7b33ee877cb68ddb8bf401777)) + +- Fix assertion ([913e3d6](913e3d6f83791037fffa32dc95d57404b639110d)) + +- Fix weigths assignment ([b80e7fe](b80e7fefedc5418dd60dc343ba333a895ceae285)) + +- Fix pathing in test ([8d52a04](8d52a048afe965dc8cd0999b2b8e72736da7df65)) + +- Fix broken notebook + +- update interact arguments in from_fdsn_client (that will be deprecated) ([20752e7](20752e769c74fc73d14ffb4beaffc27dbe46c7e5)) + +- Fix typo ([7446ab9](7446ab9ed2247c3444309cf8880ced23f809d9ed)) + +- Fix issue #385 + +- this is a workaround, the change i made should be in mt_metadata ([f7d012b](f7d012b51324b024d283e27421eeeaf309bf3957)) + + +## Documentation + +- Doc and typehints ([c81729c](c81729c27f089c8adf6fe4708c7fe7cbb0e05900)) + +- Doc / type hints ([080709c](080709cd81f218e878a1011afdeeebc4da8e1f88)) + + +## Testing + +- Test new workflow with mtpy-v2 ([9b23539](9b235392a249837f00fe04437f571400f2628397)) + +- Test new workflow with mtpy-v2 ([741004c](741004c4276792132714ebb9f2d6ee99c46df4ec)) + +- Test new workflow with mtpy-v2 ([fc4b4e0](fc4b4e0dc445d62795e6c9a3a4bc17c6aceecacc)) + +- Test new workflow with mtpy-v2 ([c0ce42a](c0ce42a936f25e480d9d9118bc1ac9892d0500e6)) + + +## Bug Fixes + +- Fix typos ([4db8c40](4db8c40e3f694e52200aa57021efaa3a8ed398e8)) + +- Fix some typos, add a TODO ([50235c8](50235c8e8016d01b47e75131779d087564eef72a)) + +- Fixing station_id to station ([bb321d8](bb321d80d519a175ba75d33bc260d7f591b1e014)) + +- Fix some typos ([d048b47](d048b47eec54d05a2060b7ed83e31fc59cfc2aa3)) + +- Fix references and affiliations ([6c48e99](6c48e99c7a1714dbd5f5f6fd69f344d5d6bf9a45)) + +- Fix typo ([8676399](86763994ca9fe85f1039532cf79b904a4979d222)) + +- Fix argument name + +- was incorrectly called figure_path, not figures_path +- never noticed before because it was a kwarg ([fed4894](fed4894f9bc3cb24a145c7f3c447692d2bfc3767)) + +- Fix sphinx warning blank line end of list ([3d1c803](3d1c8033ff421544918663f18d164737ab0f6aab)) + +- Fix docstring warning ([55b12b9](55b12b9adefacfd5ea6d228c1ce50546bc1548f5)) + +- Fix sphinx warning ([50d1eff](50d1effe87bf24867487bf95db39bf8376eb5022)) + +- Fix sphinx warning ([109d5b3](109d5b3baf746ce54b1589ce10bd46207eb1f8cb)) + +- Fix sphinx warning ([bc9eb0d](bc9eb0d5ef855349092c9ad723073379d17e9c44)) + +- Fix docstrings ([7998e3f](7998e3f98bb9fe84e4befc47250db67402795ab7)) + +- Fix sphinx warning ([c3581a3](c3581a38bbba2a453f730aa98faa6e0ad42573c9)) + +- Fix sphinx warning ([3c91036](3c910363647ae18a361beafc1fafabd56c3eade8)) + +- Fix sphinx warning ([f3a132a](f3a132a2730d54283c4c4886e2ed3ee3dd1ff57d)) + +- Fix sphinx errors ([ae9bcef](ae9bcefe231a919b52983105dc7f94c0538fbac8)) + +- Fix sphinx errors/warnings ([46e8a53](46e8a535d50064127c8ef4f017bd2ce9e7b1ce13)) + +- Fix sphinx warning ([35ed995](35ed99522ff0f0194b3e58fbc7b674bf6c519da1)) + +- Fix figure named for sphinx docs ([f7f8a53](f7f8a53ff33b29e21989dadfdd189c42547b8327)) + +- Fix some sphinx warnings ([7953204](795320407f8cd89b5fbafddfa575785004b717d6)) + +- Fix typo in figure caption (#333) ([3ca63e2](3ca63e26ccceafed3dcd62d363f750604f790f69)) + + +## Documentation + +- Doc #337 ([f4e925a](f4e925a92e9ba66ff0edf5b0dfff575af5baf941)) + +- Doc #337 ([a620d32](a620d32ee8ee0a71e051b783b5e104ef60255880)) + +- Doc #337 ([4a4bf07](4a4bf073800fe380f2c68e1446f8f9b1ccabd370)) + + +## Testing + +- Tested locally ([c33ab01](c33ab014bae2e140ba076d3549ec7b2496c33ba6)) + + +## Bug Fixes + +- Fix typo ([3c8aacc](3c8aaccaba1e8f4273a340e9d2dae6957490a014)) + +- Fix indentation ([f5245ae](f5245aee7a452202ac8160e5f9abc01a809a60a2)) + +- Fix missing __init__ + +- Without this PyPI won't add needed file io_helpers/emtf_band_setup.py ([e0e6437](e0e6437be9d6e096fda4f60b56cbeed6907963ee)) + +- Fix typo in filename ([e306109](e3061092472f96744f19c1379aae34610912b7e7)) + +- Fix typo in filename ([913de9f](913de9f49e54714323b68464de2cc98aea1fb736)) + +- Fix typo in filename ([6838608](683860854a109378974b5464dc4b4df32d33b0fc)) + +- Fix typo ([c807e4a](c807e4af5b34fe44e50fb99044f0b43bb4570690)) + +- Fix paths ([bed6e99](bed6e9939609477a892fdbc55911b2103b69ff53)) + + +## Revert + +- Revert to main branches for release ([9fe5610](9fe56101e9c53d8cea8c52be66d318006bd5be44)) + + +## Testing + +- Test all 4 python versions ([46cf89e](46cf89e5a173f52eb7c9272bc9dd2e3640e8f3e6)) + +- Test can use loguru + +- since mth5 already has loguru we can just use that +- also fix one minor bug with variable name change on nan_to_mean ([d411c11](d411c1150e2807e9d2701574cdc495194fab3aae)) + +- Test precommit on new machine ([571369e](571369eda1a4e56b5440504b3599baaf106c4b4b)) + +- Test precommit on new machine ([b74f589](b74f589733f4260c24e04b41cf38c516408145b7)) + + +## Bug Fixes + +- Fixed issue of double runs when only one given ([345523e](345523ef0159443c35732453283f1a9b8394c134)) + +- Fix bug, reference station id was being set to local in header ([6897fb6](6897fb6f494703e1c8ac9baa56ca776c1623ec13)) + +- Fix bug from last commit so synthetic driver can find its config ([1fe448f](1fe448f789637488ac0e4866c6cef5e7612046b6)) + +- Fix typo from previous commit ([82cfe0e](82cfe0e8944e1a35bf2dcdaaed829776f5ddbdc1)) + +- Fix incorrect indentation ([cfeaa73](cfeaa73db52443c16b72ab3d164581a3733585e3)) + +- Fix tab bug ([3a40892](3a40892527dac7dd0612ce8ca240fed55b64e061)) + +- Fix bug, real part of residuals not being written to z-file ([c5aa885](c5aa88536f22e5d6098030714331d2925fecaf02)) + +- Fix bug issue #161 ([448cf84](448cf84bc7bacaaddd857188a48b6a70474d1924)) + +- Fix bug #161 confirmed, added exception in test if encountered again ([5bf7db0](5bf7db05d39fb399acfb8c2bf09956eb9137dfed)) + +- Fix some typos and get synthetic matlab test comparing against expected values ([228d96d](228d96ddce2f10fb210710687df8ff240341547e)) + +- Fix typo bug ([21c5209](21c5209f9d6ca258ff8753062bcaf8530b417927)) + +- Fix logic handling both styles of processing config to only support new style ([798d544](798d544805a56fad5e9b35d543cc99df6a72fc67)) + +- Fix unclosed else statement bug ([8c716da](8c716dadde5bff9b150a4bca490bba7065857dec)) + +- Fix typo bug ([9c3c601](9c3c6011cef6d4b31a65fe4b8e988b863b218c3f)) + +- Fix bug: pass tfk_dataset, not tfk_dataset.df to create_test_run_config ([792d7ef](792d7ef353b0d655ce716e70c5194e8530f8c0b2)) + +- Fix type ([3fb691e](3fb691e28a38d254f94e65734b1da4a53be538a5)) + +- Fixed issue 90, simplified call to survey_dict at end of porcess_mth5 ([d37763e](d37763ef967c1c1c5f93481c5250caccc8f23906)) + +- Fix bug ([ae9e128](ae9e1287c3b72fdc211cdc284c4126972e304e30)) + +- Fix bug ([e30d6ae](e30d6aedff382951729688311c3a75de54f741fb)) + +- Fix bug (read mode for initialize_mth5s) ([c1e8a31](c1e8a31165e7ea0a3396f3cc736d9344f2aa1c47)) + +- Fix merge conflict ([223ad2a](223ad2a02d6d5cbe375e724e356683b1234150d3)) + +- Fix path bug in io test ([ed4db4b](ed4db4be05883bf4e4dafcd176a47a8e858bc7cc)) + +- Fix path bug in io test ([71d1578](71d1578a8143476ba72f4b2fdae1162396ba71f7)) + +- Fix issue #274 ([201a3f9](201a3f95918eb4659bb7e62ff3667d2f1dadddcf)) + +- Fixed bug in station_metadata from run_group ([c460df0](c460df0a48d816b48aa366574bcb87a4d3db1373)) + +- Fix typo ([6190255](61902556818452d4be33914e04fb877108974f6f)) + +- Fix typo ([851c294](851c294dca7f5657c8eb752038f3956da7b6ff65)) + +- Fix bug ([3a0fe55](3a0fe55271dec39cb508c0b14a9b75a307b2d3e4)) + +- Fix import ([c54be58](c54be5898a5a5d3318f432fe98b5d6b4b9cf2216)) + +- Fix bug + +- build_request_df was error when start/end were not None but time_period_dict unused ([adf4878](adf487841174756b6a6284bc302e0641459173ea)) + +- Fix bug - no data for test ([061e786](061e7866840672c1325af91b0878d13b42dc80f3)) + +- Fix typo ([506ebde](506ebde8e0f8471a86ab30b5744032b137452743)) + + +## Documentation + +- Doc strings ([b4975d4](b4975d4739146481ca01e0f288325089b95010b2)) + +- Documentation tests are failing - not sure why, got a warning about this class, trying remove 'object' from class definition ([be13cec](be13cec7fe10e9143dd01861327d25319b59ce7b)) + +- Docstrings ([c7f2ac9](c7f2ac912a2457e5bf2dfd44551fd9255584edeb)) + +- Docs fails without examples folder ([d7cd250](d7cd250f0a090e6be96756ec5afaa1bffa1651cb)) + +- Doc strings added ([2762ec7](2762ec75b765776ef522036921aa3e5ec3e95fd4)) + +- Docstrings, and shorten check of valid run ([6dd1074](6dd10742674e5043c887ffb92750a5c581b1a4d7)) + +- Docs ([35c2d8a](35c2d8a7340f76b6ba2cccc35419bce3ebce7f10)) + + +## Revert + +- Reverted RMS back to original numbers ([a5e7fea](a5e7fea03d9792488406960632a9a050083673dc)) + +- Revert tests to track fc branch of mth5 ([c5ff9f4](c5ff9f4b854ac8663d5aa4cd37dcd0cb4cd4cdb0)) + +- Revert tests to master branch of mth5 ([697130a](697130ad7bd8cc2bf595cf9e350ffbff46a7ae3b)) + + +## Testing + +- Tested residual_variance_method1 ([c716e12](c716e123fb1a72492ccf1d18c0de66f46b10ad91)) + +- Testing can fix doc test failure (issue #151) by manual edit of .rst ([b5f98f9](b5f98f979630c13139b4f95d22c8809c94566093)) + +- Test_stft_methods_agree using new Processing() class ([357fa87](357fa87f8aba3b5b65d49c0d907b8c686cea21ab)) + +- Test slicing, start/end run times can be used in pipeline ([5f2f625](5f2f625b224d360b7c4a28dcee24e731fa11c968)) + +- Testing to see if error from PR177 reproduces on dev branch ([66cbc1e](66cbc1e7ba4ad4b350625f8d7ffd698d2828b399)) + +- Test change python version in metadata of ipynb ([2d8d0e4](2d8d0e4cdeaea9bd817dac9601ef222b48b3f93e)) + +- Test deprecated library ([3b485a9](3b485a9eb4adf6517a419c6b33f8e94f61066157)) + +- Test all three versions 3.6, 3.7, 3.8 ([b9f535a](b9f535a0fd9e7352b6281615b5aa366cdb742991)) + +- Test precommit ([7a6798d](7a6798dfaea7f78ef90894438fbe28aabc843f3f)) + +- Test precommit ([965b0eb](965b0eba5f908d96f4d2adca626593692697688f)) + +- Test pre-commit ([1c646da](1c646da1240923afa17ecc4a3be4a967bc0840f2)) + +- Test pre-commit ([73943e9](73943e95add0f58afff41fa5b7c66c60734c3557)) + +- Tested can pass very wide time interval to get all data ([0d20049](0d20049431971a2a9fa31250ea73f6e2d316ea19)) + +- Test on all three python versions with mth5 fix_issue_105 branch merged into mth5 main ([48dc886](48dc886fab4279d2b0740b128c1e853602b6d04f)) + +- Test on all versions ([95b4630](95b463038276db9f0d4c01afc48652384123c8cd)) + +- Test on v3.9 ([06a3ab7](06a3ab72557bf6de677d46db317c785f1da84b5a)) + +- Test fix_issue_133 on mt_metadata ([6541f94](6541f943d4460ee20336ccd2cdd10e9b789c373a)) + +- Test build with dask ([9fe13e3](9fe13e30334c68e2eb1f3d783c64976bde80d58f)) + +- Test fix on mt_metadata branch ([e63e9d6](e63e9d6f5148b4bd5416a8204616dc118f9cd083)) + + +## Bug Fixes + +- Fixing mess ([4936516](4936516e01c9c071ce00ec99dc30ccf9d5a4b69a)) + +- Fixing mess ([ad824b5](ad824b556f7dd6cfaafd93f2b1f1e9f76a578a8d)) + +- Fix a path that depended on rename from emtf_synthetic to synthetic ([eb36778](eb367781eed308dfd02205b1ebf22f623ff6bd1c)) + +- Fixed units in pkd processing test ([067c038](067c0389c5708469d882406ea0ee081f571262c2)) + +- Fix bug that RR processing was only occurring for decimation level 1 ([5238c69](5238c69f1fe039c26604fe8ba6fd86a7509f651f)) + +- Fix bug in synthetic data with nan creation ([d635d4a](d635d4a86a0dcf5e499ddf9f36e45404b4e011db)) + +- Fix reversed order periods in adhoc matlab reader ([5bb7c6a](5bb7c6abfd36b5a641fbaa7d9ff6e78d152e4a2b)) + +- Fix a bug with default band setup,add kwarg to config creator ([effc815](effc8152e8da5e65dbe182f9534f36e18e6fe623)) + +- Fix bug in filter stage name -- mt_metadata issue #58 ([56748b4](56748b4e5555f96ce75c42048744b7156c083da1)) + +- Fix introduced bug in calibration test -- new PKD data from NCEDC has 287999 points, not 288000, change window length to come from run_ts_obj.dataset instead of hard coded ([8fb708b](8fb708be72ab32b99c07de42576e06afc8e2ba48)) + +- Fix typos in plot legend and axes labels ([92948b4](92948b49cb6ccce920780f6826f52c30574567ae)) + +- Fix typo ([916318b](916318b9f6d20f39c67b96de5732e2cdf8acee29)) + +- Fix conda install line ([bcf4532](bcf45329d7a9736c8694b6023a289a09b9aaef9a)) + + +## Documentation + +- Document aurora, make readme point to docs ([f948acb](f948acb4475108b0ca71ba5a8f7e774a2f970da6)) + + +## Testing + +- Tests still passing - ready to start swapping out ProcessingConfig for RunConfig objects ([63686bd](63686bd3c12fd906b7089f2f67287b5b8176027c)) + +- Test minor change with black ([9be250c](9be250c04cd68b2e3cde0206aa49d1a5cbacf876)) + +- Test use branch for git actions.see mth5 issue#38 ([7388ba5](7388ba592892e73eada0f54721c2a7aa20ab41fc)) + +- Test hop off fix2 branch ([ebae900](ebae9002d76c6969b341af7c35f7fbd6b5889900)) + +- Test hop off @channel_nomenclature_38 branch ([722c38e](722c38e344a9af07fd20977bbfdf73c7465cd064)) + +- Testing gitactions on main branch of mt_metadata ([e1505aa](e1505aaac55a35a6e9ce1d349c8e974085c1c48d)) + + +## Bug Fixes + +- Fix broken test, prepping to handle detrend options with scipy.signal, change from boolean to linear, constant, False ([1e8468d](1e8468dc1182a239d290374da3fbea40021e8e40)) + diff --git a/aurora/__init__.py b/aurora/__init__.py index 065baaa3..f16c310e 100644 --- a/aurora/__init__.py +++ b/aurora/__init__.py @@ -1,4 +1,4 @@ -__version__ = "0.5.2" +__version__ = "0.6.0" import sys from loguru import logger diff --git a/aurora/config/config_creator.py b/aurora/config/config_creator.py index 76df065a..2076aebc 100644 --- a/aurora/config/config_creator.py +++ b/aurora/config/config_creator.py @@ -5,20 +5,21 @@ Development Notes: - Initially, we only supported EMTF style band specification, but this should be updated - to allow logarithmically spaced bands of constant Q. + to allow logarithmically spaced bands of constant Q. """ +import pathlib from typing import Optional, Union + from loguru import logger +from mt_metadata.processing.window import Window +from mth5.processing.kernel_dataset import KernelDataset from aurora.config import BANDS_DEFAULT_FILE from aurora.config.metadata import Processing from aurora.sandbox.io_helpers.emtf_band_setup import EMTFBandSetupFile -from mth5.processing.kernel_dataset import KernelDataset -from mt_metadata.processing.window import Window -import pathlib SUPPORTED_BAND_SPECIFICATION_STYLES = ["EMTF", "band_edges"] @@ -136,29 +137,22 @@ def create_from_kernel_dataset( **kwargs, ) -> Processing: """ - This creates a processing config from a kernel dataset. + Creates a processing config from a kernel dataset. TODO: Make this a method of kernel_dataset. **Development Notes:** - 1. 2022-09-10 - - The reading-in from EMTF band setup file used to be very terse, carried - some baked in assumptions about decimation factors, and did not acknowledge - specific frequency bands in Hz. I am adding some complexity to the method - that populates bands from EMTF band setup file but am now explict about the - assumption of decimation factors, and do provide the frequency bands in Hz. + 1. The number of decimation levels must be defined either by: - 2. The number of decimation levels must be defined either by: + - decimation_factors argument (normally accompanied by a bands_dict) + - number of decimations implied by EMTF band setup file. - - decimation_factors argument (normally accompanied by a bands_dict) - - number of decimations implied by EMTF band setup file. + Theoretically, you could also use the number of decimations implied by bands_dict but this is sloppy, because it would assume the decimation factor. - Theoretically, you could also use the number of decimations implied by bands_dict but this is sloppy, because it would assume the decimation factor. + 2. 2024-12-29 Added setting of decimation_obj.stft.per_window_detrend_type = "linear" - 3. 2024-12-29 Added setting of decimation_obj.stft.per_window_detrend_type = "linear" - This makes tests pass following a refactoring of mt_metadata. Could use more testing. + This makes tests pass following a refactoring of mt_metadata. Could use more testing. Parameters ---------- @@ -180,12 +174,14 @@ def create_from_kernel_dataset( List of decimation factors, normally [1, 4, 4, 4, ... 4] num_samples_window: Optional[Union[int, None]] The size of the window (usually for FFT) + **kwargs: Additional keyword arguments passed to Processing constructor. Could contain: + - save_fcs: bool - - If True, save Fourier coefficients during processing. + If True, save Fourier coefficients during processing. - save_fcs_type: str - - File type for saving Fourier coefficients. Options are "h5" or "csv". + File type for saving Fourier coefficients. Options are "h5" or "csv". Returns ------- @@ -209,7 +205,6 @@ def create_from_kernel_dataset( # Set Frequency Bands if self.band_specification_style == "EMTF": - # see note 1 emtf_band_setup_file = EMTFBandSetupFile( filepath=self._emtf_band_file, sample_rate=kernel_dataset.sample_rate, diff --git a/aurora/config/metadata/processing.py b/aurora/config/metadata/processing.py index ce59aacf..9550a797 100644 --- a/aurora/config/metadata/processing.py +++ b/aurora/config/metadata/processing.py @@ -8,29 +8,27 @@ # Imports # ============================================================================= -from aurora.time_series.windowing_scheme import window_scheme_from_decimation -from loguru import logger -from mt_metadata.processing.aurora.processing import ( - Processing as AuroraProcessing, -) -from mt_metadata.common.list_dict import ListDict +import json +import pathlib from typing import Optional, Union -import json import pandas as pd -import pathlib +from loguru import logger +from mt_metadata.common.list_dict import ListDict +from mt_metadata.processing.aurora.processing import Processing as AuroraProcessing + +from aurora.time_series.windowing_scheme import window_scheme_from_decimation class Processing(AuroraProcessing): def __init__(self, **kwargs): """ - Constructor + Constructor. Parameters ---------- kwargs """ - # super().__init__(attr_dict=attr_dict, **kwargs) super().__init__(**kwargs) def window_scheme(self, as_type="df"): @@ -192,6 +190,7 @@ class EMTFTFHeader(ListDict): def __init__(self, **kwargs): """ Parameters + ---------- _local_station : mt_metadata.processing.tf.station.Station() Station metadata object for the station to be estimated ( location, channel_azimuths, etc.) diff --git a/aurora/general_helper_functions.py b/aurora/general_helper_functions.py index c176fbfa..b01a62f1 100644 --- a/aurora/general_helper_functions.py +++ b/aurora/general_helper_functions.py @@ -2,18 +2,15 @@ This module contains from miscellaneous functions and some global paths. """ import inspect - -# import os import pathlib -import scipy.io as sio import subprocess +from pathlib import Path +import mt_metadata +import scipy.io as sio from loguru import logger -from pathlib import Path import aurora -import mt_metadata -import mth5 init_file = inspect.getfile(aurora) @@ -61,30 +58,6 @@ def get_test_path() -> pathlib.Path: MT_METADATA_DATA = Path(mt_metadata_init).parent.parent.joinpath("data") -def count_lines(file_name): - """ - acts like wc -l in unix, - raise FileNotFoundError: if file_name does not exist. - - Parameters - ---------- - file_name: str or pathlib.Path - The file to apply line counting to - - Returns - ------- - num_lines: int - Number of lines present in fileName or -1 if file does not exist - - """ - i = -1 - with open(file_name) as f: - for i, l in enumerate(f): - pass - num_lines = i + 1 - return num_lines - - def execute_subprocess(cmd, **kwargs): """ A wrapper for subprocess.call @@ -104,8 +77,8 @@ def execute_subprocess(cmd, **kwargs): def replace_in_file(file_path: pathlib.Path, old: str, new: str) -> None: """ - Replace all instances of 'old' with 'new' in the given file. - :param file_path: Path to the file where replacements should be made. + Replace all instances of 'old' with 'new' in the given file. + :param file_path: Path to the file where replacements should be made. """ if not file_path.exists(): diff --git a/aurora/pipelines/process_mth5.py b/aurora/pipelines/process_mth5.py index ce86af0a..85ee9040 100644 --- a/aurora/pipelines/process_mth5.py +++ b/aurora/pipelines/process_mth5.py @@ -72,7 +72,7 @@ def process_tf_decimation_level( Processing pipeline for a single decimation_level TODO: Add a check that the processing config sample rates agree with the data - sampling rates otherwise raise Exception + sampling rates otherwise raise Exception TODO: Add units to local_stft_obj, remote_stft_obj TODO: This is the method that should be accessing weights This method can be single station or remote based on the process cfg diff --git a/aurora/pipelines/run_summary.py b/aurora/pipelines/run_summary.py index fc833c25..9ff299a9 100644 --- a/aurora/pipelines/run_summary.py +++ b/aurora/pipelines/run_summary.py @@ -1,38 +1,30 @@ """ +Functionality of RunSummary() -Note 1: Functionality of RunSummary() -1. User can get a list of local_station options, which correspond to unique pairs -of values: (survey_id, station_id) - +1. User can get a list of local_station options, which correspond to unique pairs: (survey_id, station_id) 2. User can see all possible ways of processing the data: -- one list per (survey_id, station_id) pair in the run_summary + - one list per (survey_id, station_id) pair in the run_summary Some of the following functionalities may end up in KernelDataset: + 3. User can select local_station --this can trigger a reduction of runs to only those that are from the local staion -and simultaneous runs at other stations + - this can trigger a reduction of runs to only those that are from the local station and simultaneous runs at other stations 4. Given a local station, a list of possible reference stations can be generated -5. Given a remote reference station, a list of all relevent runs, truncated to -maximize coverage of the local station runs is generated +5. Given a remote reference station, a list of all relevant runs, truncated to maximize coverage of the local station runs is generated 6. Given such a "restricted run list", runs can be dropped 7. Time interval endpoints can be changed - - """ import copy -import pandas as pd - +import mth5 +import pandas as pd +from loguru import logger from mt_metadata.transfer_functions import ( ALLOWED_INPUT_CHANNELS, -) -from mt_metadata.transfer_functions import ( ALLOWED_OUTPUT_CHANNELS, ) -import mth5 from mth5.utils.helpers import initialize_mth5 -from loguru import logger RUN_SUMMARY_COLUMNS = [ @@ -51,13 +43,12 @@ class RunSummary: """ - The dependencies aren't clear yet. - Maybe still Dataset: - Could have methods - "drop_runs_shorter_than" - "fill_gaps_by_time_interval" - "fill_gaps_by_run_names" - " + ** Development Notes** + The relationship between this class, MTH5's run_summary method, and the KernelDataset aren't clear yet. + This could have methods for + - "drop_runs_shorter_than" + - "fill_gaps_by_time_interval" + - "fill_gaps_by_run_names" For the full MMT case this may need modification to a channel based summary. @@ -126,9 +117,7 @@ def check_runs_are_valid(self, drop=False, **kwargs): run_obj = m.get_run(row.station_id, row.run_id, row.survey) runts = run_obj.to_runts() if runts.dataset.to_array().data.__abs__().sum() == 0: - logger.critical( - "CRITICAL: Detected a run with all zero values" - ) + logger.critical("CRITICAL: Detected a run with all zero values") self.df["valid"].at[i_row] = False # load each run, and take the median of the sum of the absolute values if drop: @@ -159,57 +148,33 @@ def channel_summary_to_run_summary( allowed_input_channels=ALLOWED_INPUT_CHANNELS, allowed_output_channels=ALLOWED_OUTPUT_CHANNELS, sortby=["station_id", "start"], -): +) -> pd.DataFrame: """ + Create a run summary dataframe from a channel summary dataframe. The run summary has one row per acquisition run, and includes columns for station_id, run_id, start, end, sample_rate, input_channels, output_channels, and channel_scale_factors. + + ** Development Notes** TODO: replace station_id with station, and run_id with run Note will need to modify: aurora/tests/config$ more test_dataset_dataframe.py - TODO: Add logic for handling input and output channels based on channel - summary. Specifically, consider the case where there is no vertical magnetic - field, this information is available via ch_summary, and output channels should - then not include hz. - TODO: Just inherit all the run-level and higher el'ts of the channel_summary, - including n_samples? - - When creating the dataset dataframe, make it have these columns: - [ - "station_id", - "run_id", - "start", - "end", - "mth5_path", - "sample_rate", - "input_channels", - "output_channels", - "remote", - "channel_scale_factors", - ] + TODO: Add logic for handling input and output channels based on channel summary. Specifically, consider the case where there is no vertical magnetic field, this information is available via ch_summary, and output channels should then not include hz. + TODO: Just inherit all the run-level and higher el'ts of the channel_summary, including n_samples? Parameters ---------- ch_summary: mth5.tables.channel_table.ChannelSummaryTable or pandas DataFrame - If its a dataframe it is a representation of an mth5 channel_summary. - Maybe restricted to only have certain stations and runs before being passed to - this method + If it's a dataframe it is a representation of an mth5 channel_summary. Maybe restricted to only have certain stations and runs before being passed to this method. allowed_input_channels: list of strings Normally ["hx", "hy", ] - These are the allowable input channel names for the processing. See further - note under allowed_output_channels. + These are the allowable input channel names for the processing. See further note under allowed_output_channels. allowed_output_channels: list of strings Normally ["ex", "ey", "hz", ] - These are the allowable output channel names for the processing. - A global list of these is kept at the top of this module. The purpose of - this is to distinguish between runs that have different layouts, for example - some runs will have hz and some will not, and we cannot process for hz the - runs that do not have it. By making this a kwarg we sort of prop the door - open for more general names (see issue #74). + These are the allowable output channel names for the processing. A global list of these is kept at the top of this module. The purpose of this is to distinguish between runs that have different layouts, for example some runs will have hz and some will not, and we cannot process for hz the runs that do not have it. By making this a kwarg we sort of prop the door open for more general names (see issue #74). sortby: bool or list Default: ["station_id", "start"] Returns ------- - run_summary_df: pd.Dataframe - A table with one row per "acquistion run" that was in the input channel - summary table + run_summary_df: pd.DataFrame + A table with one row per "acquisition run" that was in the input channel summary table. """ if isinstance(ch_summary, mth5.tables.channel_table.ChannelSummaryTable): ch_summary_df = ch_summary.to_dataframe() @@ -229,9 +194,7 @@ def channel_summary_to_run_summary( channel_scale_factors = n_station_runs * [None] i = 0 for group_values, group in grouper: - group_info = dict( - zip(group_by_columns, group_values) - ) # handy for debug + group_info = dict(zip(group_by_columns, group_values)) # handy for debug # for k, v in group_info.items(): # print(f"{k} = {v}") survey_ids[i] = group_info["survey"] @@ -242,15 +205,9 @@ def channel_summary_to_run_summary( sample_rates[i] = group.sample_rate.iloc[0] channels_list = group.component.to_list() num_channels = len(channels_list) - input_channels[i] = [ - x for x in channels_list if x in allowed_input_channels - ] - output_channels[i] = [ - x for x in channels_list if x in allowed_output_channels - ] - channel_scale_factors[i] = dict( - zip(channels_list, num_channels * [1.0]) - ) + input_channels[i] = [x for x in channels_list if x in allowed_input_channels] + output_channels[i] = [x for x in channels_list if x in allowed_output_channels] + channel_scale_factors[i] = dict(zip(channels_list, num_channels * [1.0])) i += 1 data_dict = {} @@ -302,9 +259,7 @@ def extract_run_summary_from_mth5(mth5_obj, summary_type="run"): return out_df -def extract_run_summaries_from_mth5s( - mth5_list, summary_type="run", deduplicate=True -): +def extract_run_summaries_from_mth5s(mth5_list, summary_type="run", deduplicate=True): """ ToDo: Move this method into mth5? or mth5_helpers? ToDo: Make this a class so that the __repr__ is a nice visual representation of the diff --git a/aurora/sandbox/io_helpers/emtf_band_setup.py b/aurora/sandbox/io_helpers/emtf_band_setup.py index bf060a0f..3b0d8ff4 100644 --- a/aurora/sandbox/io_helpers/emtf_band_setup.py +++ b/aurora/sandbox/io_helpers/emtf_band_setup.py @@ -2,11 +2,12 @@ In this module is a class that emulates the old EMTF Band Setup File """ -from loguru import logger -from typing import Union, Optional +import pathlib +from typing import Optional, Union + import numpy as np import pandas as pd -import pathlib +from loguru import logger class EMTFBandSetupFile: @@ -72,7 +73,7 @@ def load(self, filepath: Optional[Union[str, pathlib.Path, None]] = None) -> Non df = pd.read_csv( filepath, skiprows=1, - sep="\s+", + sep=r"\s+", names=["decimation_level", "lower_bound_index", "upper_bound_index"], ) if len(df) != self.num_bands: diff --git a/aurora/sandbox/io_helpers/make_mth5_helpers.py b/aurora/sandbox/io_helpers/make_mth5_helpers.py index 3efee6c8..aa3929e4 100644 --- a/aurora/sandbox/io_helpers/make_mth5_helpers.py +++ b/aurora/sandbox/io_helpers/make_mth5_helpers.py @@ -52,15 +52,6 @@ def create_from_server_multistation( the electric or magnetic field channels. triage_missing_coil: bool - - Returns - ------- - - target_folder: - run_id : string - This is a temporary workaround. A more robust program that assigns run - numbers, and/or gets run labels from StationXML is needed - Returns ------- h5_path: pathlib.Path diff --git a/aurora/sandbox/io_helpers/zfile_murphy.py b/aurora/sandbox/io_helpers/zfile_murphy.py deleted file mode 100644 index 8d397e9b..00000000 --- a/aurora/sandbox/io_helpers/zfile_murphy.py +++ /dev/null @@ -1,440 +0,0 @@ -""" -This module contains a class that was contributed by Ben Murphy for working with EMTF "Z-files" -""" - -import pathlib -import re -from typing import Optional, Union - -import numpy as np - - -class ZFile: - def __init__(self, filename: Union[str, pathlib.Path]): - """ - Constructor - - Parameters - ---------- - filename: Union[str, pathlib.Path] - The path to the z-file. - - """ - self.filename = filename - self.station = "" - self.decimation_levels = None - self.lower_harmonic_indices = None - self.upper_harmonic_indices = None - self.f = None - - def open_file(self) -> None: - """attempt to open file""" - try: - self.f = open(self.filename, "r") - except IOError: - raise IOError("File not found.") - return - - def skip_header_lines(self) -> None: - """Skip over the header when reading""" - f = self.f - f.readline() - f.readline() - f.readline() - return - - def get_station_id(self) -> None: - """get station ID from zfile""" - f = self.f - line = f.readline() - if line.lower().startswith("station"): - station = line.strip().split(":", 1)[1] - else: - station = line.strip() - self.station = station - - def read_coordinates_and_declination(self) -> None: - """read coordinates and declination""" - f = self.f - line = f.readline().strip().lower() - match = re.match( - r"\s*coordinate\s+(-?\d+\.?\d*)\s+(-?\d+\.?\d*)\s+" - r"declination\s+(-?\d+\.?\d*)", - line, - ) - self.coordinates = (float(match.group(1)), float(match.group(2))) - self.declination = float(match.group(3)) - return - - def read_number_of_channels_and_number_of_frequencies(self): - """read_number_of_channels_and_number_of_frequencies""" - f = self.f - line = f.readline().strip().lower() - match = re.match( - r"\s*number\s+of\s+channels\s+(\d+)\s+number\s+of" - r"\s+frequencies\s+(\d+)", - line, - ) - nchannels = int(match.group(1)) - nfreqs = int(match.group(2)) - self.nchannels = nchannels - self.nfreqs = nfreqs - - def read_channel_information(self): - """read_channel_information""" - f = self.f - self.orientation = np.zeros((self.nchannels, 2)) - self.channels = [] - for i in range(self.nchannels): - line = f.readline().strip() - match = re.match( - r"\s*\d+\s+(-?\d+\.?\d*)\s+(-?\d+\.?\d*)\s+\w*\s+" r"(\w+)", line - ) - self.orientation[i, 0] = float(match.group(1)) - self.orientation[i, 1] = float(match.group(2)) - if len(match.group(3)) > 2: - # sometimes the channel ID comes out with extra stuff - self.channels.append(match.group(3)[:2].title()) - else: - self.channels.append(match.group(3).title()) - return - - def load(self): - """load Z-file and populate attributes of class""" - self.open_file() - self.skip_header_lines() - self.get_station_id() - self.read_coordinates_and_declination() - self.read_number_of_channels_and_number_of_frequencies() - - f = self.f - - # skip line - f.readline() - self.read_channel_information() - f.readline() - - # initialize empty arrays for transfer functions - # note that EMTF, and consequently this code, assumes two predictor - # channels (horizontal magnetics) - # nchannels - 2 therefore is the number of predicted channels - self.decimation_levels = np.zeros(self.nfreqs) - self.periods = np.zeros(self.nfreqs) - self.lower_harmonic_indices = np.zeros(self.nfreqs) - self.upper_harmonic_indices = np.zeros(self.nfreqs) - self.transfer_functions = np.zeros( - (self.nfreqs, self.nchannels - 2, 2), dtype=np.complex64 - ) - - # residual covariance -- square matrix with dimension as number of - # predicted channels - self.sigma_e = np.zeros( - (self.nfreqs, self.nchannels - 2, self.nchannels - 2), dtype=np.complex64 - ) - - # inverse coherent signal power -- square matrix, with dimension as the - # number of predictor channels - # since EMTF and this code assume N predictors is 2, - # this dimension is hard-coded - self.sigma_s = np.zeros((self.nfreqs, 2, 2), dtype=np.complex64) - - # now read data for each period - for i in range(self.nfreqs): - # extract period - line = f.readline().strip() - match = re.match( - r"\s*period\s*:\s+(\d+\.?\d*)\s+" r"decimation\s+level", line - ) - self.periods[i] = float(match.group(1)) - - splitted_line1 = line.split("level") - splitted_line2 = splitted_line1[1].split("freq") - self.decimation_levels[i] = int(splitted_line2[0].strip()) - splitted_line1 = line.split("from") - splitted_line2 = splitted_line1[1].split("to") - self.lower_harmonic_indices[i] = int(splitted_line2[0].strip()) - self.upper_harmonic_indices[i] = int(splitted_line2[1].strip()) - # skip two lines - f.readline() - f.readline() - - # read transfer functions - for j in range(self.nchannels - 2): - comp1_r, comp1_i, comp2_r, comp2_i = f.readline().strip().split() - self.transfer_functions[i, j, 0] = float(comp1_r) + 1.0j * float( - comp1_i - ) - self.transfer_functions[i, j, 1] = float(comp2_r) + 1.0j * float( - comp2_i - ) - - # skip label line - f.readline() - - # read inverse coherent signal power matrix - val1_r, val1_i = f.readline().strip().split() - val2_r, val2_i, val3_r, val3_i = f.readline().strip().split() - self.sigma_s[i, 0, 0] = float(val1_r) + 1.0j * float(val1_i) - self.sigma_s[i, 1, 0] = float(val2_r) + 1.0j * float(val2_i) - self.sigma_s[i, 0, 1] = float(val2_r) - 1.0j * float(val2_i) - self.sigma_s[i, 1, 1] = float(val3_r) + 1.0j * float(val3_i) - - # skip label line - f.readline() - - # read residual covariance - for j in range(self.nchannels - 2): - values = f.readline().strip().split() - for k in range(j + 1): - if j == k: - self.sigma_e[i, j, k] = float(values[2 * k]) + 1.0j * float( - values[2 * k + 1] - ) - else: - self.sigma_e[i, j, k] = float(values[2 * k]) + 1.0j * float( - values[2 * k + 1] - ) - self.sigma_e[i, k, j] = float(values[2 * k]) - 1.0j * float( - values[2 * k + 1] - ) - - f.close() - - def impedance(self, angle: Optional[float] = 0.0): - """ - Compute the Impedance and errors from the transfer function. - - note u,v are identity matrices if angle=0 - - Parameters - ---------- - angle: float - The angle about the vertical axis by which to rotate the Z tensor. - - Returns - ------- - z: np.ndarray - The impedance tensor - error: np.ndarray - Errors for the impedance tensor - """ - # check to see if there are actually electric fields in the TFs - if "Ex" not in self.channels and "Ey" not in self.channels: - raise ValueError( - "Cannot return apparent resistivity and phase " - "data because these TFs do not contain electric " - "fields as a predicted channel." - ) - - # transform the TFs first... - # build transformation matrix for predictor channels - # (horizontal magnetic fields) - hx_index = self.channels.index("Hx") - hy_index = self.channels.index("Hy") - u = np.eye(2, 2) - u[hx_index, hx_index] = np.cos( - (self.orientation[hx_index, 0] - angle) * np.pi / 180.0 - ) - u[hx_index, hy_index] = np.sin( - (self.orientation[hx_index, 0] - angle) * np.pi / 180.0 - ) - u[hy_index, hx_index] = np.cos( - (self.orientation[hy_index, 0] - angle) * np.pi / 180.0 - ) - u[hy_index, hy_index] = np.sin( - (self.orientation[hy_index, 0] - angle) * np.pi / 180.0 - ) - u = np.linalg.inv(u) # Identity if angle=0 - - # build transformation matrix for predicted channels (electric fields) - ex_index = self.channels.index("Ex") - ey_index = self.channels.index("Ey") - v = np.eye(self.transfer_functions.shape[1], self.transfer_functions.shape[1]) - v[ex_index - 2, ex_index - 2] = np.cos( - (self.orientation[ex_index, 0] - angle) * np.pi / 180.0 - ) - v[ey_index - 2, ex_index - 2] = np.sin( - (self.orientation[ex_index, 0] - angle) * np.pi / 180.0 - ) - v[ex_index - 2, ey_index - 2] = np.cos( - (self.orientation[ey_index, 0] - angle) * np.pi / 180.0 - ) - v[ey_index - 2, ey_index - 2] = np.sin( - (self.orientation[ey_index, 0] - angle) * np.pi / 180.0 - ) - - # matrix multiplication... - rotated_transfer_functions = np.matmul( - v, np.matmul(self.transfer_functions, u.T) - ) - rotated_sigma_s = np.matmul(u, np.matmul(self.sigma_s, u.T)) - rotated_sigma_e = np.matmul(v, np.matmul(self.sigma_e, v.T)) - - # now pull out the impedance tensor - z = np.zeros((self.periods.size, 2, 2), dtype=np.complex64) - z[:, 0, 0] = rotated_transfer_functions[:, ex_index - 2, hx_index] # Zxx - z[:, 0, 1] = rotated_transfer_functions[:, ex_index - 2, hy_index] # Zxy - z[:, 1, 0] = rotated_transfer_functions[:, ey_index - 2, hx_index] # Zyx - z[:, 1, 1] = rotated_transfer_functions[:, ey_index - 2, hy_index] # Zyy - - # and the variance information - var = np.zeros((self.periods.size, 2, 2)) - var[:, 0, 0] = np.real( - rotated_sigma_e[:, ex_index - 2, ex_index - 2] - * rotated_sigma_s[:, hx_index, hx_index] - ) - - var[:, 0, 1] = np.real( - rotated_sigma_e[:, ex_index - 2, ex_index - 2] - * rotated_sigma_s[:, hy_index, hy_index] - ) - var[:, 1, 0] = np.real( - rotated_sigma_e[:, ey_index - 2, ey_index - 2] - * rotated_sigma_s[:, hx_index, hx_index] - ) - var[:, 1, 1] = np.real( - rotated_sigma_e[:, ey_index - 2, ey_index - 2] - * rotated_sigma_s[:, hy_index, hy_index] - ) - - error = np.sqrt(var) - - return z, error - - def tippers(self, angle=0.0): - """compute the tipper from transfer function""" - - # check to see if there is a vertical magnetic field in the TFs - if "Hz" not in self.channels: - raise ValueError( - "Cannot return tipper data because the TFs do not " - "contain the vertical magnetic field as a " - "predicted channel." - ) - - # transform the TFs first... - # build transformation matrix for predictor channels - # (horizontal magnetic fields) - hx_index = self.channels.index("Hx") - hy_index = self.channels.index("Hy") - u = np.eye(2, 2) - u[hx_index, hx_index] = np.cos( - (self.orientation[hx_index, 0] - angle) * np.pi / 180.0 - ) - u[hx_index, hy_index] = np.sin( - (self.orientation[hx_index, 0] - angle) * np.pi / 180.0 - ) - u[hy_index, hx_index] = np.cos( - (self.orientation[hy_index, 0] - angle) * np.pi / 180.0 - ) - u[hy_index, hy_index] = np.sin( - (self.orientation[hy_index, 0] - angle) * np.pi / 180.0 - ) - u = np.linalg.inv(u) - - # don't need to transform predicated channels (assuming no tilt in Hz) - hz_index = self.channels.index("Hz") - v = np.eye(self.transfer_functions.shape[1], self.transfer_functions.shape[1]) - - # matrix multiplication... - rotated_transfer_functions = np.matmul( - v, np.matmul(self.transfer_functions, u.T) - ) - rotated_sigma_s = np.matmul(u, np.matmul(self.sigma_s, u.T)) - rotated_sigma_e = np.matmul(v, np.matmul(self.sigma_e, v.T)) - - # now pull out tipper information - tipper = np.zeros((self.periods.size, 2), dtype=np.complex64) - tipper[:, 0] = rotated_transfer_functions[:, hz_index - 2, hx_index] # Tx - tipper[:, 1] = rotated_transfer_functions[:, hz_index - 2, hy_index] # Ty - - # and the variance/error information - var = np.zeros((self.periods.size, 2)) - var[:, 0] = np.real( - rotated_sigma_e[:, hz_index - 2, hz_index - 2] - * rotated_sigma_s[:, hx_index, hx_index] - ) # Tx - var[:, 1] = np.real( - rotated_sigma_e[:, hz_index - 2, hz_index - 2] - * rotated_sigma_s[:, hy_index, hy_index] - ) # Ty - error = np.sqrt(var) - - return tipper, error - - def apparent_resistivity(self, angle: float = 0.0): - """compute the apparent resistivity from the impedance.""" - z_tensor, error = self.impedance(angle=angle) - Zxy = z_tensor[:, 0, 1] - Zyx = z_tensor[:, 1, 0] - T = self.periods - self.rxy = T * (abs(Zxy) ** 2) / 5.0 - self.ryx = T * (abs(Zyx) ** 2) / 5.0 - self.pxy = np.rad2deg(np.arctan(np.imag(Zxy) / np.real(Zxy))) - self.pyx = np.rad2deg(np.arctan(np.imag(Zyx) / np.real(Zyx))) - return - - def rho(self, mode): - """ - Return the apparent resistivity for the given mode. - - Convenience function to help with streamlining synthetic tests - to be - eventually replaced by functionality in mt_metadata.tf - - Parameters - ---------- - mode: str - "xy" or "yx" - - Returns - ------- - rho - """ - if mode == "xy": - return self.rxy - if mode == "yx": - return self.ryx - - def phi(self, mode): - """ - Return the phase for the given mode. - - Convenience function to help with streamlining synthetic tests - to be - eventually replaced by functionality in mt_metadata.tf - Parameters - ---------- - mode: str - "xy" or "yx" - - Returns - ------- - phi - """ - if mode == "xy": - return self.pxy - if mode == "yx": - return self.pyx - - -def read_z_file(z_file_path, angle=0.0) -> ZFile: - """ - Reads a zFile and returns a ZFile object. - - Parameters - ---------- - z_file_path: string or pathlib.Path - The name of the EMTF-style z-file to operate on - angle: float - How much rotation to apply. This is a kludge variable used to help compare - legacy SPUD results which are rotated onto a cardinal grid, vs aurora which - store the TF in the coordinate system of acquisition - - Returns - ------- - z_obj: ZFile - The zFile as an object. - - """ - z_obj = ZFile(z_file_path) - z_obj.load() - z_obj.apparent_resistivity(angle=angle) - return z_obj diff --git a/aurora/test_utils/synthetic/plot_helpers_synthetic.py b/aurora/test_utils/synthetic/plot_helpers_synthetic.py index e8a9a44a..c12ccaf8 100644 --- a/aurora/test_utils/synthetic/plot_helpers_synthetic.py +++ b/aurora/test_utils/synthetic/plot_helpers_synthetic.py @@ -8,6 +8,9 @@ def make_subtitle( ): """ + Development Notes: + compute_rms is now deprecated. If that tool is desired it should be reimplemented + in aurora/transfer_function/compare.py Parameters ---------- rho_rms_aurora: float diff --git a/aurora/test_utils/synthetic/rms_helpers.py b/aurora/test_utils/synthetic/rms_helpers.py deleted file mode 100644 index 7a8f26b8..00000000 --- a/aurora/test_utils/synthetic/rms_helpers.py +++ /dev/null @@ -1,125 +0,0 @@ -""" - This module contains methods associated with RMS calculations that are used in testing - aurora processing on synthetic data. - -""" -import numpy as np -from loguru import logger - - -def compute_rms(rho, phi, model_rho_a=100.0, model_phi=45.0, verbose=False): - """ - Computes the RMS between processing results (rho, phi) and model (rho, phi). - - It is used to make annotations for comparative plots for synthetic data. Could be - used in general to compare different processing results. For example by replacing - model_rho_a and model_phi with other processing results, or other (non-uniform) - model results. - - Parameters - ---------- - rho: numpy.ndarray - 1D array of computed apparent resistivities (expected in Ohmm) - phi: numpy.ndarrayx - 1D array of computed phases (expected in degrees) - model_rho_a: float or numpy array - if numpy array must be the same shape as rho - model_phi: float or numpy array - if numpy array must be the same shape as phi. - Returns - ------- - rho_rms: float - rms misfit between the model apparent resistivity and the computed resistivity - phi_rms: float - rms misfit between the model phase (or phases) and the computed phase - """ - rho_rms = np.sqrt(np.mean((rho - model_rho_a) ** 2)) - phi_rms = np.sqrt(np.mean((phi - model_phi) ** 2)) - if verbose: - logger.info(f"rho_rms = {rho_rms}") - logger.info(f"phi_rms = {phi_rms}") - return rho_rms, phi_rms - - -def get_expected_rms_misfit(test_case_id: str, emtf_version=None) -> dict: - """ - Returns hard-coded expected results from synthetic data processing. - These results are a benchmark against which test results are compared on push to - github. - - Parameters - ---------- - test_case_id - emtf_version - - Returns - ------- - - """ - expected_rms_misfit = {} - expected_rms_misfit["rho"] = {} - expected_rms_misfit["phi"] = {} - if test_case_id == "test1": - if emtf_version == "fortran": - expected_rms_misfit["rho"]["xy"] = 4.433905 - expected_rms_misfit["phi"]["xy"] = 0.910484 - expected_rms_misfit["rho"]["yx"] = 3.658614 - expected_rms_misfit["phi"]["yx"] = 0.844645 - elif emtf_version == "matlab": - expected_rms_misfit["rho"]["xy"] = 2.706098 - expected_rms_misfit["phi"]["xy"] = 0.784229 - expected_rms_misfit["rho"]["yx"] = 3.745280 - expected_rms_misfit["phi"]["yx"] = 1.374938 - - elif test_case_id == "test2r1": - expected_rms_misfit["rho"]["xy"] = 3.971313 - expected_rms_misfit["phi"]["xy"] = 0.982613 - expected_rms_misfit["rho"]["yx"] = 3.967259 - expected_rms_misfit["phi"]["yx"] = 1.62881 - return expected_rms_misfit - - -def assert_rms_misfit_ok( - expected_rms_misfit, - xy_or_yx, - rho_rms_aurora, - phi_rms_aurora, - rho_tol=1e-4, - phi_tol=1e-4, -) -> None: - """ - Compares actual RMS misfit from processing against expected values. - Raises Assertion errors if test processing results different from expected. - - Parameters - ---------- - expected_rms_misfit: dictionary - precomputed RMS misfits for test data in rho and phi - xy_or_yx: str - mode - rho_rms_aurora: float - phi_rms_aurora: float - """ - expected_rms_rho = expected_rms_misfit["rho"][xy_or_yx] - expected_rms_phi = expected_rms_misfit["phi"][xy_or_yx] - logger.info(f"expected_rms_rho_{xy_or_yx} {expected_rms_rho}") - logger.info(f"expected_rms_phi_{xy_or_yx} {expected_rms_phi}") - if not np.isclose(rho_rms_aurora - expected_rms_rho, 0, atol=rho_tol): - logger.error("==== AURORA ====\n") - logger.error(rho_rms_aurora) - logger.error("==== EXPECTED ====\n") - logger.error(expected_rms_rho) - logger.error("==== DIFFERENCE ====\n") - logger.error(rho_rms_aurora - expected_rms_rho) - raise AssertionError("Expected misfit for resistivity is not correct") - - if not np.isclose(phi_rms_aurora - expected_rms_phi, 0, atol=phi_tol): - logger.error("==== AURORA ====\n") - logger.error(phi_rms_aurora) - logger.error("==== EXPECTED ====\n") - logger.error(expected_rms_phi) - logger.error("==== DIFFERENCE ====\n") - logger.error(phi_rms_aurora - expected_rms_phi) - raise AssertionError("Expected misfit for phase is not correct") - - return diff --git a/aurora/time_series/apodization_window.py b/aurora/time_series/apodization_window.py index 7789bda5..100a5183 100644 --- a/aurora/time_series/apodization_window.py +++ b/aurora/time_series/apodization_window.py @@ -50,10 +50,10 @@ """ +from typing import Optional, Union + import numpy as np import scipy.signal as ssig -from loguru import logger -from typing import Optional, Union class ApodizationWindow(object): @@ -90,7 +90,6 @@ def __init__( taper_additional_args: Optional[Union[dict, None]] = None, **kwargs, # needed as a passthrough argument to WindowingScheme ): - """ Constructor @@ -137,7 +136,6 @@ def summary(self) -> str: out_str: str String comprised of the taper_family, number_of_samples, and True/False if self.taper is not None """ - self.test_linear_spectral_density_factor() string1 = f"{self.taper_family} {self.num_samples_window}" string1 += f" taper_exists = {bool(self.taper.any())}" string2 = f"NENBW = {self.nenbw:.3f}, CG = {self.coherent_gain:.3f}, " @@ -231,37 +229,6 @@ def enbw(self, fs) -> float: """ return fs * self.S2 / (self.S1**2) - def test_linear_spectral_density_factor(self) -> None: - """T - This is just a test to verify some algebra - - TODO Move this into tests - - Claim: - The lsd_calibration factors - A (1./coherent\_gain)\*np.sqrt((2\*dt)/(nenbw\*N)) - and - B np.sqrt(2/(sample\_rate\*self.S2)) - - are identical. - - Note sqrt(2\*dt)==sqrt(2\*sample_rate) so we can cancel these terms and - A=B IFF - - (1./coherent_gain) \* np.sqrt(1/(nenbw\*N)) == 1/np.sqrt(S2) - which I show in githib aurora issue #3 via . - (CG\*\*2) \* NENBW \*N = S2 - - """ - lsd_factor1 = (1.0 / self.coherent_gain) * np.sqrt( - 1.0 / (self.nenbw * self.num_samples_window) - ) - lsd_factor2 = 1.0 / np.sqrt(self.S2) - if not np.isclose(lsd_factor1, lsd_factor2): - logger.error(f"factor1 {lsd_factor1} vs factor2 {lsd_factor2}") - logger.error("Incompatible spectral density factors") - raise Exception - @property def taper(self) -> np.ndarray: """Returns the values of the taper as a numpy array""" diff --git a/aurora/time_series/decorators.py b/aurora/time_series/decorators.py index 6f61c7df..e908f822 100644 --- a/aurora/time_series/decorators.py +++ b/aurora/time_series/decorators.py @@ -1,22 +1,24 @@ -""" +r""" This module is a place to put decorators used by methods in aurora. It is a work in progress. Development notes: -| Here is the decorator pattern -| def decorator(func): -| @functools.wraps(func) -| def wrapper_decorator(\*args, \*\*kwargs): -| # Do something before -| value = func\*args, \*\*kwargs) -| # Do something after -| return value -| return wrapper_decorator +Here is the decorator pattern:: + + def decorator(func): + @functools.wraps(func) + def wrapper_decorator(*args, **kwargs): + # Do something before + value = func(*args, **kwargs) + # Do something after + return value + return wrapper_decorator """ import functools + import xarray as xr @@ -43,7 +45,6 @@ def can_use_xr_dataarray(func): @functools.wraps(func) def wrapper_decorator(*args, **kwargs): - if isinstance(kwargs["data"], xr.DataArray): kwargs["data"] = kwargs["data"].to_dataset("channel") input_was_dataarray = True diff --git a/aurora/transfer_function/TTFZ.py b/aurora/transfer_function/TTFZ.py index 835a1322..4b9097fb 100644 --- a/aurora/transfer_function/TTFZ.py +++ b/aurora/transfer_function/TTFZ.py @@ -458,23 +458,13 @@ def set_rho_limits(self): def set_lims(self) -> list: """ - Set limits for the plotting axes - - TODO: Add doc or start using MTpy - - Matlab Notes: - set default limits for plotting; QD, derived from ZPLT use max/min limits of periods, rho to set limits - - function[lims, orient] = set_lims(obj) - Returns - lims : list - x_max, x_min, y_max, y_min, 0, 90 - orient: 0 + Set limits for the plotting axes - to be depreacted for MTpy plotting functions. Returns ------- lims: list The plotting limits for period, rho and phi. + [x_max, x_min, y_max, y_min, 0, 90] """ period_min, period_max = self.set_period_limits() # get limits for the x-axis rho_min, rho_max = self.set_rho_limits() @@ -486,8 +476,7 @@ def set_lims(self) -> list: rho_max = 1e4 lims = [period_min, period_max, rho_min, rho_max, phi_min, phi_max] - # orient = 0.0 - return lims # , orient + return lims def plot(self, station_id="Transfer Function", out_filename=None, **kwargs): """ diff --git a/aurora/transfer_function/compare.py b/aurora/transfer_function/compare.py index 5e6cb814..2c458616 100644 --- a/aurora/transfer_function/compare.py +++ b/aurora/transfer_function/compare.py @@ -77,6 +77,11 @@ def plot_two_transfer_functions( label_01="emtf", label_02="aurora", save_plot_path=None, + rho_xx_ylims=None, + rho_xy_ylims=None, + rho_yx_ylims=None, + rho_yy_ylims=None, + phi_ylims=None, ): """ Plots two transfer functions for comparison. @@ -94,10 +99,12 @@ def plot_two_transfer_functions( ------- """ + xy = "xy" fig = plt.figure(figsize=(12, 6)) for ii in range(2): for jj in range(2): + component = f"{xy[ii]}{xy[jj]}" plot_num_res = 1 + ii * 2 + jj plot_num_phase = 5 + ii * 2 + jj ax = fig.add_subplot(2, 4, plot_num_res) @@ -124,9 +131,17 @@ def plot_two_transfer_functions( ax.set_title(self._comp_dict[plot_num_res]) # ax.set_xlabel("Period (s)") if plot_num_res == 1: - ax.set_ylabel("Apparent Resistivity ($\Omega \cdot m$)") + ax.set_ylabel("Apparent Resistivity ($\\Omega \\cdot m$)") ax.legend() ax.grid(True, which="both", ls="--", lw=0.5, color="gray") + if component == "xx" and rho_xx_ylims is not None: + ax.set_ylim(rho_xx_ylims) + if component == "xy" and rho_xy_ylims is not None: + ax.set_ylim(rho_xy_ylims) + if component == "yx" and rho_yx_ylims is not None: + ax.set_ylim(rho_yx_ylims) + if component == "yy" and rho_yy_ylims is not None: + ax.set_ylim(rho_yy_ylims) ax2 = fig.add_subplot(2, 4, plot_num_phase) ax2.semilogx( @@ -149,6 +164,7 @@ def plot_two_transfer_functions( if plot_num_phase == 5: ax2.set_ylabel("Phase (degrees)") ax2.legend() + ax2.grid(True, which="both", ls="--", lw=0.5, color="gray") fig.tight_layout() @@ -286,6 +302,7 @@ def compare_transfer_functions( self, rtol: float = 1, atol: float = 1, + atol_phase: float = 4.0, ) -> dict: """ Compare transfer functions between two transfer_functions objects. @@ -299,6 +316,8 @@ def compare_transfer_functions( Relative tolerance for np.allclose, defaults to 1e-2 atol: float Absolute tolerance for np.allclose, defaults to 1e-2 + atol_phase: float + Absolute tolerance for phase comparison, defaults to 4.0 degrees Returns ------- @@ -357,7 +376,7 @@ def compare_transfer_functions( np.angle(z1[:, ii, jj]), np.angle(z2[:, ii, jj]), rtol=rtol, - atol=atol, + atol=atol_phase, ) result["impedance_error_close"] = np.allclose( diff --git a/aurora/transfer_function/emtf_z_file_helpers.py b/aurora/transfer_function/emtf_z_file_helpers.py index a4fca770..031abdde 100644 --- a/aurora/transfer_function/emtf_z_file_helpers.py +++ b/aurora/transfer_function/emtf_z_file_helpers.py @@ -8,14 +8,10 @@ """ import pathlib +from typing import Optional, Union -import numpy as np -from aurora.transfer_function.transfer_function_collection import ( - TransferFunctionCollection, -) -from aurora.sandbox.io_helpers.zfile_murphy import ZFile from loguru import logger -from typing import Optional, Union + EMTF_CHANNEL_ORDER = ["hx", "hy", "hz", "ex", "ey"] @@ -47,65 +43,6 @@ def get_default_orientation_block(n_ch: int = 5) -> list: return orientation_strs -def merge_tf_collection_to_match_z_file( - aux_data: ZFile, tf_collection: TransferFunctionCollection -) -> dict: - """ - method to merge tf data from a tf_collection with a Z-file when there are potentially - multiple estimates of TF at the same periods for different decimation levels. - - Development Notes: - Currently this is only used for the the synthetic test where aurora results - are compared against a stored legacy Z-file. Given data from a z_file, and a - tf_collection, the tf_collection may have several TF estimates at the same - frequency from multiple decimation levels. This tries to make a single array as - a function of period for all rho and phi. - - Parameters - ---------- - aux_data: aurora.sandbox.io_helpers.zfile_murphy.ZFile - Object representing a z-file - tf_collection: aurora.transfer_function.transfer_function_collection - .TransferFunctionCollection - Object representing the transfer function returned from the aurora processing - - - Returns - ------- - result: dict of dicts - Keyed by ["rho", "phi"], below each of these is an ["xy", "yx",] entry. The - lowest level entries are numpy arrays. - """ - rxy = np.full(len(aux_data.decimation_levels), np.nan) - ryx = np.full(len(aux_data.decimation_levels), np.nan) - pxy = np.full(len(aux_data.decimation_levels), np.nan) - pyx = np.full(len(aux_data.decimation_levels), np.nan) - dec_levels = list(set(aux_data.decimation_levels)) - dec_levels = [int(x) for x in dec_levels] - dec_levels.sort() - - for dec_level in dec_levels: - aurora_tf = tf_collection.tf_dict[dec_level - 1] - indices = np.where(aux_data.decimation_levels == dec_level)[0] - for ndx in indices: - period = aux_data.periods[ndx] - # find the nearest period in aurora_tf - aurora_ndx = np.argmin(np.abs(aurora_tf.periods - period)) - rxy[ndx] = aurora_tf.rho[aurora_ndx, 0] - ryx[ndx] = aurora_tf.rho[aurora_ndx, 1] - pxy[ndx] = aurora_tf.phi[aurora_ndx, 0] - pyx[ndx] = aurora_tf.phi[aurora_ndx, 1] - - result = {} - result["rho"] = {} - result["phi"] = {} - result["rho"]["xy"] = rxy - result["phi"]["xy"] = pxy - result["rho"]["yx"] = ryx - result["phi"]["yx"] = pyx - return result - - def clip_bands_from_z_file( z_path: Union[str, pathlib.Path], n_bands_clip: int, diff --git a/aurora/transfer_function/regression/RME.py b/aurora/transfer_function/regression/RME.py index 4efe8b4b..b5445dce 100644 --- a/aurora/transfer_function/regression/RME.py +++ b/aurora/transfer_function/regression/RME.py @@ -72,9 +72,10 @@ import numpy as np import xarray as xr -from aurora.transfer_function.regression.m_estimator import MEstimator from scipy.linalg import solve_triangular +from aurora.transfer_function.regression.m_estimator import MEstimator + class RME(MEstimator): def __init__(self, **kwargs): @@ -96,15 +97,15 @@ def __init__(self, **kwargs): Psi' is something like 1 between -r0, and r0 Psi' is zero outside So the expectation value of psi' is the number of points outside - Its the number of points that didn't get weighted /total number of points + Its the number of points that didn't get weighted divided by the total number of points. """ super(RME, self).__init__(**kwargs) self.qr_input = "X" def update_y_hat(self) -> None: - """ - Update the predicted_data $\hat{Y}$ + r""" + Update the predicted_data :math:`\hat{Y}` """ self._Y_hat = self.Q @ self.QHYc @@ -115,9 +116,9 @@ def update_residual_variance(self, correction_factor=1): return self._residual_variance def update_b(self) -> None: - """ + r""" Solve the regression problem using the latest cleaned data via QR-decompostion. - matlab was: b = R\QTY; + The matlab code for this was: b = R\\QTY; """ self.b = solve_triangular(self.R, self.QHYc) diff --git a/aurora/transfer_function/regression/RME_RR.py b/aurora/transfer_function/regression/RME_RR.py index 1b282ee8..a2582282 100644 --- a/aurora/transfer_function/regression/RME_RR.py +++ b/aurora/transfer_function/regression/RME_RR.py @@ -10,9 +10,10 @@ """ import numpy as np import xarray as xr -from aurora.transfer_function.regression.m_estimator import MEstimator from loguru import logger +from aurora.transfer_function.regression.m_estimator import MEstimator + class RME_RR(MEstimator): def __init__(self, **kwargs): @@ -59,12 +60,12 @@ def check_reference_data_shape(self) -> None: raise Exception def update_y_hat(self) -> None: - """updates the predicted data""" + """Updates the predicted data.""" self._Y_hat = self.X @ self.b def update_b(self) -> None: - """ - Updates the tf estimate data + r""" + Updates the tf estimate data. matlab code was: b = QTX\QTY """ @@ -87,7 +88,7 @@ def update_residual_variance(self, correction_factor=1) -> np.ndarray: return self._residual_variance def compute_inverse_signal_covariance(self) -> xr.DataArray: - """ + r""" Computes the inverse signal covariance matrix of the input channels. Development Notes: diff --git a/aurora/transfer_function/regression/base.py b/aurora/transfer_function/regression/base.py index ce9da4fb..6b0d8e20 100644 --- a/aurora/transfer_function/regression/base.py +++ b/aurora/transfer_function/regression/base.py @@ -1,20 +1,19 @@ """ -This module contains the base class for regression functions. - It follows Gary Egbert's EMTF Matlab code TRegression.m in - which can be found in +This module contains the base class for regression functions. It follows Gary Egbert's EMTF Matlab code TRegression.m in which can be found in iris_mt_scratch/egbert_codes-20210121T193218Z-001/egbert_codes/matlabPrototype_10-13-20/TF/classes This class originally used numpy arrays to make adapting the Matlab easier, but experimental support for xarray is now added (2024). """ +from typing import Optional, Union + import numpy as np import xarray as xr +from loguru import logger from aurora.transfer_function.regression.iter_control import IterControl -from loguru import logger -from typing import Optional, Union class RegressionEstimator(object): @@ -100,10 +99,10 @@ def __init__( iter_control: IterControl = IterControl(), input_channel_names: Optional[Union[list, None]] = None, output_channel_names: Optional[Union[list, None]] = None, - **kwargs, # n + **kwargs, ): """ - Constructor + Constructor. Parameters ---------- @@ -118,8 +117,6 @@ def __init__( output_channel_names: Optional[Union[list, None]] If Y is np.ndarray, this allows associating channel names to Y's columns - ---------- - kwargs """ self._X = X self._Y = Y @@ -144,14 +141,10 @@ def _set_up_regression_variables(self) -> None: """ Initialize arrays needed for regression and cast any xarray to numpy - Development Notes: + **Development Notes**: - When xr.Datasets are X, Y we cast to array (num channels x num observations) and then transpose them - When xr.DataArrays are X, Y extract the array -- but how do we know whether or not to transpose? - - it would be nice to have a helper function that applies the logic of getting the data from the - xarray and transposing or not appropriately. - - for now we assume that the input data are organized so that input arrays are (n_ch x n_observations). - This assumption is OK for xr.Dataset where the datavars are the MT components ("hx", "hy", etc) + When self._X,_Y are xr.Datasets they are cast to arrays (num channels x num observations) and then transposed. + When self._X,_Y are xr.DataArrays we extract the array -- but how do we know whether or not to transpose? It would be nice to have a helper function that applies the logic of getting the data from the xarray and transposing or not appropriately. For now we assume that the input data are organized so that input arrays are (n_ch x n_observations). This assumption is OK for xr.Dataset where the datavars are the MT components ("hx", "hy", etc) """ self.X = _input_to_numpy_array(self._X) @@ -303,8 +296,10 @@ def n_channels_out(self) -> int: @property def degrees_of_freedom(self) -> int: """ - gets the number of degress of freedom in the dataset. + Gets the number of degrees of freedom in the dataset. + Returns + ------- int The total number of multivariate observations minus the number of input channels. """ diff --git a/aurora/transfer_function/regression/iter_control.py b/aurora/transfer_function/regression/iter_control.py index dea77021..91995f28 100644 --- a/aurora/transfer_function/regression/iter_control.py +++ b/aurora/transfer_function/regression/iter_control.py @@ -5,8 +5,8 @@ Based on Gary's IterControl.m in iris_mt_scratch/egbert_codes-20210121T193218Z-001/egbert_codes/matlabPrototype_10-13-20/TF/classes """ -from loguru import logger import numpy as np +from loguru import logger from aurora.transfer_function.regression.helper_functions import rme_beta @@ -14,15 +14,12 @@ class IterControl(object): """ Notes: - Here is a class to hold variables that are used to control the regression - - currently this is used for variations on RME (Robust M-Estimator) - - TODO: in the original matlab code there was a class property called `epsilon`, but - it was unused; There was no documentation about it's purpose, except possibly that - the abstract base class solved Y = X*b + epsilon for b, complex-valued. - Perhaps this was intended as an intrinsic tolerated noise level. The value of - epsilon was set to 1000. - - TODO The return covariance boolean just initializes arrays of zeros. Needs to be - made functional or removed + + Here is a class to hold variables that are used to control the regression: + + - currently this is used for variations on RME (Robust M-Estimator) + - TODO: in the original matlab code there was a class property called `epsilon`, but it was unused; There was no documentation about its purpose, except possibly that the abstract base class solved Y = X*b + epsilon for b, complex-valued. Perhaps this was intended as an intrinsic tolerated noise level. The value of epsilon was set to 1000. + - TODO: The return covariance boolean just initializes arrays of zeros. Needs to be made functional or removed. """ def __init__( diff --git a/aurora/transfer_function/regression/m_estimator.py b/aurora/transfer_function/regression/m_estimator.py index 8b4d1fc2..bb23da4a 100644 --- a/aurora/transfer_function/regression/m_estimator.py +++ b/aurora/transfer_function/regression/m_estimator.py @@ -5,12 +5,13 @@ See Notes in RME, RME_RR for more details. """ +from copy import deepcopy + import numpy as np import xarray as xr +from loguru import logger from aurora.transfer_function.regression.base import RegressionEstimator -from copy import deepcopy -from loguru import logger class MEstimator(RegressionEstimator): @@ -246,13 +247,16 @@ def apply_redecending_influence_function(self) -> None: def estimate(self) -> None: """ - Executes the regression + Executes the regression. Development Notes: - Here is a comment from the matlab codes: - "need to look at how we should compute adjusted residual cov to make - consistent with tranmt" - See issue#69 aurora github repo addresses this + + Here is a comment from the matlab codes:: + + need to look at how we should compute adjusted residual cov to make + consistent with tranmt + + See issue #69 in the aurora GitHub repo addresses this. """ if self.is_underdetermined: self.solve_underdetermined() @@ -310,23 +314,25 @@ def update_y_cleaned_via_redescend_weights(self) -> None: def compute_squared_coherence(self) -> None: """ - Updates the array self.R2 + Updates the array self.R2. Here is taken the ratio of the energy in the residuals with the energy in the cleaned data. This metric can be interpreted as how much of the signal (Y) is "explained" by the regression. Development Notes: + The matlab code (TRME_RR) claimed: + % R2 is squared coherence (top row is using raw data, bottom % cleaned, with crude correction for amount of down-weighted data) - TODO: There seem to be other valid metrics for this sort of quantity. In particular, we may want to - consider SSY (the sum of squares of the observed data) over SSR. + TODO: There seem to be other valid metrics for this sort of quantity. In particular, we may want to consider SSY (the sum of squares of the observed data) over SSR. - TODO: consider renaming self.R2. That name invokes the idea of the squared residuals. That is not what - is being stored in self.R2. This is more like a CMRR. + TODO: consider renaming self.R2. That name invokes the idea of the squared residuals. That is not what is being stored in self.R2. This is more like a CMRR. - res: Residuals, the original data minus the predicted data. + Parameters + ---------- + res : Residuals, the original data minus the predicted data. SSR : Sum of squares of the residuals, per channel """ diff --git a/aurora/transfer_function/weights/edf_weights.py b/aurora/transfer_function/weights/edf_weights.py index ce1fe4a3..b6b4a693 100644 --- a/aurora/transfer_function/weights/edf_weights.py +++ b/aurora/transfer_function/weights/edf_weights.py @@ -9,10 +9,11 @@ """ +from typing import Optional + import numpy as np import xarray as xr from loguru import logger -from typing import Optional, Union class EffectiveDegreesOfFreedom(object): @@ -108,8 +109,7 @@ def compute_weights(self, X: np.ndarray, use: np.ndarray) -> np.ndarray: - The covariance matrix `S` is computed as the outer product of the selected data. - The covariance matrix is normalized by the number of selected observations. - The inverse covariance matrix `H` is computed. - - The edf weights are then calculated using the diagonal and off-diagonal elements - of the inverse covariance matrix `H` and the selected data. + - The edf weights are then calculated using the diagonal and off-diagonal elements of the inverse covariance matrix `H` and the selected data. ---Statistical Context--- The inverse covariance matrix (H = S^{-1}) is used here because it "whitens" the data: @@ -128,8 +128,8 @@ def compute_weights(self, X: np.ndarray, use: np.ndarray) -> np.ndarray: the EDF weights reflect true statistical leverage. A note on usage of real vs complex data: - The terms ( X[0, :] * \conj{X[0, :]} ) and ( X[1, :] * \conj{X[1, :]} ) are always real - and non-negative (they are squared magnitudes). The cross term ( \conj{X[1, :]} * X[0, :] ) + The terms ( X[0, :] * \\conj{X[0, :]} ) and ( X[1, :] * \\conj{X[1, :]} ) are always real + and non-negative (they are squared magnitudes). The cross term ( \\conj{X[1, :]} * X[0, :] ) can be complex, but in the context of covariance and quadratic forms you want the real part, not the absolute value. Why? diff --git a/cliff.toml b/cliff.toml new file mode 100644 index 00000000..ab27d749 --- /dev/null +++ b/cliff.toml @@ -0,0 +1,36 @@ +[changelog] +header = """# Changelog + +All notable changes to this project will be documented in this file. +""" + +body = """ +{% for group, commits in commits | group_by(attribute="group") %} + ## {{ group | upper_first }} + {% for commit in commits %} + - {% if commit.breaking %}[**BREAKING**] {% endif %}{{ commit.message | upper_first }} ([{{ commit.id | truncate(length=7, end="") }}]({{ commit.id }})) + {% endfor %} +{% endfor %} +""" + +trim = true + +[git] +conventional_commits = true +filter_unconventional = false +commit_parsers = [ + { message = "^feat", group = "Features" }, + { message = "^fix", group = "Bug Fixes" }, + { message = "^doc", group = "Documentation" }, + { message = "^perf", group = "Performance" }, + { message = "^refactor", group = "Refactor" }, + { message = "^style", group = "Styling" }, + { message = "^test", group = "Testing" }, + { message = "^chore\\(release\\): prepare for", skip = true }, + { message = "^chore\\(deps\\)", skip = true }, + { message = "^chore\\(pr\\)", skip = true }, + { message = "^chore\\(pull\\)", skip = true }, + { message = "^chore", group = "Miscellaneous Tasks" }, + { body = ".*security", group = "Security" }, + { message = "^revert", group = "Revert" }, +] diff --git a/data/synthetic/emtf_results/from_matlab_256_26.zss b/data/synthetic/emtf_results/from_matlab_256_26.zss index e7e40edb..da98c7b0 100644 --- a/data/synthetic/emtf_results/from_matlab_256_26.zss +++ b/data/synthetic/emtf_results/from_matlab_256_26.zss @@ -2,7 +2,7 @@ ********** WITH FULL ERROR COVARINCE********** Robust Single station station :test1 -coordinate 1007.996 0.0 declination 0.0 +coordinate 37.996 0.0 declination 0.0 number of channels 5 number of frequencies 26 orientations and tilts of each channel 1 0.00 0.00 tes Hx diff --git a/data/synthetic/emtf_results/test2.zss b/data/synthetic/emtf_results/test2.zss index dd9e65da..6978cf02 100644 --- a/data/synthetic/emtf_results/test2.zss +++ b/data/synthetic/emtf_results/test2.zss @@ -2,7 +2,7 @@ ********** WITH FULL ERROR COVARINCE********** Robust Single station station :test2 -coordinate 2951.999 0.000 declination 0.00 +coordinate 31.999 0.000 declination 0.00 number of channels 5 number of frequencies 25 orientations and tilts of each channel 1 0.00 0.00 tes Hx diff --git a/data/synthetic/emtf_results/test2r1.zrr b/data/synthetic/emtf_results/test2r1.zrr index 1bed707e..a51781d7 100644 --- a/data/synthetic/emtf_results/test2r1.zrr +++ b/data/synthetic/emtf_results/test2r1.zrr @@ -2,7 +2,7 @@ ********** WITH FULL ERROR COVARINCE********** Robust Remote Reference EMAP station :test2r1 -coordinate 1007.996 102.190 declination 0.00 +coordinate 37.996 102.190 declination 0.00 number of channels 5 number of frequencies 25 orientations and tilts of each channel 1 0.00 0.00 tes Hx diff --git a/docs/api/aurora.config.emtf_band_setup.rst b/docs/api/aurora.config.emtf_band_setup.rst deleted file mode 100644 index cbd6d23c..00000000 --- a/docs/api/aurora.config.emtf_band_setup.rst +++ /dev/null @@ -1,10 +0,0 @@ -aurora.config.emtf\_band\_setup package -======================================= - -Module contents ---------------- - -.. automodule:: aurora.config.emtf_band_setup - :members: - :undoc-members: - :show-inheritance: diff --git a/docs/api/aurora.config.rst b/docs/api/aurora.config.rst index f08a5a89..7b963248 100644 --- a/docs/api/aurora.config.rst +++ b/docs/api/aurora.config.rst @@ -7,7 +7,6 @@ Subpackages .. toctree:: :maxdepth: 4 - aurora.config.emtf_band_setup aurora.config.metadata Submodules diff --git a/docs/api/aurora.pipelines.rst b/docs/api/aurora.pipelines.rst index 4d26934f..751dba48 100644 --- a/docs/api/aurora.pipelines.rst +++ b/docs/api/aurora.pipelines.rst @@ -4,14 +4,6 @@ aurora.pipelines package Submodules ---------- -aurora.pipelines.fourier\_coefficients module ---------------------------------------------- - -.. automodule:: aurora.pipelines.fourier_coefficients - :members: - :undoc-members: - :show-inheritance: - aurora.pipelines.helpers module ------------------------------- diff --git a/docs/api/aurora.rst b/docs/api/aurora.rst index 95e0dba6..636ee1a2 100644 --- a/docs/api/aurora.rst +++ b/docs/api/aurora.rst @@ -17,14 +17,6 @@ Subpackages Submodules ---------- -aurora.general\_helper\_functions module ----------------------------------------- - -.. automodule:: aurora.general_helper_functions - :members: - :undoc-members: - :show-inheritance: - Module contents --------------- diff --git a/docs/api/aurora.sandbox.io_helpers.rst b/docs/api/aurora.sandbox.io_helpers.rst index c1154b4d..d6b99e11 100644 --- a/docs/api/aurora.sandbox.io_helpers.rst +++ b/docs/api/aurora.sandbox.io_helpers.rst @@ -36,14 +36,6 @@ aurora.sandbox.io\_helpers.make\_mth5\_helpers module :undoc-members: :show-inheritance: -aurora.sandbox.io\_helpers.zfile\_murphy module ------------------------------------------------ - -.. automodule:: aurora.sandbox.io_helpers.zfile_murphy - :members: - :undoc-members: - :show-inheritance: - Module contents --------------- diff --git a/docs/api/aurora.sandbox.rst b/docs/api/aurora.sandbox.rst index ec222145..e2060879 100644 --- a/docs/api/aurora.sandbox.rst +++ b/docs/api/aurora.sandbox.rst @@ -12,14 +12,6 @@ Subpackages Submodules ---------- -aurora.sandbox.debug\_mt\_metadata\_issue\_85 module ----------------------------------------------------- - -.. automodule:: aurora.sandbox.debug_mt_metadata_issue_85 - :members: - :undoc-members: - :show-inheritance: - aurora.sandbox.mth5\_channel\_summary\_helpers module ----------------------------------------------------- diff --git a/docs/api/aurora.test_utils.synthetic.rst b/docs/api/aurora.test_utils.synthetic.rst index 022f790c..6173d1d5 100644 --- a/docs/api/aurora.test_utils.synthetic.rst +++ b/docs/api/aurora.test_utils.synthetic.rst @@ -4,14 +4,6 @@ aurora.test\_utils.synthetic package Submodules ---------- -aurora.test\_utils.synthetic.make\_mth5\_from\_asc module ---------------------------------------------------------- - -.. automodule:: aurora.test_utils.synthetic.make_mth5_from_asc - :members: - :undoc-members: - :show-inheritance: - aurora.test\_utils.synthetic.make\_processing\_configs module ------------------------------------------------------------- @@ -36,22 +28,6 @@ aurora.test\_utils.synthetic.processing\_helpers module :undoc-members: :show-inheritance: -aurora.test\_utils.synthetic.rms\_helpers module ------------------------------------------------- - -.. automodule:: aurora.test_utils.synthetic.rms_helpers - :members: - :undoc-members: - :show-inheritance: - -aurora.test\_utils.synthetic.station\_config module ---------------------------------------------------- - -.. automodule:: aurora.test_utils.synthetic.station_config - :members: - :undoc-members: - :show-inheritance: - Module contents --------------- diff --git a/docs/api/aurora.time_series.rst b/docs/api/aurora.time_series.rst index 67fcd43e..64f244d8 100644 --- a/docs/api/aurora.time_series.rst +++ b/docs/api/aurora.time_series.rst @@ -28,14 +28,6 @@ aurora.time\_series.frequency\_band\_helpers module :undoc-members: :show-inheritance: -aurora.time\_series.spectrogram module --------------------------------------- - -.. automodule:: aurora.time_series.spectrogram - :members: - :undoc-members: - :show-inheritance: - aurora.time\_series.time\_axis\_helpers module ---------------------------------------------- diff --git a/docs/api/aurora.transfer_function.plot.rst b/docs/api/aurora.transfer_function.plot.rst deleted file mode 100644 index 663ed701..00000000 --- a/docs/api/aurora.transfer_function.plot.rst +++ /dev/null @@ -1,45 +0,0 @@ -aurora.transfer\_function.plot package -====================================== - -Submodules ----------- - -aurora.transfer\_function.plot.comparison\_plots module -------------------------------------------------------- - -.. automodule:: aurora.transfer_function.plot.comparison_plots - :members: - :undoc-members: - :show-inheritance: - -aurora.transfer\_function.plot.error\_bar\_helpers module ---------------------------------------------------------- - -.. automodule:: aurora.transfer_function.plot.error_bar_helpers - :members: - :undoc-members: - :show-inheritance: - -aurora.transfer\_function.plot.rho\_phi\_helpers module -------------------------------------------------------- - -.. automodule:: aurora.transfer_function.plot.rho_phi_helpers - :members: - :undoc-members: - :show-inheritance: - -aurora.transfer\_function.plot.rho\_plot module ------------------------------------------------ - -.. automodule:: aurora.transfer_function.plot.rho_plot - :members: - :undoc-members: - :show-inheritance: - -Module contents ---------------- - -.. automodule:: aurora.transfer_function.plot - :members: - :undoc-members: - :show-inheritance: diff --git a/docs/api/aurora.transfer_function.regression.rst b/docs/api/aurora.transfer_function.regression.rst index 8c75b213..94b4a964 100644 --- a/docs/api/aurora.transfer_function.regression.rst +++ b/docs/api/aurora.transfer_function.regression.rst @@ -20,13 +20,6 @@ aurora.transfer\_function.regression.RME\_RR module :undoc-members: :show-inheritance: -aurora.transfer\_function.regression.base module ------------------------------------------------- - -.. automodule:: aurora.transfer_function.regression.base - :members: - :undoc-members: - :show-inheritance: aurora.transfer\_function.regression.helper\_functions module ------------------------------------------------------------- diff --git a/docs/api/aurora.transfer_function.rst b/docs/api/aurora.transfer_function.rst index 6352add2..7ba312b1 100644 --- a/docs/api/aurora.transfer_function.rst +++ b/docs/api/aurora.transfer_function.rst @@ -7,7 +7,6 @@ Subpackages .. toctree:: :maxdepth: 4 - aurora.transfer_function.plot aurora.transfer_function.regression aurora.transfer_function.weights @@ -38,14 +37,6 @@ aurora.transfer\_function.emtf\_z\_file\_helpers module :undoc-members: :show-inheritance: -aurora.transfer\_function.kernel\_dataset module ------------------------------------------------- - -.. automodule:: aurora.transfer_function.kernel_dataset - :members: - :undoc-members: - :show-inheritance: - aurora.transfer\_function.transfer\_function\_collection module --------------------------------------------------------------- diff --git a/docs/conf.py b/docs/conf.py index 28772325..a003e347 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -14,9 +14,9 @@ # import sys # sys.path.insert(0, os.path.abspath('.')) -import aurora from sphinx_gallery.sorting import FileNameSortKey + # -- Project information ----------------------------------------------------- project = "aurora" @@ -38,6 +38,7 @@ "sphinx.ext.autosummary", "sphinx.ext.coverage", "sphinx.ext.doctest", + "sphinx.ext.mathjax", "sphinx.ext.extlinks", "sphinx.ext.intersphinx", "sphinx.ext.mathjax", @@ -76,14 +77,21 @@ # a list of builtin themes. # try: - import sphinx_rtd_theme + pass html_theme = "sphinx_rtd_theme" - html_theme_path = [sphinx_rtd_theme.get_html_theme_path()] - pass + except Exception: html_theme = "default" +# Napoleon settings for docstring parsing +napoleon_use_rtype = True +napoleon_use_param = True +napoleon_include_init_with_doc = True +napoleon_include_private_with_doc = True +napoleon_include_special_with_doc = True + + # Add any paths that contain custom static files (such as style sheets) here, # relative to this directory. They are copied after the builtin static files, # so a file named "default.css" will overwrite the builtin "default.css". diff --git a/docs/tutorials/earthscope_magnetic_data_tutorial.ipynb b/docs/tutorials/earthscope_magnetic_data_tutorial.ipynb index 04a41902..a224794c 100644 --- a/docs/tutorials/earthscope_magnetic_data_tutorial.ipynb +++ b/docs/tutorials/earthscope_magnetic_data_tutorial.ipynb @@ -181,24 +181,25 @@ "source": [ "## FDSN formats\n", "\n", - "Earthscope is one of several organizations that use standards defined by the International Federation of Digital Seismograph Networks (FDSN) which has specific conventions for specifying data time series. \n", + "Earthscope is one of several organizations that use standards defined by the International Federation of Digital Seismograph Networks (FDSN) which has specific conventions for specifying data time series.\n", "\n", - "The FDSN gives detailed information of their \"SEED\" convention in the [SEED manual](http://www.fdsn.org/pdf/SEEDManual_V2.4_Appendix-A.pdf) and the interested reader is referred to that document for details on conventions, including an Appendix on channel naming conventions. \n", + "The FDSN gives detailed information of their \"SEED\" convention in the [SEED manual](http://www.fdsn.org/pdf/SEEDManual_V2.4_Appendix-A.pdf) and the interested reader is referred to that document for details on conventions, including an Appendix on channel naming conventions.\n", "\n", "For the example here, we need only know two things:\n", - "1. The Earthscope data collection is organized as a heirarchy, with the top level container being a \"network\", which in turn contains \"stations\" (physical observatories), which in turn contain \"channels\" which correspond to individual data streams from instruments. \n", - "2. that a 3-character code is used to specify a data stream, where the first character specifies sample rate, the second corresponds to the field being measured, and the third character indicates the orientation of the field.\n", + "\n", + "1. The Earthscope data collection is organized as a hierarchy, with the top level container being a \"network\", which in turn contains \"stations\" (physical observatories), which in turn contain \"channels\" which correspond to individual data streams from instruments.\n", + "2. A 3-character code is used to specify a data stream, where the first character specifies sample rate, the second corresponds to the field being measured, and the third character indicates the orientation of the field.\n", "\n", "The main features of the interface will be illustrated for the purpose of magnetometer data.\n", - " - Magnetic data is associated with the SEED Instrument Code F, \n", - " - The SEED manual specifies that orientations for magnetometer data are [\"N\", \"E\", \"Z\"], however, in recent years the codes [\"1\", \"2\", \"3\"] have also been used.\n", - " - Many sample rates for magnetic field data are used, so we will use the wildcard character for sample rate ('\\*', '?').\n", - " - The '?' wildcard replaces a single character, whereas the '\\*' wildcard can replace mulitple values. We use '?' here.\n", + "- Magnetic data is associated with the SEED Instrument Code F.\n", + "- The SEED manual specifies that orientations for magnetometer data are [\"N\", \"E\", \"Z\"], however, in recent years the codes [\"1\", \"2\", \"3\"] have also been used.\n", + "- Many sample rates for magnetic field data are used, so we will use the wildcard character for sample rate ('*', '?').\n", + " - The '?' wildcard replaces a single character, whereas the '*' wildcard can replace multiple values. We use '?' here.\n", "\n", - " Notes:\n", - " - Using wildcards in the third (field component) position is not advised as the \"F\" designator can also be associated with station health metrics such \"VFP\" as Packet Buffer Usage.\n", - " - The \"Channel\" field in the map accepts lists. Thus, to get a sense of all the magnetometer data, one can place in the the channel field \"?F1, ?F2, ?F3, ?FN, ?FE, ?FZ\" as in the example below.\n", - " - Most Magnetometer measurements are 3-component, so it is normally sufficient to query one numeric style and on letter style component, for example [\"?F3, ?FZ\"] to get a sense of data coverage." + "Notes:\n", + "- Using wildcards in the third (field component) position is not advised as the \"F\" designator can also be associated with station health metrics such as \"VFP\" (Packet Buffer Usage).\n", + "- The \"Channel\" field in the map accepts lists. Thus, to get a sense of all the magnetometer data, one can place in the channel field \"?F1, ?F2, ?F3, ?FN, ?FE, ?FZ\" as in the example below.\n", + "- Most Magnetometer measurements are 3-component, so it is normally sufficient to query one numeric style and one letter style component, for example [\"?F3, ?FZ\"] to get a sense of data coverage.\n" ] }, { diff --git a/docs/tutorials/synthetic_data_processing.ipynb b/docs/tutorials/synthetic_data_processing.ipynb index a69e8cac..5610d90a 100644 --- a/docs/tutorials/synthetic_data_processing.ipynb +++ b/docs/tutorials/synthetic_data_processing.ipynb @@ -1671,7 +1671,7 @@ "\n", "# Feature Storage **Experimental work in progress**\n", "\n", - "### Make an example of a feature embedded in mth5 file from the FCs.\n", + "## Make an example of a feature embedded in mth5 file from the FCs.\n", "\n", "Any spectrogram (FC array at some decimation level) can be passed to a feature extraction method.\n", "These features can, in turn be stored in the mth5.\n", diff --git a/pyproject.toml b/pyproject.toml index b89e313e..3c467c23 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,7 +5,7 @@ build-backend = "setuptools.build_meta" [project] name = "aurora" -version = "0.5.2" +version = "0.6.0" description = "Processing Codes for Magnetotelluric Data" readme = "README.rst" requires-python = ">=3.10" @@ -74,6 +74,21 @@ dev = [ "sphinx_rtd_theme", ] +[tool.bumpversion] +current_version = "0.6.0" +parse = "(?P\\d+)\\.(?P\\d+)\\.(?P\\d+)" +serialize = ["{major}.{minor}.{patch}"] +commit = true +tag = true + +[[tool.bumpversion.files]] +filename = "pyproject.toml" +search = 'version = "{current_version}"' +replace = 'version = "{new_version}"' + +[[tool.bumpversion.files]] +filename = "aurora/__init__.py" + diff --git a/pytest.ini b/pytest.ini index 2fce28ba..d27c71b5 100644 --- a/pytest.ini +++ b/pytest.ini @@ -5,8 +5,6 @@ testpaths = tests python_files = test_*.py python_classes = Test* python_functions = test_* -timeout = 300 -timeout_method = thread filterwarnings = ignore:Pydantic serializer warnings:UserWarning ignore:.*Jupyter is migrating its paths to use standard platformdirs.*:DeprecationWarning diff --git a/tests/cas04/01_make_cas04_mth5.py b/tests/cas04/01_make_cas04_mth5.py index 3a27debc..57e55eae 100644 --- a/tests/cas04/01_make_cas04_mth5.py +++ b/tests/cas04/01_make_cas04_mth5.py @@ -25,14 +25,12 @@ """ import pandas as pd - -from aurora.general_helper_functions import get_test_path -from aurora.general_helper_functions import execute_subprocess -from aurora.sandbox.mth5_helpers import build_request_df - +from loguru import logger from mth5.clients import FDSN from mth5.utils.helpers import read_back_data -from loguru import logger + +from aurora.general_helper_functions import execute_subprocess, get_test_path +from aurora.sandbox.mth5_helpers import build_request_df # Define paths @@ -103,7 +101,6 @@ def make_all_stations_individually( mth5_version: """ for station_id in STATION_IDS: - # request_df = build_request_df(NETWORK_ID, station_id, channels=["*F*", "*Q*", ], start=None, end=None) request_df = build_request_df( NETWORK_ID, station_id, channels=CHANNELS, start=None, end=None ) diff --git a/tests/synthetic/test_compare_aurora_vs_archived_emtf_pytest.py b/tests/synthetic/test_compare_aurora_vs_archived_emtf_pytest.py index 3eb85dc2..65b41d4b 100644 --- a/tests/synthetic/test_compare_aurora_vs_archived_emtf_pytest.py +++ b/tests/synthetic/test_compare_aurora_vs_archived_emtf_pytest.py @@ -1,3 +1,19 @@ +""" +Tests comparing aurora processing results against archived EMTF results. + +Development Notes: +- In the early days of development these tests were useful to check that +the same results were obtained from processing an mth5 with both stations in it +and two separate mth5 files. This could probably be made into a simpler test in mth5 +that checks that the data are the same in the two files. +- This used to make a homebrew resistivity and phase plot for comparison between +archived emtf z-files, but has been replaced with mt_metadata methods. +- TODO: Check phases in these plots -- they are off by 180 so there may be a sign error +in the data, or maybe the emtf results are using a different convention. Need to investigate. +- The comparison with the matlab emtf results uses a slighly different windowing method. + +""" +import numpy as np import pytest from loguru import logger from mth5.helpers import close_open_files @@ -5,17 +21,8 @@ from aurora.general_helper_functions import DATA_PATH from aurora.pipelines.process_mth5 import process_mth5 -from aurora.sandbox.io_helpers.zfile_murphy import read_z_file from aurora.test_utils.synthetic.make_processing_configs import create_test_run_config -from aurora.test_utils.synthetic.plot_helpers_synthetic import plot_rho_phi -from aurora.test_utils.synthetic.rms_helpers import ( - assert_rms_misfit_ok, - compute_rms, - get_expected_rms_misfit, -) -from aurora.transfer_function.emtf_z_file_helpers import ( - merge_tf_collection_to_match_z_file, -) +from aurora.transfer_function.compare import CompareTF # Path to baseline EMTF results in source tree @@ -29,9 +36,8 @@ def aurora_vs_emtf( auxilliary_z_file, z_file_base, tfk_dataset, + atol_phase=4.0, make_rho_phi_plot=True, - show_rho_phi_plot=False, - use_subtitle=True, ): """ Compare aurora processing results against EMTF baseline. @@ -60,52 +66,27 @@ def aurora_vs_emtf( test_case_id, tfk_dataset, matlab_or_fortran=emtf_version ) - expected_rms_misfit = get_expected_rms_misfit(test_case_id, emtf_version) z_file_path = AURORA_RESULTS_PATH.joinpath(z_file_base) - tf_collection = process_mth5( + process_mth5( processing_config, tfk_dataset=tfk_dataset, z_file_path=z_file_path, return_collection=True, ) - - aux_data = read_z_file(auxilliary_z_file) - aurora_rho_phi = merge_tf_collection_to_match_z_file(aux_data, tf_collection) - data_dict = {} - data_dict["period"] = aux_data.periods - data_dict["emtf_rho_xy"] = aux_data.rxy - data_dict["emtf_phi_xy"] = aux_data.pxy - for xy_or_yx in ["xy", "yx"]: - aurora_rho = aurora_rho_phi["rho"][xy_or_yx] - aurora_phi = aurora_rho_phi["phi"][xy_or_yx] - aux_rho = aux_data.rho(xy_or_yx) - aux_phi = aux_data.phi(xy_or_yx) - rho_rms_aurora, phi_rms_aurora = compute_rms( - aurora_rho, aurora_phi, verbose=True + comparator = CompareTF(tf_01=z_file_path, tf_02=auxilliary_z_file) + result = comparator.compare_transfer_functions(atol_phase=4.0) + assert np.isclose(result["impedance_ratio"]["Z_10"], 1.0, atol=1e-2) + assert np.isclose(result["impedance_ratio"]["Z_01"], 1.0, atol=1e-2) + assert result["impedance_phase_close"] + if make_rho_phi_plot: + comparator.plot_two_transfer_functions( + save_plot_path=AURORA_RESULTS_PATH.joinpath( + f"{test_case_id}_aurora_vs_emtf_{emtf_version}_tf_compare.png" + ), + rho_xy_ylims=(10.0, 1000.0), + rho_yx_ylims=(10.0, 1000.0), ) - rho_rms_emtf, phi_rms_emtf = compute_rms(aux_rho, aux_phi) - data_dict["aurora_rho_xy"] = aurora_rho - data_dict["aurora_phi_xy"] = aurora_phi - if expected_rms_misfit is not None: - assert_rms_misfit_ok( - expected_rms_misfit, xy_or_yx, rho_rms_aurora, phi_rms_aurora - ) - - if make_rho_phi_plot: - plot_rho_phi( - xy_or_yx, - tf_collection, - rho_rms_aurora, - rho_rms_emtf, - phi_rms_aurora, - phi_rms_emtf, - emtf_version, - aux_data=aux_data, - use_subtitle=use_subtitle, - show_plot=show_rho_phi_plot, - output_path=AURORA_RESULTS_PATH, - ) @pytest.mark.slow diff --git a/tests/test_doc_build.py b/tests/test_doc_build.py index 7e3db336..22f2d27e 100644 --- a/tests/test_doc_build.py +++ b/tests/test_doc_build.py @@ -1,13 +1,13 @@ -import subprocess -import unittest import os import platform +import subprocess +import unittest class Doc_Test(unittest.TestCase): @property def path_to_docs(self): - dirname, file_name = os.path.split(os.path.abspath(__file__)) + dirname, _ = os.path.split(os.path.abspath(__file__)) return dirname.split(os.path.sep)[:-1] + ["docs"] def test_html(self): diff --git a/tests/test_general_helper_functions.py b/tests/test_general_helper_functions.py index 91ae895e..64abdb66 100644 --- a/tests/test_general_helper_functions.py +++ b/tests/test_general_helper_functions.py @@ -1,12 +1,8 @@ - import logging import unittest -from aurora.general_helper_functions import count_lines -from aurora.general_helper_functions import DotDict -from aurora.general_helper_functions import get_test_path -from aurora.general_helper_functions import replace_in_file -from loguru import logger +from aurora.general_helper_functions import DotDict, get_test_path, replace_in_file + TEST_PATH = get_test_path() @@ -18,18 +14,6 @@ def setUp(self): logging.getLogger("matplotlib.font_manager").disabled = True logging.getLogger("matplotlib.ticker").disabled = True - def test_count_lines(self): - tmp_file = TEST_PATH.joinpath("tmp.txt") - n_lines_in = 42 - lines = n_lines_in * ["test\n"] - f = open(tmp_file, "w") - f.writelines(lines) - f.close() - n_lines_out = count_lines(tmp_file) - assert n_lines_out == n_lines_in - tmp_file.unlink() - return - def test_dot_dict(self): tmp = {} tmp["a"] = "aa" @@ -37,7 +21,7 @@ def test_dot_dict(self): dot_dict = DotDict(tmp) assert dot_dict.a == tmp["a"] assert dot_dict.b == "bb" - + def test_replace_in_file(self): # Create a temporary file tmp_file = TEST_PATH.joinpath("tmp_replace.txt") diff --git a/tests/time_series/test_apodization_window_pytest.py b/tests/time_series/test_apodization_window_pytest.py index d2a50bdd..d40fab6a 100644 --- a/tests/time_series/test_apodization_window_pytest.py +++ b/tests/time_series/test_apodization_window_pytest.py @@ -215,6 +215,37 @@ def test_window_energy_conservation(self, subtests): assert energy > 0 assert np.isfinite(energy) + def test_linear_spectral_density_factor(self, subtests) -> None: + r""" + This is just a test to verify some algebra + + Claim: + The lsd_calibration factors + A (1./coherent_gain)*np.sqrt((2*dt)/(nenbw*N)) + and + B np.sqrt(2/(sample_rate*self.S2)) + + Note sqrt(2*dt)==sqrt(2*sample_rate) so we can cancel these terms and A=B IFF + + (1./coherent_gain) * np.sqrt(1/(nenbw*N)) == 1/np.sqrt(S2) + which is shown in github aurora issue #3 via (CG**2) * NENBW *N = S2 + + """ + taper_families = ["boxcar", "hamming", "hann", "blackman"] + + for family in taper_families: + with subtests.test(taper_family=family): + window = ApodizationWindow(taper_family=family, num_samples_window=256) + lsd_factor1 = (1.0 / window.coherent_gain) * np.sqrt( + 1.0 / (window.nenbw * window.num_samples_window) + ) + lsd_factor2 = 1.0 / np.sqrt(window.S2) + if not np.isclose(lsd_factor1, lsd_factor2): + msg = f"Linear spectral density factors do not match for {family} window: \n" + msg += f"lsd_factor1 {lsd_factor1} vs lsd_factor2 {lsd_factor2}" + logger.error(msg) + raise Exception(msg) + class TestApodizationWindowParameterVariations: """Test windows with various parameter combinations."""