diff --git a/.github/ISSUE_TEMPLATE.md b/.github/ISSUE_TEMPLATE.md new file mode 100644 index 0000000..bd01eeb --- /dev/null +++ b/.github/ISSUE_TEMPLATE.md @@ -0,0 +1,15 @@ +* synthterrain version: +* Python version: +* Operating System: + +### Description + +Describe what you were trying to get done. +Tell us what happened, what went wrong, and what you expected to happen. + +### What I Did + +``` +Paste the command(s) you ran and the output. +If there was a crash, please include the traceback here. +``` diff --git a/.github/pull_request_template.md b/.github/pull_request_template.md new file mode 100644 index 0000000..5d4da34 --- /dev/null +++ b/.github/pull_request_template.md @@ -0,0 +1,46 @@ + + +## Description + + +## Related Issue + + + + +## Motivation and Context + + + +## How Has This Been Tested? + + + + + +## Types of changes + +- Bug fix (non-breaking change which fixes an issue) +- New feature (non-breaking change which adds functionality) +- Breaking change (fix or feature that would cause existing functionality to change) + +## Checklist: + + +- My change requires a change to the documentation. +- I have updated the documentation accordingly. +- I have added tests to cover my changes. +- All new and existing tests passed. + +## Licensing: + +This project is released under the [LICENSE](https://github.com/NeoGeographyToolkit/synthterrain/blob/master/LICENSE). + + +- I claim copyrights on my contributions in this pull request, and I provide those contributions via this pull request under the same license terms that the synthterrain project uses, and I have submitted a CLA. +- I dedicate any and all copyright interest in my contributions in this pull request to the public domain. I make this dedication for the benefit of the public at large and to the detriment of my heirs and successors. I intend this dedication to be an overt act of relinquishment in perpetuity of all present and future rights to this contribution under copyright law. + + + + diff --git a/.github/workflows/lint.yml b/.github/workflows/lint.yml new file mode 100644 index 0000000..e5e2292 --- /dev/null +++ b/.github/workflows/lint.yml @@ -0,0 +1,30 @@ +name: Lint + +on: [push] + +jobs: + black: + runs-on: ubuntu-latest + steps: + - name: Check out source repository + uses: actions/checkout@v4 + - name: Black Version + uses: psf/black@stable + with: + options: "--version" + - name: Black Check + uses: psf/black@stable + + flake8: + runs-on: ubuntu-latest + steps: + - name: Check out source repository + uses: actions/checkout@v4 + - name: Set up Python environment + uses: actions/setup-python@v5 + with: + python-version: "3.8" + - name: flake8 Lint + uses: py-actions/flake8@v2 + with: + path: "src/python/synthterrain tests/python" diff --git a/.github/workflows/python-test.yml b/.github/workflows/python-test.yml new file mode 100644 index 0000000..3023d0e --- /dev/null +++ b/.github/workflows/python-test.yml @@ -0,0 +1,44 @@ +# This workflow will install Python dependencies, run tests and lint with a variety of Python versions +# For more information see: https://help.github.com/actions/language-and-framework-guides/using-python-with-github-actions + +name: Python Testing + +on: [push] + +jobs: + build: + + runs-on: ubuntu-latest + strategy: + matrix: + python-version: ['3.8', '3.9', '3.10', '3.11'] + + steps: + - name: Checkout + uses: actions/checkout@v4 + with: + ref: ${{ github.event.pull_request.head.sha }} + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python-version }} + - name: Display Python version + run: python -c "import sys; print(sys.version)" + - name: Install dependencies + run: | + python -m pip install --upgrade pip + python -m pip install pytest-cov + python -m pip install -r requirements.txt + - name: Install + # run: python -m pip install -e ${{ matrix.install-target }} + run: python -m pip install -e . + - name: Test with pytest and generate coverage report + run: | + pytest --cov=./src/python/synthterrain --cov-report=xml + + # - name: Upload coverage to Codecov + # uses: codecov/codecov-action@v4 + # with: + # verbose: true + # env: + # CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} diff --git a/AUTHORS.rst b/AUTHORS.rst index 46b9bba..2b04497 100644 --- a/AUTHORS.rst +++ b/AUTHORS.rst @@ -5,11 +5,12 @@ Credits Development Lead ---------------- -* Ross Beyer +* Ross Beyer , SETI Institute & NASA Ames Contributors ------------ * Mark Allan, NASA Ames & KBR +* Scott McMichael, NASA Ames & KBR * Audrow Nash, NASA Ames & KBR -* Jennifer Nguyen, NASA Ames & KBR \ No newline at end of file +* Jennifer Nguyen, NASA Ames & KBR diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 6d184f0..9da950e 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -27,6 +27,9 @@ and the release date, in year-month-day format (see examples below). Unreleased ---------- -0.1.0 (2022-??-??) -* First release. + +0.1.0 (2024-07-02) +------------------ + +- First release. diff --git a/CONTRIBUTING.rst b/CONTRIBUTING.rst index 0a21438..6eea187 100644 --- a/CONTRIBUTING.rst +++ b/CONTRIBUTING.rst @@ -65,40 +65,37 @@ Ready to contribute? Here's how to set up `synthterrain` for local development. 1. Fork the `synthterrain` repo on GitHub. 2. Clone your fork locally:: - $ git clone git@github.com:your_name_here/synthterrain.git + $> git clone git@github.com:your_name_here/synthterrain.git 3. Install your local copy into a virtual environment of your choice (there are many to choose from like conda, etc.). We will assume conda here, but any should work:: - $ cd synthterrain/ - $ conda env create -n synthterrain - $ conda activate synthterrain - $ mamba env update --file environment_dev.yml - $ mamba env update --file environment.yml - $ pip install opensimplex - $ pip install --no-deps -e . + $> cd synthterrain/ + $> conda env create -n synthterrain + $> conda activate synthterrain + $> mamba env update --file environment_dev.yml + $> mamba env update --file environment.yml + $> pip install --no-deps -e . The last ``pip install`` installs synthterrain in "editable" mode which facilitates using the programs and testing. 4. Create a branch for local development:: - $ git checkout -b name-of-your-bugfix-or-feature + $> git checkout -b name-of-your-bugfix-or-feature Now you can make your changes locally. 5. When you're done making changes, check that your changes pass flake8 and the - tests, including testing other Python versions with tox:: + tests. - $ flake8 src/synthterrain tests - $ python setup.py test or pytest - $ tox + $> make lint + $> make test - To get flake8 and tox, just pip install them into your virtual environment. 6. Commit your changes and push your branch to GitHub:: - $ git add . - $ git commit -m "Your detailed description of your changes." - $ git push origin name-of-your-bugfix-or-feature + $> git add . + $> git commit -m "Your detailed description of your changes." + $> git push origin name-of-your-bugfix-or-feature 7. Submit a pull request through the GitHub website. @@ -111,17 +108,107 @@ Before you submit a pull request, check that it meets these guidelines: 2. If the pull request adds functionality, the docs should be updated. Put your new functionality into a function with a docstring, and add the feature to the list in CHANGELOG.rst. -3. The pull request should work for Python 3.7, 3.8, and 3.9 and optionally for PyPy. +3. The pull request should work for Python 3.8, 3.9, 3.10, 3.11 and optionally for PyPy. And make sure that the tests pass for all supported Python versions. - -Deploying ---------- - -A reminder for the maintainers on how to deploy. -Make sure all your changes are committed (including an entry in HISTORY.rst). -Then run:: - -$ bump2version patch # possible: major / minor / patch -$ git push -$ git push --tags +What to expect +-------------- + +Our development of synthterrain is not particularly continuous, +and it is entirely possible that when you submit a PR +(pull request), none of us will have the time to evaluate or integrate +your PR. If we don't, we'll try and communicate that with you via the +PR. + +For large contributions, it is likely that you, or your employer, +will be retaining your copyrights, but releasing the contributions +via an open-source license. It must be compatible with the Apache-2 +license that synthterrain is distributed with, so that we can redistribute +that contribution with synthterrain, give you credit, and make synthterrain even +better! Please contact us if you have a contribution of that nature, +so we can be sure to get all of the details right. + +For smaller contributions, where you (or your employer) are not +concerned about retaining copyright (but we will give you credit!), +you will need to fill out a Contributor License Agreement (CLA) +before we can accept your PR. The CLA assigns your copyright in +your contribution to NASA, so that our NASA copyright statement +remains true: + + Copyright (c) YEAR, United States Government as represented by the + Administrator of the National Aeronautics and Space Administration. + All rights reserved. + +There is an `Individual CLA +`_ and a +`Corporate CLA +`_. + +synthterrain People +------------------- + +- A synthterrain **Contributor** is any individual creating or commenting + on an issue or pull request. Anyone who has authored a PR that was + merged should be listed in the AUTHORS.rst file. + +- A synthterrain **Committer** is a subset of contributors, typically NASA + employees or contractors, who have been given write access to the + repository. + +Rules for Merging Pull Requests +------------------------------- + +Any change to resources in this repository must be through pull +requests (PRs). This applies to all changes to documentation, code, +binary files, etc. Even long term committers must use pull requests. + +In general, the submitter of a PR is responsible for making changes +to the PR. Any changes to the PR can be suggested by others in the +PR thread (or via PRs to the PR), but changes to the primary PR +should be made by the PR author (unless they indicate otherwise in +their comments). In order to merge a PR, it must satisfy these conditions: + +1. Have been open for 24 hours. +2. Have one approval. +3. If the PR has been open for 2 days without approval or comment, then it + may be merged without any approvals. + +Pull requests should sit for at least 24 hours to ensure that +contributors in other timezones have time to review. Consideration +should also be given to weekends and other holiday periods to ensure +active committers all have reasonable time to become involved in +the discussion and review process if they wish. + +In order to encourage involvement and review, we encourage at least +one explicit approval from committers that are not the PR author. + +However, in order to keep development moving along with our low number of +active contributors, if a PR has been open for 2 days without comment, then +it could be committed without an approval. + +The default for each contribution is that it is accepted once no +committer has an objection, and the above requirements are +satisfied. + +In the case of an objection being raised in a pull request by another +committer, all involved committers should seek to arrive at a +consensus by way of addressing concerns being expressed by discussion, +compromise on the proposed change, or withdrawal of the proposed +change. + +Exceptions to the above are minor typo fixes or cosmetic changes +that don't alter the meaning of a document. Those edits can be made +via a PR and the requirement for being open 24 h is waived in this +case. + + +.. Deploying + --------- + + A reminder for the maintainers on how to deploy. + Make sure all your changes are committed (including an entry in CHANGELOG.rst). + Then run:: + + $ bump2version patch # possible: major / minor / patch + $ git push + $ git push --tags diff --git a/LICENSE b/LICENSE index 0508b75..d645695 100644 --- a/LICENSE +++ b/LICENSE @@ -1,7 +1,202 @@ -This package is currently in development and is copyright (C) 2022 -by the authors and the United States Government as represented by -the Administrator of the National Aeronautics and Space Administration. -All rights reserved. -We anticipate release under an open source license when this software -reaches its 1.0 milestone. + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Copyright [yyyy] [name of copyright owner] + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. diff --git a/Makefile b/Makefile index 006b80d..17f7127 100644 --- a/Makefile +++ b/Makefile @@ -56,8 +56,17 @@ clean-test: ## remove test and coverage artifacts rm -fr .pytest_cache rm -fr test-resources -lint: ## check style with flake8 - flake8 --max-complexity 10 --ignore E203,E501,W503 hiproc tests +lint/flake8: ## check style with flake8 + flake8 src/python/synthterrain tests/python/ + +lint/black: ## check style with black + black --check src/python/synthterrain tests + +lint/ufmt: ## check format with ufmt + ufmt check src/python/synthterrain + ufmt check tests/python + +lint: lint/flake8 lint/black lint/ufmt test: ## run tests quickly with the default Python pytest -s @@ -82,13 +91,18 @@ coverage: ## check code coverage quickly with the default Python # servedocs: docs ## compile the docs watching for changes # watchmedo shell-command -p '*.rst' -c '$(MAKE) -C docs html' -R -D . +release-check: dist ## check state of distribution + twine check dist/* + release: dist ## package and upload a release twine upload dist/* dist: clean ## builds source and wheel package - python setup.py sdist - python setup.py bdist_wheel + python -m build ls -l dist +develop: clean ## install the package in an editable format for development + pip install --no-deps -e . + install: clean ## install the package to the active Python's site-packages - python setup.py install + pip install diff --git a/README.rst b/README.rst index abbd168..12d9154 100644 --- a/README.rst +++ b/README.rst @@ -17,7 +17,7 @@ Features -------- The synthterrain package currently offers these command-line programs -when it is installed (see the CONTRIBUTING document). Arguments +when it is installed (see below). Arguments can be found by running any program with a ``-h`` flag. ``synthcraters`` @@ -33,10 +33,10 @@ can be found by running any program with a ``-h`` flag. ``synthrocks``. Also allows a set of pre-existing craters to be added to the probabiliy maps that ``synthrocks`` uses to place rocks. -``craterconvert`` +``synthcraterconvert`` Converts between the current crater CSV and old XML MATLAB formats. -``craterplot`` +``synthcraterplot`` Generates a set of plots from the CSV output of ``synthcraters``. @@ -63,6 +63,31 @@ And then you could call ``synthterrain`` like this:: You can mix regular arguments and ampersand-arguments if you wish. +Installation +------------ + +Clone or download this repository. + +It is highly suggested to install this into a virtual Python environment. + +Change directory to where you have downloaded this repository after you have +set up your virtual environment, just do this:: + +$> pip install + + +or:: + +$> make install + +If you use conda for your virtual environment, you can do this:: + +$> conda create -n synthterrain +$> conda activate synthterrain +$> conda env update --file environment.yml + + + Contributing ------------ @@ -71,3 +96,31 @@ contributing guide for details on how to help and setup a development environment. +Credits +------- + +synthterrain was developed in the open at NASA's Ames Research Center. + +See the `AUTHORS ` +file for a complete list of developers. + + +License +------- +Copyright © 2024, United States Government, as represented by the +Administrator of the National Aeronautics and Space Administration. +All rights reserved. + +The "synthterrain" software is licensed under the Apache License, +Version 2.0 (the "License"); you may not use this file except in +compliance with the License. You may obtain a copy of the License +at http://www.apache.org/licenses/LICENSE-2.0. + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +implied. See the License for the specific language governing +permissions and limitations under the License. + + + diff --git a/THIRDPARTYLICENSES.rst b/THIRDPARTYLICENSES.rst new file mode 100644 index 0000000..375c6ae --- /dev/null +++ b/THIRDPARTYLICENSES.rst @@ -0,0 +1,16 @@ +The synthterrain module would not be possible with out the +use of third party software that is also open source. The following +software is used by synthterrain when run on your system, but is not distributed +with it: + +======================== ============== ===== + Title License URL +======================== ============== ===== +matplotlib custom https://matplotlib.org/stable/users/project/license.html +numpy BSD-3-Clause https://github.com/numpy/numpy/blob/main/LICENSE.txt +opensimplex MIT https://github.com/lmas/opensimplex/blob/master/LICENSE +pandas BSD-3-Clause https://github.com/pandas-dev/pandas/blob/main/LICENSE +rasterio BSD-3-Clause https://github.com/rasterio/rasterio/blob/main/LICENSE.txt +scipy BSD-3-Clause https://github.com/scipy/scipy/blob/main/LICENSE.txt +shapely BSD-3-Clause https://github.com/shapely/shapely/blob/main/LICENSE.txt +======================== ============== ===== diff --git a/docs/Makefile b/docs/Makefile new file mode 100644 index 0000000..d4bb2cb --- /dev/null +++ b/docs/Makefile @@ -0,0 +1,20 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line, and also +# from the environment for the first two. +SPHINXOPTS ?= +SPHINXBUILD ?= sphinx-build +SOURCEDIR = . +BUILDDIR = _build + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +.PHONY: help Makefile + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/docs/authors.rst b/docs/authors.rst new file mode 100644 index 0000000..e122f91 --- /dev/null +++ b/docs/authors.rst @@ -0,0 +1 @@ +.. include:: ../AUTHORS.rst diff --git a/docs/changelog.rst b/docs/changelog.rst new file mode 100644 index 0000000..565b052 --- /dev/null +++ b/docs/changelog.rst @@ -0,0 +1 @@ +.. include:: ../CHANGELOG.rst diff --git a/docs/conf.py b/docs/conf.py new file mode 100644 index 0000000..fa098d6 --- /dev/null +++ b/docs/conf.py @@ -0,0 +1,36 @@ +# Configuration file for the Sphinx documentation builder. +# +# For the full list of built-in configuration values, see the documentation: +# https://www.sphinx-doc.org/en/master/usage/configuration.html + +import os +import sys + +sys.path.insert(0, os.path.abspath("../src/python")) + +# -- Project information ----------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information + +project = "synthterrain" +copyright = "2024, United States Government, as represented by the Administrator of " +" the National Aeronautics and Space Administration." +author = "Ross Beyer" + +# -- General configuration --------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration + +extensions = ["sphinx.ext.autodoc", "sphinxcontrib.apidoc"] + +apidoc_module_dir = "../src/python/synthterrain" +apidoc_output_dir = "." +apidoc_excluded_paths = ["tests"] + +templates_path = ["_templates"] +exclude_patterns = ["_build", "Thumbs.db", ".DS_Store"] + + +# -- Options for HTML output ------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output + +html_theme = "alabaster" +html_static_path = ["_static"] diff --git a/docs/contributing.rst b/docs/contributing.rst new file mode 100644 index 0000000..e582053 --- /dev/null +++ b/docs/contributing.rst @@ -0,0 +1 @@ +.. include:: ../CONTRIBUTING.rst diff --git a/docs/index.rst b/docs/index.rst new file mode 100644 index 0000000..41e6025 --- /dev/null +++ b/docs/index.rst @@ -0,0 +1,25 @@ +.. synthterrain documentation master file, created by + sphinx-quickstart on Mon Jul 1 15:58:03 2024. + You can adapt this file completely to your liking, but it should at least + contain the root `toctree` directive. + +Welcome to synthterrain's documentation! +======================================== + +.. toctree:: + :maxdepth: 2 + :caption: Contents: + + readme + modules + contributing + authors + changelog + + +Indices and tables +================== + +* :ref:`genindex` +* :ref:`modindex` +* :ref:`search` diff --git a/docs/modules.rst b/docs/modules.rst new file mode 100644 index 0000000..bca8ad3 --- /dev/null +++ b/docs/modules.rst @@ -0,0 +1,7 @@ +synthterrain +============ + +.. toctree:: + :maxdepth: 4 + + synthterrain diff --git a/docs/readme.rst b/docs/readme.rst new file mode 100644 index 0000000..72a3355 --- /dev/null +++ b/docs/readme.rst @@ -0,0 +1 @@ +.. include:: ../README.rst diff --git a/docs/synthterrain.crater.rst b/docs/synthterrain.crater.rst new file mode 100644 index 0000000..c3f35dd --- /dev/null +++ b/docs/synthterrain.crater.rst @@ -0,0 +1,69 @@ +synthterrain.crater package +=========================== + +Submodules +---------- + +synthterrain.crater.age module +------------------------------ + +.. automodule:: synthterrain.crater.age + :members: + :undoc-members: + :show-inheritance: + +synthterrain.crater.cli module +------------------------------ + +.. automodule:: synthterrain.crater.cli + :members: + :undoc-members: + :show-inheritance: + +synthterrain.crater.cli\_convert module +--------------------------------------- + +.. automodule:: synthterrain.crater.cli_convert + :members: + :undoc-members: + :show-inheritance: + +synthterrain.crater.cli\_plot module +------------------------------------ + +.. automodule:: synthterrain.crater.cli_plot + :members: + :undoc-members: + :show-inheritance: + +synthterrain.crater.diffusion module +------------------------------------ + +.. automodule:: synthterrain.crater.diffusion + :members: + :undoc-members: + :show-inheritance: + +synthterrain.crater.functions module +------------------------------------ + +.. automodule:: synthterrain.crater.functions + :members: + :undoc-members: + :show-inheritance: + +synthterrain.crater.profile module +---------------------------------- + +.. automodule:: synthterrain.crater.profile + :members: + :undoc-members: + :show-inheritance: + +Module contents +--------------- + +.. automodule:: synthterrain.crater + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/synthterrain.rock.rst b/docs/synthterrain.rock.rst new file mode 100644 index 0000000..a60245a --- /dev/null +++ b/docs/synthterrain.rock.rst @@ -0,0 +1,29 @@ +synthterrain.rock package +========================= + +Submodules +---------- + +synthterrain.rock.cli module +---------------------------- + +.. automodule:: synthterrain.rock.cli + :members: + :undoc-members: + :show-inheritance: + +synthterrain.rock.functions module +---------------------------------- + +.. automodule:: synthterrain.rock.functions + :members: + :undoc-members: + :show-inheritance: + +Module contents +--------------- + +.. automodule:: synthterrain.rock + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/synthterrain.rst b/docs/synthterrain.rst new file mode 100644 index 0000000..e6c69ff --- /dev/null +++ b/docs/synthterrain.rst @@ -0,0 +1,38 @@ +synthterrain package +==================== + +Subpackages +----------- + +.. toctree:: + :maxdepth: 4 + + synthterrain.crater + synthterrain.rock + +Submodules +---------- + +synthterrain.cli module +----------------------- + +.. automodule:: synthterrain.cli + :members: + :undoc-members: + :show-inheritance: + +synthterrain.util module +------------------------ + +.. automodule:: synthterrain.util + :members: + :undoc-members: + :show-inheritance: + +Module contents +--------------- + +.. automodule:: synthterrain + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/synthterrain_ARC-18971-1_Corporate_CLA.pdf b/docs/synthterrain_ARC-18971-1_Corporate_CLA.pdf new file mode 100644 index 0000000..c221f0b Binary files /dev/null and b/docs/synthterrain_ARC-18971-1_Corporate_CLA.pdf differ diff --git a/docs/synthterrain_ARC-18971-1_Individual_CLA.pdf b/docs/synthterrain_ARC-18971-1_Individual_CLA.pdf new file mode 100644 index 0000000..35fc562 Binary files /dev/null and b/docs/synthterrain_ARC-18971-1_Individual_CLA.pdf differ diff --git a/environment.yml b/environment.yml index cefb2fb..ebcc101 100644 --- a/environment.yml +++ b/environment.yml @@ -13,5 +13,9 @@ dependencies: - numpy - pandas - rasterio + - scikit-image - scipy - shapely + - pip + - pip: + - opensimplex diff --git a/environment_dev.yml b/environment_dev.yml index fe8c2ec..fc4865f 100644 --- a/environment_dev.yml +++ b/environment_dev.yml @@ -3,9 +3,10 @@ channels: - conda-forge - defaults dependencies: - - bump2version - - cookiecutter + - black + - bump-my-version - flake8 + - icecream - ipympl - jupyterlab - pytest @@ -14,5 +15,7 @@ dependencies: - sphinxcontrib - sphinxcontrib-apidoc - sphinxcontrib-autoprogram - - tox - twine + - pip + - pip: + - ufmt diff --git a/pyproject.toml b/pyproject.toml index 374b58c..4f74e52 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,68 @@ [build-system] -requires = [ - "setuptools>=42", - "wheel" -] +requires = ["setuptools >= 61.0"] build-backend = "setuptools.build_meta" + +[project] +name = "synthterrain" +version = "0.1.0" +# dynamic = ["version", "dependencies"] +dynamic = ["dependencies"] +description = "The synthterrain package is software to support the creation of synthetic terrain in the solar system." +maintainers = [ + {name = "Ross Beyer", email = "Ross.A.Beyer@nasa.gov"} +] +readme = "README.rst" +requires-python = ">=3.8" + +classifiers = [ + "Programming Language :: Python :: 3.8", + "Programming Language :: Python :: 3.9", + "Intended Audience :: Science/Research", + "License :: OSI Approved :: Apache Software License", + "Development Status :: 2 - Pre-Alpha", + "Natural Language :: English", +] + +[project.scripts] +synthcraters = "synthterrain.crater.cli:main" +synthrocks = "synthterrain.rock.cli:main" +synthterrain = "synthterrain.cli:main" +synthcraterplot = "synthterrain.crater.cli_plot:main" +synthcraterconvert = "synthterrain.crater.cli_convert:main" + +[project.urls] +Repository = "https://github.com/NeoGeographyToolkit/synthterrain" + +[tools.setup.dynamic] +# version = {attr = "synthterrain.__version__"} +dependencies = {file = ["requirements.txt"]} + +[tool.setuptools] +package-dir = {"" = "src/python"} + +[tool.bumpversion] +current_version = "0.1.0" +parse = "(?P\\d+)\\.(?P\\d+)\\.(?P\\d+)(?:-(?P[a-z]+))?" +serialize = ["{major}.{minor}.{patch}-{release}", "{major}.{minor}.{patch}"] +search = "{current_version}" +replace = "{new_version}" +regex = false +ignore_missing_version = false +ignore_missing_files = false +tag = false +sign_tags = false +tag_name = "v{new_version}" +tag_message = "Bump version: {current_version} → {new_version}" +allow_dirty = false +commit = false +message = "Bump version: {current_version} → {new_version}" +commit_args = "" + +[tool.bumpversion.parts.release] +values = ["dev", "released"] +optional_value = "released" + +[[tool.bumpversion.files]] +filename = "src/python/synthterrain/__init__.py" +search = "__version__ = {current_version}" +replace = "__version__ = {new_version}" diff --git a/requirements.txt b/requirements.txt index 6c2c1d5..557a68f 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,5 +3,6 @@ numpy opensimplex pandas rasterio +scikit-image scipy shapely diff --git a/requirements_dev.txt b/requirements_dev.txt index 2a80230..360ed0d 100644 --- a/requirements_dev.txt +++ b/requirements_dev.txt @@ -1,12 +1,9 @@ -pip==19.2.3 -bump2version==0.5.11 -wheel==0.33.6 -watchdog==0.9.0 -flake8==3.7.8 -tox==3.14.0 -coverage==4.5.4 -Sphinx==1.8.5 -twine==1.14.0 - -pytest==4.6.5 - +black +coverage +pip +bump-my-version +flake8 +pytest +Sphinx +twine +ufmt diff --git a/setup.cfg b/setup.cfg index d66b532..dbe8990 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,64 +1,4 @@ -[bumpversion] -current_version = 0.1.0 -commit = True -tag = False -parse = (?P\d+)\.(?P\d+)\.(?P\d+)(\-(?P[a-z]+))? -serialize = - {major}.{minor}.{patch}-{release} - {major}.{minor}.{patch} - -[bumpversion:part:release] -optional_value = production -values = - dev - production - -[bumpversion:file:setup.py] -search = version="{current_version}" -replace = version="{new_version}" - -[bumpversion:file:src/vipersci/__init__.py] -search = __version__ = "{current_version}" -replace = __version__ = "{new_version}" - -[metadata] -name = synthterrain -version = attr:synthterrain.__version__ -author = attr:synthterrain.__version__ -author_email = rbeyer@seti.org -description = The synthterrain package is software to support the creation of synthetic terrain in the solar system. -long_description = file: README.rst, CHANGELOG.rst -url = https://github.com/NeoGeographyToolkit/synthterrain -classifiers = - Programming Language :: Python :: 3.7 - Programming Language :: Python :: 3.8 - Programming Language :: Python :: 3.9 - License :: OSI Approved :: BSD License - Operating System :: OS Independent - Development Status :: 2 - Pre-Alpha - Natural Language :: English - -[options] -include_package_data = True -package_dir = src/python -packages = find: -# python_requires = >=3.7 - -[bdist_wheel] -universal = 1 - [flake8] exclude = docs - -[aliases] -# Define setup.py command aliases here -test = pytest - -[tool:pytest] -addopts = --ignore=setup.py - -[develop] -script_dir=$base/lib/synthterrain - -[install] -install_scripts=$base/lib/synthterrain \ No newline at end of file +max_line_length = 88 +ignore = E203,E701,W503 diff --git a/setup.py b/setup.py deleted file mode 100644 index 104325e..0000000 --- a/setup.py +++ /dev/null @@ -1,48 +0,0 @@ -#!/usr/bin/env python - -"""The setup script.""" - -from setuptools import setup, find_packages - -package_name = 'synthterrain' -author = 'Ross Beyer' -author_email = 'ross.a.beyer@nasa.gov' - -requirements = [ - "setuptools" -] - -setup( - name=package_name, - version='0.1.0', - author=author, - author_email=author_email, - maintainer=author, - maintainer_email=author_email, - entry_points={ - 'console_scripts': [ - 'synthcraters=synthterrain.crater.cli:main', - 'synthrocks=synthterrain.rock.cli:main', - 'synthterrain=synthterrain.cli:main', - 'craterplot=synthterrain.crater.cli_plot:main', - 'craterconvert=synthterrain.crater.cli_convert:main', - ], - }, - data_files=[ - ('share/ament_index/resource_index/packages', - ['resource/' + package_name]), - ('share/' + package_name, ['package.xml']), - ], - install_requires=requirements, - include_package_data=True, - # package_data={ - # "vipersci": ["data/*"], - # }, - packages=find_packages( - where="src/python", - ), - test_suite="tests/python", - zip_safe=False, - package_dir={"": "src/python"}, - tests_require=["pytest"], -) diff --git a/src/python/synthterrain/__init__.py b/src/python/synthterrain/__init__.py index c5d0c76..9d202f8 100644 --- a/src/python/synthterrain/__init__.py +++ b/src/python/synthterrain/__init__.py @@ -1,5 +1,5 @@ """Top-level package for synthterrain.""" __author__ = """synthterrain Developers""" -__email__ = "rbeyer@seti.org" +__email__ = "Ross.A.Beyer@nasa.gov" __version__ = "0.1.0" diff --git a/src/python/synthterrain/cli.py b/src/python/synthterrain/cli.py index 8e6bc51..eaf31f1 100644 --- a/src/python/synthterrain/cli.py +++ b/src/python/synthterrain/cli.py @@ -1,27 +1,35 @@ -#!/usr/bin/env python # -*- coding: utf-8 -*- """Generates synthetic crater and rock populations. """ -# Copyright 2022, synthterrain developers. +# Copyright © 2024, United States Government, as represented by the +# Administrator of the National Aeronautics and Space Administration. +# All rights reserved. # -# Reuse is permitted under the terms of the license. -# The AUTHORS file and the LICENSE file are at the -# top level of this library. +# The “synthterrain” software is licensed under the Apache License, +# Version 2.0 (the "License"); you may not use this file except in +# compliance with the License. You may obtain a copy of the License +# at http://www.apache.org/licenses/LICENSE-2.0. + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +# implied. See the License for the specific language governing +# permissions and limitations under the License. import logging -from pathlib import Path import sys +from pathlib import Path import pandas as pd from shapely.geometry import box import synthterrain.crater as cr +import synthterrain.rock as rk +from synthterrain import util from synthterrain.crater import functions as cr_dist from synthterrain.crater.cli import csfd_dict -import synthterrain.rock as rk from synthterrain.rock import functions as rk_func -import synthterrain.util as util logger = logging.getLogger(__name__) @@ -37,73 +45,76 @@ def arg_parser(): nargs=4, type=float, default=[0, 1000, 1000, 0], - metavar=('MINX', 'MAXY', 'MAXX', 'MINY'), + metavar=("MINX", "MAXY", "MAXX", "MINY"), help="The coordinates of the bounding box, expressed in meters, to " - "evaluate in min-x, max-y, max-x, min-y order (which is ulx, " - "uly, lrx, lry, the GDAL pattern). " - "Default: %(default)s" + "evaluate in min-x, max-y, max-x, min-y order (which is ulx, " + "uly, lrx, lry, the GDAL pattern). " + "Default: %(default)s", ) parser.add_argument( - "-c", "--craters", + "-c", + "--craters", type=Path, help="Crater CSV or XML file of pre-existing craters. This option is usually " - "used as follows: A set of 'real' craters are identified from a target " - "area above a certain diameter (say 5 m/pixel) and given to this option. " - "Then --cr-mind and --cr_maxd are set to some range less than 5 m/pixel. " - "This generates synthetic craters in the specified range, and then uses " - "those synthetic craters in addition to the craters from --craters when " - "building the rock probability map." + "used as follows: A set of 'real' craters are identified from a target " + "area above a certain diameter (say 5 m/pixel) and given to this option. " + "Then --cr-mind and --cr_maxd are set to some range less than 5 m/pixel. " + "This generates synthetic craters in the specified range, and then uses " + "those synthetic craters in addition to the craters from --craters when " + "building the rock probability map.", ) parser.add_argument( "--csfd", default="VIPER_Env_Spec", choices=csfd_dict.keys(), help="The name of the crater size-frequency distribution to use. " - f"Options are: {', '.join(csfd_dict.keys())}. " - "Default: %(default)s" + f"Options are: {', '.join(csfd_dict.keys())}. " + "Default: %(default)s", ) parser.add_argument( "--cr_maxd", default=1000, type=float, - help="Maximum crater diameter in meters. Default: %(default)s" + help="Maximum crater diameter in meters. Default: %(default)s", ) parser.add_argument( "--cr_mind", default=1, type=float, - help="Minimum crater diameter in meters. Default: %(default)s" + help="Minimum crater diameter in meters. Default: %(default)s", ) parser.add_argument( "--rk_maxd", default=2, type=float, - help="Maximum crater diameter in meters. Default: %(default)s" + help="Maximum crater diameter in meters. Default: %(default)s", ) parser.add_argument( "--rk_mind", default=0.1, type=float, - help="Minimum crater diameter in meters. Default: %(default)s" + help="Minimum crater diameter in meters. Default: %(default)s", ) parser.add_argument( - "-p", "--plot", + "-p", + "--plot", action="store_true", help="This will cause a matplotlib windows to open with some summary " - "plots after the program has generated craters and then rocks." + "plots after the program has generated craters and then rocks.", ) parser.add_argument( - "-x", "--xml", + "-x", + "--xml", action="store_true", help="Default output is in CSV format, but if given this will result " - "in XML output that conforms to the old MATLAB code." + "in XML output that conforms to the old MATLAB code.", ) parser.add_argument( "--probability_map_gsd", type=float, default=1, help="This program builds a probability map to generate locations, and this " - "sets the ground sample distance in the units of --bbox for that map.", + "sets the ground sample distance in the units of --bbox for that map.", ) group = parser.add_mutually_exclusive_group(required=True) group.add_argument( @@ -111,18 +122,11 @@ def arg_parser(): action=util.PrintDictAction, dict=csfd_dict, help="If given, will list detailed information about each of the " - "available CSFDs and exit." - ) - group.add_argument( - "--cr_outfile", - type=Path, - help="Path to crater output file." + "available CSFDs and exit.", ) + group.add_argument("--cr_outfile", type=Path, help="Path to crater output file.") parser.add_argument( - "--rk_outfile", - type=Path, - required=True, - help="Path to crater output file." + "--rk_outfile", type=Path, required=True, help="Path to crater output file." ) return parser @@ -165,8 +169,6 @@ def main(): df, pmap, [poly.bounds[0], poly.bounds[2], poly.bounds[1], poly.bounds[3]] ) - return - if __name__ == "__main__": sys.exit(main()) diff --git a/src/python/synthterrain/crater/__init__.py b/src/python/synthterrain/crater/__init__.py index e3bfbfe..78fffd7 100644 --- a/src/python/synthterrain/crater/__init__.py +++ b/src/python/synthterrain/crater/__init__.py @@ -1,25 +1,34 @@ -#!/usr/bin/env python # -*- coding: utf-8 -*- """Generates synthetic crater populations. """ -# Copyright 2022-2023, synthterrain developers. +# Copyright © 2024, United States Government, as represented by the +# Administrator of the National Aeronautics and Space Administration. +# All rights reserved. # -# Reuse is permitted under the terms of the license. -# The AUTHORS file and the LICENSE file are at the -# top level of this library. +# The “synthterrain” software is licensed under the Apache License, +# Version 2.0 (the "License"); you may not use this file except in +# compliance with the License. You may obtain a copy of the License +# at http://www.apache.org/licenses/LICENSE-2.0. + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +# implied. See the License for the specific language governing +# permissions and limitations under the License. import logging import math -from pathlib import Path import random +from pathlib import Path -from matplotlib.collections import PatchCollection -from matplotlib.patches import Circle import matplotlib.pyplot as plt -from matplotlib.ticker import ScalarFormatter import numpy as np import pandas as pd + +from matplotlib.collections import PatchCollection +from matplotlib.patches import Circle +from matplotlib.ticker import ScalarFormatter from shapely.geometry import Point, Polygon from synthterrain.crater import functions @@ -37,9 +46,10 @@ def synthesize( by_bin=True, min_d=None, max_d=None, + return_surfaces=False, ): """Return a pandas DataFrame which contains craters and their properties - synthesized from the input parameters. + synthesized from the input parameters. """ if production_fn is None: production_fn = determine_production_function(crater_dist.a, crater_dist.b) @@ -64,15 +74,6 @@ def synthesize( f"{math.ceil(df['age'].max()):,} years." ) - # Generate depth to diameter ratio - if by_bin: - df = diffuse_d_over_D_by_bin(df, start_dd_mean="Stopar step") - else: - df["d/D"] = df.apply( - lambda crater: diffuse_d_over_D(crater["diameter"], crater["age"]), - axis=1, - ) - # Randomly generate positions within the polygon for the locations of # the craters. logger.info("Generating center positions.") @@ -82,29 +83,51 @@ def synthesize( df["x"] = xlist df["y"] = ylist + # Generate depth to diameter ratio + if by_bin: + df = diffuse_d_over_D_by_bin( + df, start_dd_mean="Stopar step", return_surfaces=return_surfaces + ) + else: + if return_surfaces: + df["surface"] = None + df["surface"].astype(object) + df["d/D", "surface"] = df.apply( + lambda crater: diffuse_d_over_D( + crater["diameter"], crater["age"], return_surface=True + ), + axis=1, + result_type="expand", + ) + else: + df["d/D"] = df.apply( + lambda crater: diffuse_d_over_D(crater["diameter"], crater["age"]), + axis=1, + ) + return df def determine_production_function(a: float, b: float): if a >= 10: return functions.NPF(a, b) - elif b <= 2.76: + if b <= 2.76: return functions.Grun(a, b) - else: - return functions.GNPF(a, b) + + return functions.GNPF(a, b) def random_points(poly: Polygon, num_points: int): """Returns two lists, the first being the x coordinates, and the second - being the y coordinates, each *num_points* long that represent - random locations within the provided *poly*. + being the y coordinates, each *num_points* long that represent + random locations within the provided *poly*. """ # We could have returned a list of shapely Point objects, but that's # not how we need the data later. min_x, min_y, max_x, max_y = poly.bounds # points = [] - x_list = list() - y_list = list() + x_list = [] + y_list = [] # while len(points) < num_points: while len(x_list) < num_points: random_point = Point( @@ -126,13 +149,11 @@ def generate_diameters(crater_dist, area, min, max): larger than *max* will be returned. """ size = crater_dist.count(area, min) - crater_dist.count(area, max) - diameters = list() + diameters = [] while len(diameters) != size: d = crater_dist.rvs(size=(size - len(diameters))) - diameters += d[ - np.logical_and(min <= d, d <= max) - ].tolist() + diameters += d[np.logical_and(min <= d, d <= max)].tolist() return np.array(diameters) @@ -177,13 +198,13 @@ def plot(df): cumulative=-1, log=True, bins=50, - histtype='stepfilled', - label="Craters" + histtype="stepfilled", + label="Craters", ) ax_csfd.set_ylabel("Count") ax_csfd.yaxis.set_major_formatter(ScalarFormatter()) ax_csfd.set_xlabel("Diameter (m)") - ax_csfd.legend(loc='best', frameon=False) + ax_csfd.legend(loc="best", frameon=False) ax_age.scatter(df["diameter"], df["age"], alpha=0.2, edgecolors="none", s=10) ax_age.set_xscale("log") @@ -199,9 +220,8 @@ def plot(df): ax_dd.set_xlabel("Diameter (m)") patches = [ - Circle((x_, y_), s_) for x_, y_, s_ in np.broadcast( - df["x"], df["y"], df["diameter"] / 2 - ) + Circle((x_, y_), s_) + for x_, y_, s_ in np.broadcast(df["x"], df["y"], df["diameter"] / 2) ] collection = PatchCollection(patches) collection.set_array(df["d/D"]) # Sets circle color to this data property. @@ -211,11 +231,9 @@ def plot(df): fig.colorbar(collection, ax=ax_location) plt.show() - return def to_file(df: pd.DataFrame, outfile: Path, xml=False): - if xml: # Write out the dataframe in the XML style of the old MATLAB # program. @@ -245,8 +263,6 @@ def to_file(df: pd.DataFrame, outfile: Path, xml=False): columns=["x", "y", "diameter", "age", "d/D"], ) - return - def from_file(infile: Path): """Load previously written crater information from disk""" diff --git a/src/python/synthterrain/crater/age.py b/src/python/synthterrain/crater/age.py index 4a9475d..98487da 100644 --- a/src/python/synthterrain/crater/age.py +++ b/src/python/synthterrain/crater/age.py @@ -1,12 +1,22 @@ -#!/usr/bin/env python # -*- coding: utf-8 -*- """ Functions for estimating crater ages. """ -# Copyright 2022-2023, United States Government as represented by the +# Copyright © 2024, United States Government, as represented by the # Administrator of the National Aeronautics and Space Administration. # All rights reserved. +# +# The “synthterrain” software is licensed under the Apache License, +# Version 2.0 (the "License"); you may not use this file except in +# compliance with the License. You may obtain a copy of the License +# at http://www.apache.org/licenses/LICENSE-2.0. + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +# implied. See the License for the specific language governing +# permissions and limitations under the License. import bisect import logging @@ -47,25 +57,31 @@ def estimate_age(diameter, dd, max_age): d/D value by attempting to match the given d/D value to the diffusion shape of a crater of the given *diameter* and the *max_age* of that crater, which could be estimated via the equilibrium_ages() function. + + This function returns estimated ages in multiples of a million years. More + precision than that is not accurate for this approach. """ fresh_dd = stopar_fresh_dd(diameter) if dd > fresh_dd: return 0 - dd_rev_list = list(reversed(diffuse_d_over_D( - diameter, - max_age, - start_dd_adjust=fresh_dd, - return_steps=True - ))) + dd_rev_list = list( + reversed( + diffuse_d_over_D( + diameter, max_age, start_dd_adjust=fresh_dd, return_steps=True + ) + ) + ) nsteps = len(dd_rev_list) age_step = bisect.bisect_left(dd_rev_list, dd) years_per_step = max_age / nsteps - return (nsteps - age_step) * years_per_step + age = (nsteps - age_step) * years_per_step + + return round(age / 1e6) * 1e6 def estimate_age_by_bin( @@ -114,14 +130,16 @@ def estimate_age_by_bin( ) df["diameter_bin"] = pd.cut( - df["diameter"], bins=bin_edges, include_lowest=True, + df["diameter"], + bins=bin_edges, + include_lowest=True, ) # df["equilibrium_age"] = equilibrium_ages(df["diameter"], pd_csfd, eq_csfd) df["age"] = 0 for i, (interval, count) in enumerate( - df["diameter_bin"].value_counts(sort=False).items() + df["diameter_bin"].value_counts(sort=False).items() ): logger.info( f"Processing bin {i + 1}/{total_bins}, interval: {interval}, count: {count}" @@ -137,25 +155,25 @@ def estimate_age_by_bin( else: age = 4.5e9 - dd_rev_list = list(reversed(diffuse_d_over_D( - interval.mid, - age, - start_dd_adjust=fresh_dd, - return_steps=True - ))) + dd_rev_list = list( + reversed( + diffuse_d_over_D( + interval.mid, age, start_dd_adjust=fresh_dd, return_steps=True + ) + ) + ) nsteps = len(dd_rev_list) years_per_step = age / nsteps def guess_age(dd): age_step = bisect.bisect_left(dd_rev_list, dd) - return int((nsteps - age_step) * years_per_step) + return round(int((nsteps - age_step) * years_per_step) / 1e6) * 1e6 df.loc[df["diameter_bin"] == interval, "age"] = df.loc[ df["diameter_bin"] == interval - ].apply( - lambda row: guess_age(row["d/D"]), - axis=1 - ) + ].apply(lambda row: guess_age(row["d/D"]), axis=1) + + df["age"] = df["age"].astype("int64") logger.info("estimate_age_by_bin complete.") diff --git a/src/python/synthterrain/crater/cli.py b/src/python/synthterrain/crater/cli.py index c9a2062..80934ef 100644 --- a/src/python/synthterrain/crater/cli.py +++ b/src/python/synthterrain/crater/cli.py @@ -1,29 +1,42 @@ -#!/usr/bin/env python # -*- coding: utf-8 -*- """Generates synthetic crater populations. """ -# Copyright 2022, synthterrain developers. +# Copyright © 2024, United States Government, as represented by the +# Administrator of the National Aeronautics and Space Administration. +# All rights reserved. # -# Reuse is permitted under the terms of the license. -# The AUTHORS file and the LICENSE file are at the -# top level of this library. +# The “synthterrain” software is licensed under the Apache License, +# Version 2.0 (the "License"); you may not use this file except in +# compliance with the License. You may obtain a copy of the License +# at http://www.apache.org/licenses/LICENSE-2.0. + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +# implied. See the License for the specific language governing +# permissions and limitations under the License. import logging -from pathlib import Path +import math import sys +from pathlib import Path +import numpy as np +import rasterio +from rasterio.transform import from_origin +from rasterio.windows import from_bounds from shapely.geometry import box -import synthterrain.crater as crater +from synthterrain import crater, util from synthterrain.crater import functions -import synthterrain.util as util +from synthterrain.crater.diffusion import make_crater_field logger = logging.getLogger(__name__) # Assemble this global variable value. -csfd_dict = dict() +csfd_dict = {} for fn in functions.equilibrium_functions: csfd_dict[fn.__name__] = fn.__doc__ @@ -44,37 +57,43 @@ def arg_parser(): nargs=4, type=float, default=[0, 1000, 1000, 0], - metavar=('MINX', 'MAXY', 'MAXX', 'MINY'), + metavar=("MINX", "MAXY", "MAXX", "MINY"), help="The coordinates of the bounding box, expressed in meters, to " - "evaluate in min-x, max-y, max-x, min-y order (which is ulx, " - "uly, lrx, lry, the GDAL pattern). " - "Default: %(default)s" + "evaluate in min-x, max-y, max-x, min-y order (which is ulx, " + "uly, lrx, lry, the GDAL pattern). " + "Default: %(default)s", ) parser.add_argument( "--csfd", default="VIPER_Env_Spec", choices=csfd_dict.keys(), help="The name of the crater size-frequency distribution to use. " - f"Options are: {', '.join(csfd_dict.keys())}. " - "Default: %(default)s" + f"Options are: {', '.join(csfd_dict.keys())}. " + "Default: %(default)s", ) parser.add_argument( "--maxd", default=1000, type=float, - help="Maximum crater diameter in meters. Default: %(default)s" + help="Maximum crater diameter in meters. Default: %(default)s", ) parser.add_argument( "--mind", default=1, type=float, - help="Minimum crater diameter in meters. Default: %(default)s" + help="Minimum crater diameter in meters. Default: %(default)s", ) parser.add_argument( - "-p", "--plot", + "-p", + "--plot", action="store_true", help="This will cause a matplotlib window to open with some summary " - "plots after the program has generated the data." + "plots after the program has generated the data.", + ) + parser.add_argument( + "-proj", + help="If -t is given this is needed to determine the proj string for the " + "output GeoTIFF. e.g. '+proj=eqc +R=1737400 +units=m'", ) parser.add_argument( "--run_individual", @@ -82,18 +101,28 @@ def arg_parser(): # given to the by_bin parameter of synthesize. action="store_false", help="If given, this will run a diffusion model for each synthetic " - "crater individually and depending on the area provided and the " - "crater range could cause this program to run for many hours as " - "it tried to calculate tens of thousands of diffusion models. " - "The default behavior is to gather the craters into diameter bins " - "and only run a few representative diffusion models to span the " - "parameter space." + "crater individually and depending on the area provided and the " + "crater range could cause this program to run for many hours as " + "it tried to calculate tens of thousands of diffusion models. " + "The default behavior is to gather the craters into diameter bins " + "and only run a few representative diffusion models to span the " + "parameter space.", ) parser.add_argument( - "-x", "--xml", + "-t", + "--terrain_gsd", + type=float, + help="If provided, will trigger creation of a terrain model based on bbox, and " + "the value given here will set the ground sample distance (GSD) of that model. " + "The terrain model output file will be the same as --outfile, but with a .tif " + "ending.", + ) + parser.add_argument( + "-x", + "--xml", action="store_true", help="Default output is in CSV format, but if given this will result " - "in XML output that conforms to the old MATLAB code." + "in XML output that conforms to the old MATLAB code.", ) group = parser.add_mutually_exclusive_group(required=True) group.add_argument( @@ -101,13 +130,9 @@ def arg_parser(): action=util.PrintDictAction, dict=csfd_dict, help="If given, will list detailed information about each of the " - "available CSFDs and exit." - ) - group.add_argument( - "-o", "--outfile", - type=Path, - help="Path to output file." + "available CSFDs and exit.", ) + group.add_argument("-o", "--outfile", type=Path, help="Path to output file.") return parser @@ -132,7 +157,8 @@ def main(): polygon=poly, by_bin=args.run_individual, min_d=args.mind, - max_d=args.maxd + max_d=args.maxd, + return_surfaces=bool(args.terrain_gsd), ) if args.plot: @@ -141,7 +167,28 @@ def main(): # Write out results. crater.to_file(df, args.outfile, args.xml) - return + if args.terrain_gsd is not None: + transform = from_origin( + poly.bounds[0], poly.bounds[3], args.terrain_gsd, args.terrain_gsd + ) + window = from_bounds(*poly.bounds, transform=transform) + + tm = make_crater_field( + df, np.zeros((math.ceil(window.height), math.ceil(window.width))), transform + ) + + with rasterio.open( + args.outfile.with_suffix(".tif"), + "w", + driver="GTiff", + height=tm.shape[0], + width=tm.shape[1], + count=1, + dtype=tm.dtype, + crs=args.proj, + transform=transform, + ) as dst: + dst.write(tm, 1) if __name__ == "__main__": diff --git a/src/python/synthterrain/crater/cli_convert.py b/src/python/synthterrain/crater/cli_convert.py index 12271c5..256a029 100644 --- a/src/python/synthterrain/crater/cli_convert.py +++ b/src/python/synthterrain/crater/cli_convert.py @@ -1,25 +1,33 @@ -#!/usr/bin/env python # -*- coding: utf-8 -*- """Converts between the crater CSV and XML formats. """ -# Copyright 2022, synthterrain developers. +# Copyright © 2024, United States Government, as represented by the +# Administrator of the National Aeronautics and Space Administration. +# All rights reserved. # -# Reuse is permitted under the terms of the license. -# The AUTHORS file and the LICENSE file are at the -# top level of this library. +# The “synthterrain” software is licensed under the Apache License, +# Version 2.0 (the "License"); you may not use this file except in +# compliance with the License. You may obtain a copy of the License +# at http://www.apache.org/licenses/LICENSE-2.0. + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +# implied. See the License for the specific language governing +# permissions and limitations under the License. import argparse import logging -from pathlib import Path import sys +from pathlib import Path import pandas as pd -import synthterrain.crater as crater -from synthterrain.crater.age import estimate_age_by_bin import synthterrain.crater.functions as crater_func -import synthterrain.util as util + +from synthterrain import crater, util +from synthterrain.crater.age import estimate_age_by_bin logger = logging.getLogger(__name__) @@ -34,31 +42,29 @@ def arg_parser(): "--estimate_ages", action="store_true", help="When given, craters in the input file with no age specified, or an age " - "of zero, will have an age estimated based on their diameter and d/D " - "ratio using the Grun/Neukum production function and the VIPER " - "Environmental Specification equilibrium crater function. Some resulting " - "craters may still yield a zero age if the d/D ratio was large relative " - "to the diameter, indicating a very fresh crater." + "of zero, will have an age estimated based on their diameter and d/D " + "ratio using the Grun/Neukum production function and the VIPER " + "Environmental Specification equilibrium crater function. Some resulting " + "craters may still yield a zero age if the d/D ratio was large relative " + "to the diameter, indicating a very fresh crater.", ) parser.add_argument( "--full_age", action="store_true", help="Ignored unless --estimate_ages is also given. When provided, it will " - "cause the diffusion calculation to run for the age of the solar system " - "instead of just the equilibrium age for each crater size. This may " - "provide improved age estimates, but could also cause longer run times. " - "Please use with caution." + "cause the diffusion calculation to run for the age of the solar system " + "instead of just the equilibrium age for each crater size. This may " + "provide improved age estimates, but could also cause longer run times. " + "Please use with caution.", ) parser.add_argument( - "infile", - type=Path, - help="A CSV or XML file produced by synthcraters." + "infile", type=Path, help="A CSV or XML file produced by synthcraters." ) parser.add_argument( "outfile", type=Path, help="The output file an XML or CSV file (whatever the opposite of the " - "first argument is)." + "first argument is).", ) return parser @@ -81,9 +87,7 @@ def main(): df["age"] = 0 # No age information from XML file format. else: - parser.error( - f"The input file {args.infile} did not end in .csv or .xml." - ) + parser.error(f"The input file {args.infile} did not end in .csv or .xml.") if args.estimate_ages: if args.full_age: @@ -103,14 +107,12 @@ def main(): else: df = estimate_age_by_bin(df, 50, pd_csfd, eq_csfd) except ValueError: - logger.error( - "The provided file has no craters with an age of zero." - ) + logger.error("The provided file has no craters with an age of zero.") return 1 crater.to_file(df, args.outfile, xml=(args.outfile.suffix.casefold() == ".xml")) - return + return None if __name__ == "__main__": diff --git a/src/python/synthterrain/crater/cli_plot.py b/src/python/synthterrain/crater/cli_plot.py index 26504a3..32f59b3 100644 --- a/src/python/synthterrain/crater/cli_plot.py +++ b/src/python/synthterrain/crater/cli_plot.py @@ -1,23 +1,30 @@ -#!/usr/bin/env python # -*- coding: utf-8 -*- """Generates plots from .csv files. """ -# Copyright 2022, synthterrain developers. +# Copyright © 2024, United States Government, as represented by the +# Administrator of the National Aeronautics and Space Administration. +# All rights reserved. # -# Reuse is permitted under the terms of the license. -# The AUTHORS file and the LICENSE file are at the -# top level of this library. +# The “synthterrain” software is licensed under the Apache License, +# Version 2.0 (the "License"); you may not use this file except in +# compliance with the License. You may obtain a copy of the License +# at http://www.apache.org/licenses/LICENSE-2.0. + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +# implied. See the License for the specific language governing +# permissions and limitations under the License. import argparse import logging -from pathlib import Path import sys +from pathlib import Path import pandas as pd -import synthterrain.crater as crater -import synthterrain.util as util +from synthterrain import crater, util logger = logging.getLogger(__name__) @@ -32,7 +39,7 @@ def arg_parser(): "csv", type=Path, help="A CSV file with a header row, and the following columns: " - "x, y, diameter, age, d/D." + "x, y, diameter, age, d/D.", ) return parser @@ -46,8 +53,6 @@ def main(): crater.plot(df) - return - if __name__ == "__main__": sys.exit(main()) diff --git a/src/python/synthterrain/crater/diffusion.py b/src/python/synthterrain/crater/diffusion.py index 5124c39..aebe902 100644 --- a/src/python/synthterrain/crater/diffusion.py +++ b/src/python/synthterrain/crater/diffusion.py @@ -1,4 +1,3 @@ -#!/usr/bin/env python # -*- coding: utf-8 -*- """ Performs diffusion of a crater shape. @@ -10,19 +9,33 @@ https://scipython.com/book2/chapter-7-matplotlib/examples/the-two-dimensional-diffusion-equation/ """ -# Copyright 2022, synthterrain developers. +# Copyright © 2024, United States Government, as represented by the +# Administrator of the National Aeronautics and Space Administration. +# All rights reserved. # -# Reuse is permitted under the terms of the license. -# The AUTHORS file and the LICENSE file are at the -# top level of this library. +# The “synthterrain” software is licensed under the Apache License, +# Version 2.0 (the "License"); you may not use this file except in +# compliance with the License. You may obtain a copy of the License +# at http://www.apache.org/licenses/LICENSE-2.0. + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +# implied. See the License for the specific language governing +# permissions and limitations under the License. import logging import math import statistics +from typing import Union import numpy as np -from numpy.polynomial import Polynomial +import numpy.typing as npt import pandas as pd +from numpy.polynomial import Polynomial +from rasterio.transform import rowcol +from rasterio.windows import get_data_window, intersect, Window +from skimage.transform import rescale from synthterrain.crater.profile import FTmod_Crater, stopar_fresh_dd @@ -100,8 +113,8 @@ def diffuse_d_over_D( start_dd_mean=0.15, start_dd_std=0.02, return_steps=False, - return_surface=False, - crater_cls=FTmod_Crater + return_surface: Union[bool, str] = False, + crater_cls=FTmod_Crater, ): """ Returns a depth to diameter ratio of a crater of *diameter* in meters @@ -132,7 +145,8 @@ def diffuse_d_over_D( will be returned, where the nature of the zeroth element is based on *return_steps* and the final element is the numpy array of elevations which represents the final diffused surface relative to a starting flat - surface of zero. + surface of zero. If "all" is given for *return_surface*, the final tuple + element will be a list of numpy arrays, one for each time step. If a different crater profile is desired, pass a subclass (not an instance) of crater.profile.Crater to *crater_cls* that takes a depth parameter, @@ -196,18 +210,21 @@ def diffuse_d_over_D( # D * dt appears in the Hill (2020) diffusion calculation, but # we have formulated dls, such that D * dt = dls - dd_for_each_step = list() + dd_for_each_step = [ + (np.max(u) - np.min(u)) / diameter, + ] + u_for_each_step = [ + u, + ] un = np.copy(u) for step in range(nsteps): un[1:-1, 1:-1] = u[1:-1, 1:-1] + dls * ( - ( - u[2:, 1:-1] - 2 * u[1:-1, 1:-1] + u[:-2, 1:-1] - ) / dx2 + ( - u[1:-1, 2:] - 2 * u[1:-1, 1:-1] + u[1:-1, :-2] - ) / dx2 + (u[2:, 1:-1] - 2 * u[1:-1, 1:-1] + u[:-2, 1:-1]) / dx2 + + (u[1:-1, 2:] - 2 * u[1:-1, 1:-1] + u[1:-1, :-2]) / dx2 ) - dd_for_each_step.append((np.max(u) - np.min(u)) / diameter) + dd_for_each_step.append((np.max(un) - np.min(un)) / diameter) u = np.copy(un) + u_for_each_step.append(u) # Final depth to diameter: if return_steps: @@ -215,10 +232,13 @@ def diffuse_d_over_D( else: dd = dd_for_each_step[-1] + if return_surface == "all": + return dd, u_for_each_step + if return_surface: - return dd, u - else: - return dd + return dd, u_for_each_step[-1] + + return dd def diffuse_d_over_D_by_bin( @@ -227,6 +247,7 @@ def diffuse_d_over_D_by_bin( domain_size=200, start_dd_mean=0.15, start_dd_std=0.02, + return_surfaces=False, ) -> pd.DataFrame: """ Returns a pandas DataFrame identical to the input *df* but with the @@ -255,12 +276,14 @@ def diffuse_d_over_D_by_bin( by the "mean" d/D ratio curve, and a standard deviation specified by the mean of the difference of the "high" and "low" diffiusion model values with the "mean" d/D model value at that time step. + + If *return_surfaces* is True, there will be an additional column, "surface", in the + returned dataframe that contains a 2D numpy array which represents the height field + of the crater in each row. """ logger.info("diffuse_d_over_D_by_bin start.") - bin_edges = np.geomspace( - df["diameter"].min(), df["diameter"].max(), num=num + 1 - ) + bin_edges = np.geomspace(df["diameter"].min(), df["diameter"].max(), num=num + 1) # logger.info(f"{df.shape[0]} craters") logger.info( f"Divided the craters into {num} diameter bins (not all bins may have " @@ -271,21 +294,19 @@ def diffuse_d_over_D_by_bin( # This is a 3-degree fit to the data from Stopar et al. (2017) # The resulting function, stopar_dD() will return d/D ratios when # given a diameter in meters. - stopar_poly = Polynomial([ - 1.23447427e-01, 1.49135061e-04, -6.16681361e-08, 7.08449143e-12 - ]) + stopar_poly = Polynomial( + [1.23447427e-01, 1.49135061e-04, -6.16681361e-08, 7.08449143e-12] + ) def start_dd(diameter): if diameter < 850: return stopar_poly(diameter) - else: - return 0.2 + return 0.2 def start_std(diameter): if diameter < 10: return start_dd_std + 0.01 - else: - return start_dd_std + return start_dd_std elif start_dd_mean == "Stopar step": # Stopar et al. (2017) define a set of graduate d/D categories @@ -308,91 +329,223 @@ def start_std(diameter): return start_dd_std else: + def start_dd(diameter): return start_dd_mean def start_std(diameter): return start_dd_std - df["diameter_bin"] = pd.cut(df["diameter"], bins=bin_edges) + df["diameter_bin"] = pd.cut(df["diameter"], bins=bin_edges, include_lowest=True) df["d/D"] = 0.0 + df["surface"] = None + df["surface"].astype(object) # Need to convert this loop to multiprocessing. for i, (interval, count) in enumerate( df["diameter_bin"].value_counts(sort=False).items() ): - logger.info( - f"Processing bin {i}/{num}, interval: {interval}, count: {count}" - ) + logger.info(f"Processing bin {i}/{num}, interval: {interval}, count: {count}") if count == 0: continue - elif 0 < count <= 3: + + if 0 < count <= 3: # Run individual models for each crater. - df.loc[df["diameter_bin"] == interval, "d/D"] = df.loc[ - df["diameter_bin"] == interval - ].apply( - lambda row: diffuse_d_over_D( - row["diameter"], - row["age"], - domain_size=domain_size, - start_dd_adjust=True, - start_dd_mean=start_dd(row["diameter"]), - start_dd_std=start_dd(row["diameter"]), - ), - axis=1 - ) + if return_surfaces: + df.loc[df["diameter_bin"] == interval, ["d/D", "surface"]] = df.loc[ + df["diameter_bin"] == interval + ].apply( + lambda row: pd.Series( + diffuse_d_over_D( + row["diameter"], + row["age"], + domain_size=domain_size, + start_dd_adjust=True, + start_dd_mean=start_dd(row["diameter"]), + start_dd_std=start_std(row["diameter"]), + return_surface=True, + ), + index=["d/D", "surface"], + ), + axis=1, + result_type="expand", + ) + else: + df.loc[df["diameter_bin"] == interval, "d/D"] = df.loc[ + df["diameter_bin"] == interval + ].apply( + lambda row: diffuse_d_over_D( + row["diameter"], + row["age"], + domain_size=domain_size, + start_dd_adjust=True, + start_dd_mean=start_dd(row["diameter"]), + start_dd_std=start_std(row["diameter"]), + ), + axis=1, + ) else: # Run three representative models for this "bin" oldest_age = df.loc[df["diameter_bin"] == interval, "age"].max() + kappa = kappa_diffusivity(interval.mid) + dls = diffusion_length_scale(interval.mid, domain_size) + start = start_dd(interval.mid) std = start_std(interval.mid) - middle_dds = diffuse_d_over_D( - interval.mid, - oldest_age, + all_args = [interval.mid, oldest_age] + all_kwargs = dict( domain_size=domain_size, - start_dd_adjust=start, - return_steps=True + return_steps=True, ) - high_dds = diffuse_d_over_D( - interval.mid, - oldest_age, - domain_size=domain_size, - start_dd_adjust=start + std, - return_steps=True - ) - low_dds = diffuse_d_over_D( - interval.mid, - oldest_age, - domain_size=domain_size, - start_dd_adjust=start - std, - return_steps=True - ) - - kappa = kappa_diffusivity(interval.mid) - dls = diffusion_length_scale(interval.mid, domain_size) - - # Defining this in-place since it really isn't needed outside - # this function. - def dd_from_rep(age): - age_step = math.floor(age * kappa / dls) - return np.random.normal( - middle_dds[age_step], - statistics.mean([ - middle_dds[age_step] - low_dds[age_step], - high_dds[age_step] - middle_dds[age_step] - ]) + mid_kwargs = all_kwargs.copy() + mid_kwargs["start_dd_adjust"] = start + + high_kwargs = all_kwargs.copy() + high_kwargs["start_dd_adjust"] = start + std + + low_kwargs = all_kwargs.copy() + low_kwargs["start_dd_adjust"] = start - std + + if return_surfaces: + mid_kwargs["return_surface"] = "all" + + middle_dds, middle_surfs = diffuse_d_over_D(*all_args, **mid_kwargs) + high_dds = diffuse_d_over_D(*all_args, **high_kwargs) + low_dds = diffuse_d_over_D(*all_args, **low_kwargs) + + # Defining this in-place since it really isn't needed outside + # this function. + def dd_surf_from_rep(age, diam): + age_step = math.floor(age * kappa / dls) + dd = np.random.normal( + middle_dds[age_step], + statistics.mean( + [ + middle_dds[age_step] - low_dds[age_step], + high_dds[age_step] - middle_dds[age_step], + ] + ), + ) + surf = middle_surfs[age_step] + surf_depth = np.max(surf) - np.min(surf) + d_mult = dd * diam / surf_depth + + return dd, surf * d_mult + + df.loc[df["diameter_bin"] == interval, ["d/D", "surface"]] = df.loc[ + df["diameter_bin"] == interval + ].apply( + lambda row: pd.Series( + dd_surf_from_rep(row["age"], row["diameter"]), + index=["d/D", "surface"], + ), + axis=1, + result_type="expand", ) - df.loc[df["diameter_bin"] == interval, "d/D"] = df.loc[ - df["diameter_bin"] == interval - ].apply( - lambda row: dd_from_rep(row["age"]), - axis=1 - ) + else: + middle_dds = diffuse_d_over_D(*all_args, **mid_kwargs) + high_dds = diffuse_d_over_D(*all_args, **high_kwargs) + low_dds = diffuse_d_over_D(*all_args, **low_kwargs) + + # Defining this in-place since it really isn't needed outside + # this function. + def dd_from_rep(age): + age_step = math.floor(age * kappa / dls) + return np.random.normal( + middle_dds[age_step], + statistics.mean( + [ + middle_dds[age_step] - low_dds[age_step], + high_dds[age_step] - middle_dds[age_step], + ] + ), + ) + + df.loc[df["diameter_bin"] == interval, "d/D"] = df.loc[ + df["diameter_bin"] == interval + ].apply(lambda row: dd_from_rep(row["age"]), axis=1) logger.info("diffuse_d_over_D_by_bin complete.") return df + + +def make_crater_field( + df, + terrain_model, + transform, +) -> npt.NDArray: + """ + Returns a 2D numpy array whose values are heights. + + The *df* should have a "surface" column which contains 2D numpy arrays, + presumably the output of diffuse_d_over_D_by_bin() with return_surfaces=True. + + The *terrain_model* argument can either be an initial 2D Numpy Array which + contains elevation values for a surface which the craters in *df* will be + applied to and diffused over, or it can be a two-element sequence that contains + the number of rows and columns that an initial flat 2D Numpy Array will be + generated from. + + """ + logger.info("make_crater_field start.") + + # # Establish initial height field: + # if len(terrain_model.shape) == 2: + # u = terrain_model + # else: + # u = np.zeros(terrain_model) + + if abs(transform.a) != abs(transform.e): + raise ValueError("The transform does not have even spacing in X and Y.") + + gsd = transform.a + + tm_window = get_data_window(terrain_model) + + for row in df.itertuples(index=False): + surf = row.surface + if surf is None: + logger.info(f"There is no surface for this row: {row}") + continue + # So the below assumes that surf is square, and that it was built with + # diffuse_d_over_D(): + surf_gsd = row.diameter * 2 / surf.shape[0] + + new_surf = rescale( + surf, surf_gsd / gsd, preserve_range=True, anti_aliasing=True + ) + + r, c = rowcol(transform, row.x, row.y) + + surf_window = Window( + c - int(new_surf.shape[1] / 2), + r - int(new_surf.shape[0] / 2), + new_surf.shape[1], + new_surf.shape[0], + ) + tm_slices, surf_slices = to_relative_slices(tm_window, surf_window) + terrain_model[tm_slices] += new_surf[surf_slices] + + return terrain_model + + +def to_relative_slices(w1: Window, w2: Window): + if not intersect(w1, w2): + raise ValueError("The two windows do not intersect.") + + w1_r1 = Window(0, 0, w1.width, w1.height) + w2_r1 = Window( + w2.col_off - w1.col_off, w2.row_off - w1.row_off, w2.width, w2.height + ) + + w1_r2 = Window( + w1.col_off - w2.col_off, w1.row_off - w2.row_off, w1.width, w1.height + ) + w2_r2 = Window(0, 0, w2.width, w2.height) + + return w1_r1.intersection(w2_r1).toslices(), w2_r2.intersection(w1_r2).toslices() diff --git a/src/python/synthterrain/crater/functions.py b/src/python/synthterrain/crater/functions.py index 92a0ff2..394ec78 100644 --- a/src/python/synthterrain/crater/functions.py +++ b/src/python/synthterrain/crater/functions.py @@ -1,17 +1,27 @@ -#!/usr/bin/env python # -*- coding: utf-8 -*- """This module contains abstract and concrete classes for representing crater size-frequency distributions as probability distributions. """ -# Copyright 2022, United States Government as represented by the +# Copyright © 2024, United States Government, as represented by the # Administrator of the National Aeronautics and Space Administration. # All rights reserved. +# +# The “synthterrain” software is licensed under the Apache License, +# Version 2.0 (the "License"); you may not use this file except in +# compliance with the License. You may obtain a copy of the License +# at http://www.apache.org/licenses/LICENSE-2.0. + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +# implied. See the License for the specific language governing +# permissions and limitations under the License. -from abc import ABC, abstractmethod import copy import logging import math +from abc import ABC, abstractmethod from numbers import Number import numpy as np @@ -27,52 +37,51 @@ class Crater_rv_continuous(ABC, rv_continuous): """Base class for crater continuous distributions. Provides - some convenience functions common to all crater distributions. + some convenience functions common to all crater distributions. - Crater distribution terminology can be reviewed in Robbins, et al. (2018, - https://doi.org/10.1111/maps.12990), which in its - terminology section states "In crater work, the CSFD can be thought - of as a scaled version of "1−CDF." + Crater distribution terminology can be reviewed in Robbins, et al. (2018, + https://doi.org/10.1111/maps.12990), which in its + terminology section states "In crater work, the CSFD can be thought + of as a scaled version of "1−CDF." - CSFD is the crater size-frequency distribution, a widely used - representation in the scientific literature. CDF is the statistical - cumulative distribution function. + CSFD is the crater size-frequency distribution, a widely used + representation in the scientific literature. CDF is the statistical + cumulative distribution function. - Both the CSFD and the CDF are functions of d (the diameter of craters), - and are related thusly, according to Robbins et al.: + Both the CSFD and the CDF are functions of d (the diameter of craters), + and are related thusly, according to Robbins et al.: - CSFD(d) ~ 1 - CDF(d) + CSFD(d) ~ 1 - CDF(d) - For any particular count of craters, the smallest crater value - measured (d_min) gives the total number of craters in the set per - unit area, CSFD(d_min). Which implies this relation + For any particular count of craters, the smallest crater value + measured (d_min) gives the total number of craters in the set per + unit area, CSFD(d_min). Which implies this relation - CDF(d) = 1 - (CSFD(d) / CSFD(d_min)) + CDF(d) = 1 - (CSFD(d) / CSFD(d_min)) - If you scale by CSFC(d_min) which is the total number of craters, - then you get a statistical CDF. + If you scale by CSFC(d_min) which is the total number of craters, + then you get a statistical CDF. - When creating concrete classes that descend from Crater_rv_continuous, - the csfd() function must be implemented. There is a default - implementation of ._cdf(), but it is advised that it be implemented - directly for speed and efficiency. It is also heavily advised - (echoing the advice from scipy.stats.rv_continuous itself) that _ppf() - also be directly implemented. + When creating concrete classes that descend from Crater_rv_continuous, + the csfd() function must be implemented. There is a default + implementation of ._cdf(), but it is advised that it be implemented + directly for speed and efficiency. It is also heavily advised + (echoing the advice from scipy.stats.rv_continuous itself) that _ppf() + also be directly implemented. - It is assumed that the units of diameters (d) are in meters, and - that the resulting CSFD(d) is in units of number per square meter. - Implementing functions should support this. + It is assumed that the units of diameters (d) are in meters, and + that the resulting CSFD(d) is in units of number per square meter. + Implementing functions should support this. """ def __init__(self, a, **kwargs): if a <= 0: raise ValueError( - "The lower bound of the support of the distribution, a, must " - "be > 0." + "The lower bound of the support of the distribution, a, must be > 0." ) - else: - kwargs["a"] = a - super().__init__(**kwargs) + + kwargs["a"] = a + super().__init__(**kwargs) @abstractmethod def csfd(self, d): @@ -109,22 +118,24 @@ def csfd(self, d): def _cdf(self, d): """Returns an array-like which is the result of applying the - cumulative density function to *d*, the input array-like of - diameters. + cumulative density function to *d*, the input array-like of + diameters. - If the crater size-frequency distribution (CSFD) is C(d) (typically - also expressed as N_cumulative(d) ), then + If the crater size-frequency distribution (CSFD) is C(d) (typically + also expressed as N_cumulative(d) ), then - cdf(d) = 1 - (C(d) / C(d_min)) + cdf(d) = 1 - (C(d) / C(d_min)) - In the context of this class, d_min is a, the lower bound of the - support of the distribution when this class is instantiated. + In the context of this class, d_min is a, the lower bound of the + support of the distribution when this class is instantiated. - As with the parent class, rv_continuous, implementers of derived - classes are strongly encouraged to override this with more - efficient implementations, also possibly implementing _ppf(). + As with the parent class, rv_continuous, implementers of derived + classes are strongly encouraged to override this with more + efficient implementations, also possibly implementing _ppf(). """ - return np.ones_like(d, dtype=np.dtype(float)) - (self.csfd(d) / self.csfd(self.a)) + return np.ones_like(d, dtype=np.dtype(float)) - ( + self.csfd(d) / self.csfd(self.a) + ) def count(self, area, diameter=None) -> int: """Returns the number of craters based on the *area* provided @@ -142,21 +153,21 @@ def count(self, area, diameter=None) -> int: def rvs(self, *args, **kwargs): """Overrides the parent rvs() function by adding an *area* - parameter, all other arguments are identical. + parameter, all other arguments are identical. - If an *area* parameter is provided, it is interpreted as the - area in square meters which has accumulated craters. + If an *area* parameter is provided, it is interpreted as the + area in square meters which has accumulated craters. - Specifying it will cause the *size* parameter (if given) to - be overridden such that + Specifying it will cause the *size* parameter (if given) to + be overridden such that - size = CSDF(d_min) * area + size = CSDF(d_min) * area - and then the parent rvs() function will be called. + and then the parent rvs() function will be called. - Since the CSDF interpreted at the minimum crater size is - the total number of craters per square meter, multiplying - it by the desired area results in the desired number of craters. + Since the CSDF interpreted at the minimum crater size is + the total number of craters per square meter, multiplying + it by the desired area results in the desired number of craters. """ if "area" in kwargs: kwargs["size"] = self.count(kwargs["area"]) @@ -166,19 +177,18 @@ def rvs(self, *args, **kwargs): class Test_Distribution(Crater_rv_continuous): - """This is testing a simple function. - """ + """This is testing a simple function.""" def csfd(self, d): """Returns the crater cumulative size frequency distribution function - such that - CSFD(d) = N_cum(d) = 29174 / d^(1.92) + such that + CSFD(d) = N_cum(d) = 29174 / d^(1.92) """ return 29174 * np.float_power(d, -1.92) def _cdf(self, d): """Override of parent function to eliminate unnecessary division - of 29174 by itself. + of 29174 by itself. """ return np.ones_like(d, dtype=np.dtype(float)) - ( np.float_power(d, -1.92) / np.float_power(self.a, -1.92) @@ -198,13 +208,17 @@ class VIPER_Env_Spec(Crater_rv_continuous): def csfd(self, d): """ - CSFD( d <= 80 ) = (29174 / d^(1.92)) / (1000^2) - CSFD( d > 80 ) = (156228 / d^(2.389)) / (1000^2) + CSFD( d <= 80 ) = (29174 / d^(1.92)) / (1000^2) + CSFD( d > 80 ) = (156228 / d^(2.389)) / (1000^2) """ if isinstance(d, Number): # Convert to numpy array, if needed. - diam = np.array([d, ]) + diam = np.array( + [ + d, + ] + ) else: diam = d c = np.empty_like(diam, dtype=np.dtype(float)) @@ -213,8 +227,8 @@ def csfd(self, d): out = c / (1000 * 1000) if isinstance(d, Number): return out.item() - else: - return out + + return out # See comment on commented out parent isfd() function. # def isfd(self, d): @@ -231,7 +245,7 @@ def csfd(self, d): def _cdf(self, d): """Override parent function to eliminate unnecessary division - by constants. + by constants. """ c = np.empty_like(d, dtype=np.dtype(float)) c[d <= 80] = np.float_power(d[d <= 80], -1.92) / np.float_power(self.a, -1.92) @@ -240,16 +254,22 @@ def _cdf(self, d): def _ppf(self, q): """Override parent function to make things faster for .rvs().""" - q80 = float(self._cdf(np.array([80, ]))) + q80 = float( + self._cdf( + np.array( + [ + 80, + ] + ) + ) + ) ones = np.ones_like(q, dtype=np.dtype(float)) p = np.empty_like(q, dtype=np.dtype(float)) p[q <= q80] = np.float_power( - (ones[q <= q80] / (ones[q <= q80] - q[q <= q80])), - (1 / 1.92) + (ones[q <= q80] / (ones[q <= q80] - q[q <= q80])), (1 / 1.92) ) p[q > q80] = np.float_power( - (ones[q > q80] / (ones[q > q80] - q[q > q80])), - (1 / 2.389) + (ones[q > q80] / (ones[q > q80] - q[q > q80])), (1 / 2.389) ) return self.a * p @@ -314,15 +334,15 @@ def csfd(self, d): class Coef_Distribution(Crater_rv_continuous): """This class instantiates a continuous crater distribution based - on a polynomial. This notation for a crater distribution is - used by Neukum et al. (2001, https://doi.org/10.1023/A:1011989004263) - and in the craterstats package. - - The coefficients generally assume that the diameter values are in - kilometers, and the math here is based on that, but only matters for - the specification of the coefficients. The diameter values passed - to csfd() are expected to be in meters, and the returned value - is number per square meter. + on a polynomial. This notation for a crater distribution is + used by Neukum et al. (2001, https://doi.org/10.1023/A:1011989004263) + and in the craterstats package. + + The coefficients generally assume that the diameter values are in + kilometers, and the math here is based on that, but only matters for + the specification of the coefficients. The diameter values passed + to csfd() are expected to be in meters, and the returned value + is number per square meter. """ def __init__(self, *args, coef=None, poly=None, **kwargs): @@ -342,19 +362,19 @@ def __init__(self, *args, coef=None, poly=None, **kwargs): def csfd(self, d): """Returns the crater cumulative size frequency distribution function - such that - CSFD(d) = N_cum(d) = 10^x / (1000 * 1000) + such that + CSFD(d) = N_cum(d) = 10^x / (1000 * 1000) - where x is the summation of j from zero to n (typically ~ 11) of - a_j * ( lg(d/1000) )^j + where x is the summation of j from zero to n (typically ~ 11) of + a_j * ( lg(d/1000) )^j - where lg() is the base 10 logarithm, and the values a_j are the - coefficients provided via the constructor. + where lg() is the base 10 logarithm, and the values a_j are the + coefficients provided via the constructor. - Since published coefficients are typically provided for diameter - values in kilometers and areas in square kilometers, the equation - for CSFD(d) is adjusted so that diameters can be provided in units - of meters, and CSFD(d) is returned in counts per square meter. + Since published coefficients are typically provided for diameter + values in kilometers and areas in square kilometers, the equation + for CSFD(d) is adjusted so that diameters can be provided in units + of meters, and CSFD(d) is returned in counts per square meter. """ # The 1000s are to take diameters in meters, convert to kilometers # for application by the polynomial, and then division by a square @@ -430,8 +450,7 @@ def csfd(self, d): def _cdf(self, d): """Override parent function to speed up.""" return np.ones_like(d, dtype=np.dtype(float)) - np.float_power( - 10, - self.poly(np.log10(d / 1000)) - self.poly(np.log10(self.a / 1000)) + 10, self.poly(np.log10(d / 1000)) - self.poly(np.log10(self.a / 1000)) ) @@ -454,7 +473,6 @@ class NPF(Coef_Distribution): """ def __init__(self, a, b, **kwargs): - if a < 10: raise ValueError( "The lower bound of the support of the distribution, a, must " @@ -470,20 +488,31 @@ def __init__(self, a, b, **kwargs): kwargs["b"] = b super().__init__( coef=[ - -3.0768, -3.557528, 0.781027, 1.021521, -0.156012, -0.444058, - 0.019977, 0.086850, -0.005874, -0.006809, 8.25e-04, 5.54e-05 + -3.076756, + # -3.0768, + -3.557528, + 0.781027, + 1.021521, + -0.156012, + -0.444058, + 0.019977, + 0.086850, + -0.005874, + -0.006809, + 8.25e-04, + 5.54e-05, ], - **kwargs + **kwargs, ) class Interp_Distribution(Crater_rv_continuous): """This class instantiates a continuous crater distribution based - on interpolation of a set of data points. + on interpolation of a set of data points. - The input arrays assume that the diameter values are in meters - and the cumulative size frequency distribution values are in - counts per square meter. + The input arrays assume that the diameter values are in meters + and the cumulative size frequency distribution values are in + counts per square meter. """ def __init__(self, *args, diameters=None, csfds=None, func=None, **kwargs): @@ -504,15 +533,14 @@ def __init__(self, *args, diameters=None, csfds=None, func=None, **kwargs): def csfd(self, d): """Returns the crater cumulative size frequency distribution function - value for *d*. + value for *d*. """ return np.float_power(10, self.func(np.log10(d))) def _cdf(self, d): """Override parent function to speed up.""" return np.ones_like(d, dtype=np.dtype(float)) - np.float_power( - 10, - self.func(np.log10(d)) - self.func(np.log10(self.a)) + 10, self.func(np.log10(d)) - self.func(np.log10(self.a)) ) @@ -579,7 +607,7 @@ def parameters(): # Constants indicated below equation A2 from Grun et al. (1985) # First element here is arbitrarily zero so that indexes match # with printed A2 equation for easier comparison. - c = (0, 4e+29, 1.5e+44, 1.1e-2, 2.2e+3, 15.) + c = (0, 4e29, 1.5e44, 1.1e-2, 2.2e3, 15.0) gamma = (0, 1.85, 3.7, -0.52, 0.306, -4.38) def a_elem(mass, i): @@ -587,12 +615,9 @@ def a_elem(mass, i): def a2(mass): # Returns flux in m^-2 s^-1 - return ( - np.float_power( - a_elem(mass, 1) + a_elem(mass, 2) + c[3], gamma[3] - ) + - np.float_power(a_elem(mass, 4) + c[5], gamma[5]) - ) + return np.float_power( + a_elem(mass, 1) + a_elem(mass, 2) + c[3], gamma[3] + ) + np.float_power(a_elem(mass, 4) + c[5], gamma[5]) # fluxes = a2(masses) * 86400.0 * 365.25 # convert to m^-2 yr^-1 # fluxes = a2(masses) * 86400.0 * 365.25 * 1e6 # convert to km^-2 yr^-1 @@ -607,10 +632,9 @@ def a2(mass): # # r = [ (3 * m) / (4 * pi * rho) ] ^(1/3) # - rho = 2.5e+6 # g/m^-3 + rho = 2.5e6 # g/m^-3 radii = np.float_power( - (3 * masses) / (4 * math.pi * rho), - 1 / 3 + (3 * masses) / (4 * math.pi * rho), 1 / 3 ) # should be radii in meters. # Now these "impactor" radii need to be converted to crater size via @@ -645,7 +669,7 @@ def hoho_diameter( mu=0.43, # ~0.4 to 0.55 K1=0.132, K2=0.26, - Kr=(1.1 * 1.3) # Kr and KrRim + Kr=(1.1 * 1.3), # Kr and KrRim ): # This function is adapted from Caleb's research code, but is based # on Holsapple (1993, @@ -655,7 +679,7 @@ def hoho_diameter( # and a discontinuity with Neukum effvelocity = velocity * math.sin(math.radians(alpha)) - densityratio = (targdensity / rho) + densityratio = targdensity / rho # impmass=((4.0*math.pi)/3.0)*impdensity*(impradius**3.0) #impactormass pi2 = (gravity * radii) / math.pow(effvelocity, 2.0) @@ -667,12 +691,9 @@ def hoho_diameter( expthree = (2.0 + mu) / 2.0 expfour = (-3.0 * mu) / (2.0 + mu) piV = K1 * np.float_power( - (pi2 * np.float_power(densityratio, expone)) + - np.float_power( - K2 * pi3 * np.float_power(densityratio, exptwo), - expthree - ), - expfour + (pi2 * np.float_power(densityratio, expone)) + + np.float_power(K2 * pi3 * np.float_power(densityratio, exptwo), expthree), + expfour, ) V = (masses * piV) / targdensity # m3 for crater rim_radius = Kr * np.float_power(V, (1 / 3)) @@ -703,8 +724,7 @@ def __init__(self, a, b, interp="extendGrun", **kwargs): self.interp = interp else: raise ValueError( - f"The interpolation method, {interp} " - f"is not one of {interp_types}." + f"The interpolation method, {interp} " f"is not one of {interp_types}." ) # Will now construct *this* as an NPF with a Grun hidden inside. @@ -729,9 +749,7 @@ def __init__(self, a, b, interp="extendGrun", **kwargs): diameters = np.append(grun_d, n_diams) fluxes = np.append(grun_f, n_fluxes) - p = Polynomial.fit( - np.log10(diameters / 1000), np.log10(fluxes * 1e6), 11 - ) + p = Polynomial.fit(np.log10(diameters / 1000), np.log10(fluxes * 1e6), 11) self.grun = Coef_Distribution(poly=p, **grun_kwargs) elif self.interp == "interp": @@ -754,7 +772,11 @@ def __init__(self, a, b, interp="extendGrun", **kwargs): def csfd(self, d): if isinstance(d, Number): # Convert to numpy array, if needed. - diam = np.array([float(d), ]) + diam = np.array( + [ + float(d), + ] + ) else: diam = d c = np.empty_like(diam, dtype=np.dtype(float)) @@ -765,9 +787,7 @@ def csfd(self, d): c[diam < 10] = self.grun.csfd(diam[diam < 10]) elif self.interp == "linear": d_interp = np.log10((self.grunstop, 10)) - c_interp = np.log10(( - self.grun.csfd(self.grunstop), super().csfd(10) - )) + c_interp = np.log10((self.grun.csfd(self.grunstop), super().csfd(10))) # cs = CubicSpline(d_interp, c_interp) f = interp1d(d_interp, c_interp) @@ -790,8 +810,8 @@ def csfd(self, d): if isinstance(d, Number): return c.item() - else: - return c + + return c # def isfd(self, d): # if isinstance(d, Number): @@ -805,7 +825,11 @@ def csfd(self, d): def _cdf(self, d): if isinstance(d, Number): # Convert to numpy array, if needed. - diam = np.array([d, ]) + diam = np.array( + [ + d, + ] + ) else: diam = d c = np.empty_like(diam, dtype=np.dtype(float)) @@ -815,9 +839,7 @@ def _cdf(self, d): c[diam < 10] = self.grun._cdf(diam[diam < 10]) elif self.interp == "linear": d_interp = np.log10((self.grunstop, 10)) - c_interp = np.log10(( - self.grun._cdf(self.grunstop), super()._cdf(10) - )) + c_interp = np.log10((self.grun._cdf(self.grunstop), super()._cdf(10))) # cs = CubicSpline(d_interp, c_interp) f = interp1d(d_interp, c_interp) @@ -835,8 +857,8 @@ def _cdf(self, d): if isinstance(d, Number): return c.item() - else: - return c + + return c class GNPF(NPF): @@ -878,12 +900,14 @@ def __init__(self, a, b, **kwargs): grun_kwargs["csfds"] = fluxes self.grun = Interp_Distribution(**grun_kwargs) - return - def csfd(self, d): if isinstance(d, Number): # Convert to numpy array, if needed. - diam = np.array([float(d), ]) + diam = np.array( + [ + float(d), + ] + ) else: diam = d c = np.empty_like(diam, dtype=np.dtype(float)) @@ -893,13 +917,17 @@ def csfd(self, d): if isinstance(d, Number): return c.item() - else: - return c + + return c def _cdf(self, d): if isinstance(d, Number): # Convert to numpy array, if needed. - diam = np.array([d, ]) + diam = np.array( + [ + d, + ] + ) else: diam = d c = np.empty_like(diam, dtype=np.dtype(float)) @@ -959,16 +987,11 @@ def __init__(self, a, b, **kwargs): # print(np.log10(fluxes)) # np.set_printoptions(threshold=False) - p = Polynomial.fit( - np.log10(diameters), np.log10(fluxes), 11 - ) + p = Polynomial.fit(np.log10(diameters), np.log10(fluxes), 11) kwargs["a"] = a kwargs["b"] = b - super().__init__( - poly=p, - **kwargs - ) + super().__init__(poly=p, **kwargs) # If new equilibrium functions are added, add them to this list to expose them diff --git a/src/python/synthterrain/crater/profile.py b/src/python/synthterrain/crater/profile.py index 9514bc9..993937b 100644 --- a/src/python/synthterrain/crater/profile.py +++ b/src/python/synthterrain/crater/profile.py @@ -3,11 +3,20 @@ https:doi.org/10.1109/TAES.2014.120282 .""" -# Copyright 2022, synthterrain developers. +# Copyright © 2024, United States Government, as represented by the +# Administrator of the National Aeronautics and Space Administration. +# All rights reserved. # -# Reuse is permitted under the terms of the license. -# The AUTHORS file and the LICENSE file are at the -# top level of this library. +# The “synthterrain” software is licensed under the Apache License, +# Version 2.0 (the "License"); you may not use this file except in +# compliance with the License. You may obtain a copy of the License +# at http://www.apache.org/licenses/LICENSE-2.0. + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +# implied. See the License for the specific language governing +# permissions and limitations under the License. # Consider attempting to implement the Mahanti et al. (2014, # https://doi.org/10.1016/j.icarus.2014.06.023) Chebyshev polynomial approach. @@ -19,7 +28,7 @@ from numpy.polynomial import Polynomial -class Crater(): +class Crater: """A base class for establishing characteristics for a crater, in order to query its elevation at particular radial distances.""" @@ -34,46 +43,45 @@ def r(self): def profile(self, r): """Implementing classes must override this function. - This function returns a numpy array of the same shape as - *r*. - This function returns the elevation of the crater profile - at the radial distance *r* where *r* is a fraction of the - crater radius. Such that *r* = 0 is at the center, *r* = 1 - is at the crater diameter, and values of *r* greater than 1 - are distances outside the crater rim. - - Values returned in the numpy array are elevation values - in the distance units of whatever units the self.diameter - parameter is in. Values of zero are considered pre-existing - surface elevations. + This function returns a numpy array of the same shape as + *r*. + This function returns the elevation of the crater profile + at the radial distance *r* where *r* is a fraction of the + crater radius. Such that *r* = 0 is at the center, *r* = 1 + is at the crater diameter, and values of *r* greater than 1 + are distances outside the crater rim. + + Values returned in the numpy array are elevation values + in the distance units of whatever units the self.diameter + parameter is in. Values of zero are considered pre-existing + surface elevations. """ raise NotImplementedError( - f"The class {self.__name__} has not implemented elevation() as " - "required." + f"The class {self.__name__} has not implemented elevation() as " "required." ) class FT_Crater(Crater): """A crater whose profile is defined by functions described in - Fassett and Thomson (2014, https://doi.org/10.1002/2014JE004698), - equation 4. + Fassett and Thomson (2014, https://doi.org/10.1002/2014JE004698), + equation 4. """ def profile(self, r, radius_fix=True): """Returns a numpy array of elevation values based in the input numpy - array of radius fraction values, such that a radius fraction value - of 1 is at the rim, less than that interior to the crater, etc. - - A ValueError will be thrown if any values in r are < 0. - - The Fassett and Thomson (2014) paper defined equations which - placed the rim point at a radius fraction of 0.98, but that - results in a crater with a smaller diameter than specifed. - If radius_fix is True (the default) the returned profile will - extend the interior slope and place the rim at radius fraction - 1.0, but this may cause a discontinuity at the rim. If you - would like a profile with the original behavior, set radius_fix - to False. + array of radius fraction values, such that a radius fraction value + of 1 is at the rim, less than that interior to the crater, etc. + + A ValueError will be thrown if any values in r are < 0. + + The Fassett and Thomson (2014) paper defined equations which + placed the rim point at a radius fraction of 0.98, but that + results in a crater with a smaller diameter than specifed. + If radius_fix is True (the default) the returned profile will + extend the interior slope and place the rim at radius fraction + 1.0, but this may cause a discontinuity at the rim. If you + would like a profile with the original behavior, set radius_fix + to False. """ if not isinstance(r, np.ndarray): @@ -82,9 +90,7 @@ def profile(self, r, radius_fix=True): out_arr = np.zeros_like(r) if np.any(r < 0): - raise ValueError( - "The radius fraction value can't be less than zero." - ) + raise ValueError("The radius fraction value can't be less than zero.") # In F&T (2014) the boundary between inner and outer was at 0.98 # which put the rim not at r=1, Caleb's subsequent code revised @@ -98,12 +104,8 @@ def profile(self, r, radius_fix=True): inner_idx = np.logical_and(0.2 < r, r <= rim) outer_idx = np.logical_and(rim < r, r <= 1.5) - inner_poly = Polynomial( - [-0.228809953, 0.227533882, 0.083116795, -0.039499407] - ) - outer_poly = Polynomial( - [0.188253307, -0.187050452, 0.01844746, 0.01505647] - ) + inner_poly = Polynomial([-0.228809953, 0.227533882, 0.083116795, -0.039499407]) + outer_poly = Polynomial([0.188253307, -0.187050452, 0.01844746, 0.01505647]) out_arr[flat_idx] = self.diameter * -0.181 out_arr[inner_idx] = self.diameter * inner_poly(r[inner_idx]) @@ -134,11 +136,7 @@ class FTmod_Crater(Crater): thus shallower, thus larger (relative) flat floors. """ - def __init__( - self, - diameter, - depth=None - ): + def __init__(self, diameter, depth=None): super().__init__(diameter) if depth is None: self.depth = stopar_fresh_dd(self.diameter) * self.diameter @@ -147,10 +145,10 @@ def __init__( def profile(self, r): """Returns a numpy array of elevation values based in the input numpy - array of radius fraction values, such that a radius fraction value - of 1 is at the rim, less than that interior to the crater, etc. + array of radius fraction values, such that a radius fraction value + of 1 is at the rim, less than that interior to the crater, etc. - A ValueError will be thrown if any values in r are < 0. + A ValueError will be thrown if any values in r are < 0. """ if not isinstance(r, np.ndarray): @@ -159,20 +157,14 @@ def profile(self, r): out_arr = np.zeros_like(r) if np.any(r < 0): - raise ValueError( - "The radius fraction value can't be less than zero." - ) + raise ValueError("The radius fraction value can't be less than zero.") inner_idx = np.logical_and(0 <= r, r <= 0.98) rim_idx = np.logical_and(0.98 < r, r <= 1.02) outer_idx = np.logical_and(1.02 < r, r <= 1.5) - inner_poly = Polynomial( - [-0.228809953, 0.227533882, 0.083116795, -0.039499407] - ) - outer_poly = Polynomial( - [0.188253307, -0.187050452, 0.01844746, 0.01505647] - ) + inner_poly = Polynomial([-0.228809953, 0.227533882, 0.083116795, -0.039499407]) + outer_poly = Polynomial([0.188253307, -0.187050452, 0.01844746, 0.01505647]) rim_hoverd = 0.036822095 @@ -188,14 +180,14 @@ def profile(self, r): class MPD_Crater(Crater): """A crater whose profile is defined by functions described in - Martin, Parkes, and Dunstan (2014, - https:doi.org/10.1109/TAES.2014.120282). The published equations - for beta and h3 result in non-realistic profiles. For this class, - the definition of beta has been adjusted so that it is a positive value - (which we think was intended). We have also replaced the published - function for h3, with a cubic that actually matches up with h2 and h4, - although the matching with h4 is imperfect, so there is likely a better - representation for h3. + Martin, Parkes, and Dunstan (2014, + https:doi.org/10.1109/TAES.2014.120282). The published equations + for beta and h3 result in non-realistic profiles. For this class, + the definition of beta has been adjusted so that it is a positive value + (which we think was intended). We have also replaced the published + function for h3, with a cubic that actually matches up with h2 and h4, + although the matching with h4 is imperfect, so there is likely a better + representation for h3. """ def __init__( @@ -205,7 +197,7 @@ def __init__( rim_height=None, emin=0, # height of the ejecta at x = D/2 pre_rim_elevation=0, - plane_elevation=0 + plane_elevation=0, ): self.h0 = self.height_naught(diameter) self.hr0 = self.height_r_naught(diameter) @@ -245,25 +237,29 @@ def profile_x(self, x: float): if -1 <= x <= alpha: return self.h1(x, self.hr, self.h, self.hr0) - elif alpha <= x <= 0: - return self.h2( - x, self.hr, self.h, self.hr0, alpha, self.tr, self.pr - ) - elif 0 <= x <= beta: + if alpha <= x <= 0: + return self.h2(x, self.hr, self.h, self.hr0, alpha, self.tr, self.pr) + + if 0 <= x <= beta: # return self.h3( return self.h3_alt( - self.diameter, self.emin, - x, self.hr, self.h, self.hr0, alpha, beta, - self.tr, self.pr - ) - elif beta <= x: - return self.h4( - x, self.diameter, self.fc(x, self.emin, self.tr, self.pr) + self.diameter, + self.emin, + x, + self.hr, + self.h, + self.hr0, + alpha, + beta, + self.tr, + self.pr, ) - else: - # Really should not be able to get here. - raise ValueError(err_msg) + if beta <= x: + return self.h4(x, self.diameter, self.fc(x, self.emin, self.tr, self.pr)) + + # Really should not be able to get here. + raise ValueError(err_msg) @staticmethod def height_naught(diameter: float): @@ -279,16 +275,13 @@ def height_r_naught(diameter: float): @staticmethod def h1(x: float, hr: float, h: float, hr_naught: float): """Eqn 4 in Martin, Parkes, and Dunstan.""" - h_ = (hr_naught - hr + h) + h_ = hr_naught - hr + h return (h_ * math.pow(x, 2)) + (2 * h_ * x) + hr_naught @staticmethod - def h2( - x: float, hr: float, h: float, hr_naught: float, alpha: float, - tr=0, pr=0 - ): + def h2(x: float, hr: float, h: float, hr_naught: float, alpha: float, tr=0, pr=0): """Eqn 5 in Martin, Parkes, and Dunstan.""" - h_ = (hr_naught - hr + h) + h_ = hr_naught - hr + h return ((h_ * (alpha + 1)) / alpha) * math.pow(x, 2) + hr + tr - pr @staticmethod @@ -299,20 +292,33 @@ def alpha(hr: float, h: float, hr_naught: float, tr=0, pr=0): @staticmethod def h3( - x: float, hr: float, h: float, hr_naught: float, - alpha: float, beta: float, tr=0, pr=0 + x: float, + hr: float, + h: float, + hr_naught: float, + alpha: float, + beta: float, + tr=0, + pr=0, ): """Eqn 7 in Martin, Parkes, and Dunstan.""" - h_ = (hr_naught - hr + h) + h_ = hr_naught - hr + h t1 = -1 * ((2 * h_) / (3 * math.pow(beta, 2))) * math.pow(x, 3) t2 = (h_ + ((2 * h_) / beta)) * math.pow(x, 2) return t1 + t2 + hr + tr - pr @staticmethod def h3_alt( - diameter, emin, - x: float, hr: float, h: float, hr_naught: float, - alpha: float, beta: float, tr=0, pr=0 + diameter, + emin, + x: float, + hr: float, + h: float, + hr_naught: float, + alpha: float, + beta: float, + tr=0, + pr=0, ): """Improved cubic form.""" # ax^3 + bx ^ 2 + cx + d = elevation @@ -321,9 +327,7 @@ def h3_alt( # implies that c = 0. # The inflection point should be where this function meets up # with h4, so that means that the inflection point is at x = beta - h4_at_beta = MPD_Crater.h4( - beta, diameter, MPD_Crater.fc(beta, emin, tr, pr) - ) + h4_at_beta = MPD_Crater.h4(beta, diameter, MPD_Crater.fc(beta, emin, tr, pr)) a = (hr - h4_at_beta) / (2 * math.pow(beta, 3)) b = -3 * a * beta cubic = Polynomial([hr, 0, b, a]) @@ -332,7 +336,7 @@ def h3_alt( @staticmethod def beta(hr: float, h: float, hr_naught: float, tr=0, pr=0): """Eqn 8 in Martin, Parkes, and Dunstan.""" - h_ = (hr_naught - hr + h) + h_ = hr_naught - hr + h # This changes the order of hr_naught and hr from the # paper, as this ensures that this term will be positive. return (3 * (hr_naught - hr + tr - pr)) / (2 * h_) diff --git a/src/python/synthterrain/rock/__init__.py b/src/python/synthterrain/rock/__init__.py index 6a7cc83..1868689 100644 --- a/src/python/synthterrain/rock/__init__.py +++ b/src/python/synthterrain/rock/__init__.py @@ -1,32 +1,40 @@ -#!/usr/bin/env python # -*- coding: utf-8 -*- """Generates synthetic rock populations. """ -# Copyright 2022, synthterrain developers. +# Copyright © 2024, United States Government, as represented by the +# Administrator of the National Aeronautics and Space Administration. +# All rights reserved. # -# Reuse is permitted under the terms of the license. -# The AUTHORS file and the LICENSE file are at the -# top level of this library. +# The “synthterrain” software is licensed under the Apache License, +# Version 2.0 (the "License"); you may not use this file except in +# compliance with the License. You may obtain a copy of the License +# at http://www.apache.org/licenses/LICENSE-2.0. + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +# implied. See the License for the specific language governing +# permissions and limitations under the License. import logging import math -from pathlib import Path import time +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np import opensimplex +import pandas as pd from matplotlib.collections import PatchCollection from matplotlib.patches import Circle -import matplotlib.pyplot as plt from matplotlib.ticker import ScalarFormatter -import numpy as np from rasterio.features import geometry_mask from rasterio.transform import from_origin, rowcol, xy -from rasterio.windows import Window, from_bounds, intersect, shape -import pandas as pd +from rasterio.windows import from_bounds, intersect, shape, Window from shapely.geometry import Polygon - from synthterrain.crater import generate_diameters from synthterrain.crater.functions import Crater_rv_continuous @@ -346,7 +354,6 @@ def plot(df, pmap=None, extent=None): def to_file(df: pd.DataFrame, outfile: Path, xml=False): - if xml: # Write out the dataframe in the XML style of the old MATLAB # program. diff --git a/src/python/synthterrain/rock/cli.py b/src/python/synthterrain/rock/cli.py index 7ea2606..fc80b1a 100644 --- a/src/python/synthterrain/rock/cli.py +++ b/src/python/synthterrain/rock/cli.py @@ -1,24 +1,30 @@ -#!/usr/bin/env python # -*- coding: utf-8 -*- """Generates synthetic rock populations. """ -# Copyright 2022, synthterrain developers. +# Copyright © 2024, United States Government, as represented by the +# Administrator of the National Aeronautics and Space Administration. +# All rights reserved. # -# Reuse is permitted under the terms of the license. -# The AUTHORS file and the LICENSE file are at the -# top level of this library. +# The “synthterrain” software is licensed under the Apache License, +# Version 2.0 (the "License"); you may not use this file except in +# compliance with the License. You may obtain a copy of the License +# at http://www.apache.org/licenses/LICENSE-2.0. + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +# implied. See the License for the specific language governing +# permissions and limitations under the License. import logging -from pathlib import Path import sys +from pathlib import Path from shapely.geometry import box -import synthterrain.crater as crater -import synthterrain.rock as rock +from synthterrain import crater, rock, util from synthterrain.rock import functions -import synthterrain.util as util logger = logging.getLogger(__name__) @@ -121,8 +127,6 @@ def main(): # Write out results. rock.to_file(df, args.outfile, args.xml) - return - if __name__ == "__main__": sys.exit(main()) diff --git a/src/python/synthterrain/rock/functions.py b/src/python/synthterrain/rock/functions.py index e4dbf60..1d7c82a 100644 --- a/src/python/synthterrain/rock/functions.py +++ b/src/python/synthterrain/rock/functions.py @@ -1,14 +1,22 @@ -#!/usr/bin/env python # -*- coding: utf-8 -*- """This module contains abstract and concrete classes for representing rock size-frequency distributions as probability distributions. """ -# Copyright 2022, synthterrain developers. +# Copyright © 2024, United States Government, as represented by the +# Administrator of the National Aeronautics and Space Administration. +# All rights reserved. # -# Reuse is permitted under the terms of the license. -# The AUTHORS file and the LICENSE file are at the -# top level of this library. +# The “synthterrain” software is licensed under the Apache License, +# Version 2.0 (the "License"); you may not use this file except in +# compliance with the License. You may obtain a copy of the License +# at http://www.apache.org/licenses/LICENSE-2.0. + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +# implied. See the License for the specific language governing +# permissions and limitations under the License. import logging diff --git a/src/python/synthterrain/util.py b/src/python/synthterrain/util.py index af7c322..ac968ff 100644 --- a/src/python/synthterrain/util.py +++ b/src/python/synthterrain/util.py @@ -1,24 +1,32 @@ -#!/usr/bin/env python """This module contains viss utility functions.""" -# Copyright 2022, synthterrain developers. +# Copyright © 2024, United States Government, as represented by the +# Administrator of the National Aeronautics and Space Administration. +# All rights reserved. # -# Reuse is permitted under the terms of the license. -# The AUTHORS file and the LICENSE file are at the -# top level of this library. +# The “synthterrain” software is licensed under the Apache License, +# Version 2.0 (the "License"); you may not use this file except in +# compliance with the License. You may obtain a copy of the License +# at http://www.apache.org/licenses/LICENSE-2.0. + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +# implied. See the License for the specific language governing +# permissions and limitations under the License. import argparse import logging -import textwrap import sys +import textwrap import synthterrain class FileArgumentParser(argparse.ArgumentParser): """This argument parser sets the fromfile_prefix_chars to the - at-symbol (@), treats lines that begin with the octothorpe (#) - as comments, and allows multiple argument elements per line. + at-symbol (@), treats lines that begin with the octothorpe (#) + as comments, and allows multiple argument elements per line. """ def __init__(self, *args, **kwargs): @@ -43,19 +51,19 @@ def __init__(self, *args, **kwargs): def convert_arg_line_to_args(self, arg_line): if arg_line.startswith("#"): - return list() - else: - return arg_line.split() + return [] + + return arg_line.split() class PrintDictAction(argparse.Action): """A custom action that interrupts argument processing, prints - the contents of the *dict* argument, and then exits the - program. + the contents of the *dict* argument, and then exits the + program. - It may need to be placed in a mutually exclusive argument - group (see argparse documentation) with any required arguments - that your program should have. + It may need to be placed in a mutually exclusive argument + group (see argparse documentation) with any required arguments + that your program should have. """ def __init__(self, *args, dict=None, **kwargs): @@ -82,20 +90,20 @@ def parent_parser() -> argparse.ArgumentParser: "--verbose", action="count", default=0, - help="Displays additional information." + help="Displays additional information.", ) parent.add_argument( - '--version', - action='version', + "--version", + action="version", version=f"synthterrain Software version {synthterrain.__version__}", - help="Show library version number." + help="Show library version number.", ) return parent def set_logger(verblvl=None) -> None: """Sets the log level and configuration for applications.""" - logger = logging.getLogger(__name__.split(".")[0]) + logger = logging.getLogger(__name__.split(".", maxsplit=1)[0]) lvl_dict = {0: logging.WARNING, 1: logging.INFO, 2: logging.DEBUG} if verblvl in lvl_dict: lvl = lvl_dict[verblvl] @@ -114,5 +122,3 @@ def set_logger(verblvl=None) -> None: ch.setFormatter(formatter) logger.addHandler(ch) - - return diff --git a/tests/python/crater/test_age.py b/tests/python/crater/test_age.py index 0311744..9d3c754 100644 --- a/tests/python/crater/test_age.py +++ b/tests/python/crater/test_age.py @@ -1,9 +1,19 @@ -#!/usr/bin/env python """This module has tests for the synthterrain crater age functions.""" -# Copyright 2023, United States Government as represented by the +# Copyright © 2024, United States Government, as represented by the # Administrator of the National Aeronautics and Space Administration. # All rights reserved. +# +# The “synthterrain” software is licensed under the Apache License, +# Version 2.0 (the "License"); you may not use this file except in +# compliance with the License. You may obtain a copy of the License +# at http://www.apache.org/licenses/LICENSE-2.0. + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +# implied. See the License for the specific language governing +# permissions and limitations under the License. import unittest @@ -15,7 +25,6 @@ class Test_Ages(unittest.TestCase): - def test_equilibrium_age(self): diameters = np.array([1, 2, 3, 4, 5, 10, 20, 50, 100]) a = 1 @@ -24,29 +33,75 @@ def test_equilibrium_age(self): eq_func = fns.VIPER_Env_Spec(a=a, b=b) eq_ages = age.equilibrium_age(diameters, pd_func.csfd, eq_func.csfd) - np.testing.assert_allclose(eq_ages, np.array([ - 1.35931206e+06, 6.91870352e+06, 2.21362552e+07, 2.89919160e+07, - 3.57407099e+07, 6.55084979e+07, 1.38188561e+08, 4.22354647e+08, - 7.26330624e+08 - ])) + np.testing.assert_allclose( + eq_ages, + np.array( + [ + 1.359312e06, + 6.918704e06, + 2.213654e07, + 2.899164e07, + 3.573974e07, + 6.550186e07, + 1.381746e08, + 4.223119e08, + 7.262570e08, + ] + ), + rtol=1e-6, + ) def test_estimate_age(self): a = age.estimate_age(10, 0.09, 5e7) - self.assertAlmostEqual(a, 25322581, places=0) + self.assertAlmostEqual(a, 25000000, places=0) def test_estimate_age_by_bin(self): pd_func = fns.GNPF(a=1, b=1000) eq_func = fns.VIPER_Env_Spec(a=1, b=1000) - df = pd.DataFrame(data={ - 'diameter': [ - 1., 1., 2., 2., 3., 3., 4., 4., 5., 5., 10., 10., - 20., 20., 50., 50., 100., 100. - ], - 'd/D': [ - 0.05, 0.1, 0.05, 0.1, 0.05, 0.1, 0.05, 0.1, 0.05, 0.1, 0.05, 0.1, - 0.05, 0.1, 0.05, 0.12, 0.05, 0.13 - ]} + df = pd.DataFrame( + data={ + "diameter": [ + 1.0, + 1.0, + 2.0, + 2.0, + 3.0, + 3.0, + 4.0, + 4.0, + 5.0, + 5.0, + 10.0, + 10.0, + 20.0, + 20.0, + 50.0, + 50.0, + 100.0, + 100.0, + ], + "d/D": [ + 0.05, + 0.1, + 0.05, + 0.1, + 0.05, + 0.1, + 0.05, + 0.1, + 0.05, + 0.1, + 0.05, + 0.1, + 0.05, + 0.1, + 0.05, + 0.12, + 0.05, + 0.13, + ], + } ) df_out = age.estimate_age_by_bin( df, @@ -55,24 +110,44 @@ def test_estimate_age_by_bin(self): eq_func.csfd, ) - age_series = pd.Series([ - 1511054, 3538, 5768173, 12861, 12054097, 26876, 25171337, 56123, 35854662, - 0, 63160550, 294454, 137994962, 12513397, 424194730, 68788334, 702317463, - 8072614], name="age") + age_series = pd.Series( + [ + 2000000, + 0, + 6000000, + 0, + 12000000, + 0, + 25000000, + 0, + 36000000, + 0, + 63000000, + 0, + 138000000, + 12000000, + 424000000, + 68000000, + 702000000, + 8000000, + ], + name="age", + ) pd.testing.assert_series_equal(age_series, df_out["age"]) - - df2 = pd.DataFrame(data={ - 'diameter': [100., 100., 100., 100.], - 'd/D': [0.01, 0.06, 0.10, 0.17] - }) + df2 = pd.DataFrame( + data={ + "diameter": [100.0, 100.0, 100.0, 100.0], + "d/D": [0.01, 0.06, 0.10, 0.17], + } + ) df2_out = age.estimate_age_by_bin( df2, 50, # With only one diameter, num is irrelevant ) - age2_series = pd.Series([4500000000, 4500000000, 1388888888, 0], name="age") + age2_series = pd.Series([4500000000, 4500000000, 1388000000, 0], name="age") - pd.testing.assert_series_equal(age2_series, df2_out["age"]) \ No newline at end of file + pd.testing.assert_series_equal(age2_series, df2_out["age"]) diff --git a/tests/python/crater/test_functions.py b/tests/python/crater/test_functions.py index f7b9522..a16297c 100644 --- a/tests/python/crater/test_functions.py +++ b/tests/python/crater/test_functions.py @@ -1,27 +1,33 @@ -#!/usr/bin/env python """This module has tests for the synthterrain crater distribution functions.""" -# Copyright 2022, synthterrain developers. +# Copyright © 2024, United States Government, as represented by the +# Administrator of the National Aeronautics and Space Administration. +# All rights reserved. # -# Reuse is permitted under the terms of the license. -# The AUTHORS file and the LICENSE file are at the -# top level of this library. +# The “synthterrain” software is licensed under the Apache License, +# Version 2.0 (the "License"); you may not use this file except in +# compliance with the License. You may obtain a copy of the License +# at http://www.apache.org/licenses/LICENSE-2.0. + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +# implied. See the License for the specific language governing +# permissions and limitations under the License. import unittest import numpy as np from numpy.polynomial import Polynomial -import synthterrain.crater.functions as fns +import synthterrain.crater.functions as fns # usort: skip class Test_Crater_rv_continuous(unittest.TestCase): - def test_abstract(self): self.assertRaises(TypeError, fns.Crater_rv_continuous) class Test_Dist(fns.Crater_rv_continuous): - def csfd(self, d): return np.power(d, -2.0) @@ -35,10 +41,27 @@ def csfd(self, d): def test_VIPER_Env_Spec(self): rv = fns.VIPER_Env_Spec(a=1, b=100) self.assertAlmostEqual(rv.csfd(10), 0.0003507) - np.testing.assert_allclose(rv._cdf(np.array([10,])), np.array([0.98797736])) + np.testing.assert_allclose( + rv._cdf( + np.array( + [ + 10, + ] + ) + ), + np.array([0.98797736]), + ) np.testing.assert_allclose( - rv._ppf(np.array([0.5, 0.99,])), np.array([1.43478377, 11.00694171]) + rv._ppf( + np.array( + [ + 0.5, + 0.99, + ] + ) + ), + np.array([1.43478377, 11.00694171]), ) def test_Trask(self): @@ -48,16 +71,44 @@ def test_Trask(self): def test_Coef_Distribution(self): self.assertRaises(ValueError, fns.Coef_Distribution, a=10, b=300000) - rv = fns.Coef_Distribution(a=10, b=300000, poly=Polynomial( - [ - -3.0768, -3.557528, 0.781027, 1.021521, -0.156012, -0.444058, - 0.019977, 0.086850, -0.005874, -0.006809, 8.25e-04, 5.54e-05 - ])) + rv = fns.Coef_Distribution( + a=10, + b=300000, + poly=Polynomial( + [ + -3.0768, + -3.557528, + 0.781027, + 1.021521, + -0.156012, + -0.444058, + 0.019977, + 0.086850, + -0.005874, + -0.006809, + 8.25e-04, + 5.54e-05, + ] + ), + ) self.assertEqual(rv.csfd(10), 0.003796582136635746) - self.assertEqual(rv._cdf(np.array([10, ])), np.array([0,])) + self.assertEqual( + rv._cdf( + np.array( + [ + 10, + ] + ) + ), + np.array( + [ + 0, + ] + ), + ) def test_NPF(self): self.assertRaises(ValueError, fns.NPF, a=10, b=300001) self.assertRaises(ValueError, fns.NPF, a=9, b=300) rv = fns.NPF(a=10, b=300000) - self.assertEqual(rv.csfd(10), 0.003796582136635746) + self.assertEqual(rv.csfd(10), 0.0037969668020723783) diff --git a/tests/python/crater/test_init.py b/tests/python/crater/test_init.py index 89b30a6..715b68c 100644 --- a/tests/python/crater/test_init.py +++ b/tests/python/crater/test_init.py @@ -1,22 +1,29 @@ -#!/usr/bin/env python """This module has tests for the synthterrain crater init functions.""" -# Copyright 2022, synthterrain developers. +# Copyright © 2024, United States Government, as represented by the +# Administrator of the National Aeronautics and Space Administration. +# All rights reserved. # -# Reuse is permitted under the terms of the license. -# The AUTHORS file and the LICENSE file are at the -# top level of this library. +# The “synthterrain” software is licensed under the Apache License, +# Version 2.0 (the "License"); you may not use this file except in +# compliance with the License. You may obtain a copy of the License +# at http://www.apache.org/licenses/LICENSE-2.0. + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +# implied. See the License for the specific language governing +# permissions and limitations under the License. import unittest -from shapely.geometry import Point, Polygon +from shapely.geometry import Polygon -import synthterrain.crater as cr +import synthterrain.crater as cr # usort: skip import synthterrain.crater.functions as fns class Test_Init(unittest.TestCase): - def test_generate_diameters(self): min_d = 10 max_d = 11 @@ -31,8 +38,6 @@ def test_generate_diameters(self): self.assertEqual(0, d[d > max_d].size) def test_random_points(self): - poly = Polygon(((0, 0), (1, 0), (0, 1), (0, 0))) xs, ys = cr.random_points(poly, 5) - diff --git a/tests/python/test_util.py b/tests/python/test_util.py index 5ea1c71..adbacd9 100644 --- a/tests/python/test_util.py +++ b/tests/python/test_util.py @@ -1,12 +1,20 @@ -#!/usr/bin/env python # -*- coding: utf-8 -*- """This module has tests for the util module.""" -# Copyright 2022, synthterrain developers. +# Copyright © 2024, United States Government, as represented by the +# Administrator of the National Aeronautics and Space Administration. +# All rights reserved. # -# Reuse is permitted under the terms of the license. -# The AUTHORS file and the LICENSE file are at the -# top level of this library. +# The “synthterrain” software is licensed under the Apache License, +# Version 2.0 (the "License"); you may not use this file except in +# compliance with the License. You may obtain a copy of the License +# at http://www.apache.org/licenses/LICENSE-2.0. + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +# implied. See the License for the specific language governing +# permissions and limitations under the License. import argparse import logging @@ -17,7 +25,6 @@ class TestUtil(unittest.TestCase): - def test_FileArgumentParser(self): p = util.FileArgumentParser() @@ -41,12 +48,7 @@ def test_PrintDictAction(self, m_print, m_exit): a("dummy", "dummy", "dummy") self.assertEqual( m_print.call_args_list, - [ - call("a"), - call(" a value"), - call("b"), - call(" b value") - ] + [call("a"), call(" a value"), call("b"), call(" b value")], ) m_exit.assert_called_once()