diff --git a/.github/dependabot.yml b/.github/dependabot.yml new file mode 100644 index 00000000..9c39a860 --- /dev/null +++ b/.github/dependabot.yml @@ -0,0 +1,16 @@ +# To get started with Dependabot version updates, you'll need to specify which +# package ecosystems to update and where the package manifests are located. +# Please see the documentation for all configuration options: +# https://docs.github.com/github/administering-a-repository/configuration-options-for-dependency-updates + +version: 2 +updates: + - package-ecosystem: "github-actions" # See documentation for possible values + directory: ".github/workflows" # Location of package manifests + schedule: + interval: "monthly" + groups: + actions: + patterns: + - "*" + diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 12215c1f..550eb3b0 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -8,7 +8,7 @@ on: jobs: build: - uses: OpenAstronomy/github-actions-workflows/.github/workflows/publish_pure_python.yml@v1 + uses: OpenAstronomy/github-actions-workflows/.github/workflows/publish_pure_python.yml@d68193b68216da64eafaa618f53c59f5d271c56e # v1.14.0 with: upload_to_pypi: ${{ (github.event_name == 'release') && (github.event.action == 'released') }} secrets: diff --git a/.github/workflows/changelog.yml b/.github/workflows/changelog.yml index 6d0be9fb..ffd65f44 100644 --- a/.github/workflows/changelog.yml +++ b/.github/workflows/changelog.yml @@ -9,7 +9,7 @@ jobs: name: Verify that a changelog entry exists for this pull request runs-on: ubuntu-latest steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@11bd71901bbe5b1630ceea73d27597364c9af683 # v4.2.2 with: submodules: true - run: grep -P '\[[^\]]*#${{github.event.number}}[,\]]' CHANGES.rst diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index f87919f5..766823c9 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -1,4 +1,4 @@ -name: 'CI' +name: test on: push: @@ -7,6 +7,7 @@ on: tags: - '*' pull_request: + workflow_dispatch: schedule: # Weekly Monday 9AM build # * is a special character in YAML so you have to quote this string @@ -18,75 +19,20 @@ concurrency: jobs: check: - name: - uses: OpenAstronomy/github-actions-workflows/.github/workflows/tox.yml@v1 + uses: OpenAstronomy/github-actions-workflows/.github/workflows/tox.yml@d68193b68216da64eafaa618f53c59f5d271c56e # v1.14.0 with: envs: | - linux: check-style - linux: check-security test: - uses: OpenAstronomy/github-actions-workflows/.github/workflows/tox.yml@main + uses: OpenAstronomy/github-actions-workflows/.github/workflows/tox.yml@d68193b68216da64eafaa618f53c59f5d271c56e # v1.14.0 with: envs: | - - linux: test - python-version: 3.9 - - linux: test-numpy125 - python-version: 3.9 - - linux: test-numpy123 - python-version: 3.9 - - linux: test - python-version: 3.10 - - linux: test-numpy125 - python-version: 3.10 - - linux: test-numpy123 - python-version: 3.10 - - linux: test - python-version: 3.11 - pytest-results-summary: true - - macos: test - python-version: 3.11 - pytest-results-summary: true - - linux: test-dev - - linux: test-pyargs - python-version: 3.11 - - linux: test-cov - python-version: 3.11 + - linux: py310-oldestdeps + - linux: py311-cov coverage: codecov pytest-results-summary: true - crds: - name: retrieve current CRDS context - runs-on: ubuntu-latest - env: - OBSERVATORY: jwst - CRDS_PATH: /tmp/crds_cache - CRDS_SERVER_URL: https://jwst-crds.stsci.edu - steps: - - id: context - run: > - echo "pmap=$( - curl -s -X POST -d '{"jsonrpc": "1.0", "method": "get_default_context", "params": ["${{ env.OBSERVATORY }}"], "id": 1}' ${{ env.CRDS_SERVER_URL }}/json/ | - python -c "import sys, json; print(json.load(sys.stdin)['result'])" - )" >> $GITHUB_OUTPUT - # Get default CRDS_CONTEXT without installing crds client - # See https://hst-crds.stsci.edu/static/users_guide/web_services.html#generic-request - - id: path - run: echo "path=${{ env.CRDS_PATH }}" >> $GITHUB_OUTPUT - - id: server - run: echo "url=${{ env.CRDS_SERVER_URL }}" >> $GITHUB_OUTPUT - outputs: - context: ${{ steps.context.outputs.pmap }} - path: ${{ steps.path.outputs.path }} - server: ${{ steps.server.outputs.url }} - test_downstream: - uses: OpenAstronomy/github-actions-workflows/.github/workflows/tox.yml@main - needs: [ crds ] - with: - setenv: | - CRDS_PATH: ${{ needs.crds.outputs.path }} - CRDS_SERVER_URL: ${{ needs.crds.outputs.server }} - CRDS_CLIENT_RETRY_COUNT: 3 - CRDS_CLIENT_RETRY_DELAY_SECONDS: 20 - cache-path: ${{ needs.crds.outputs.path }} - cache-key: crds-${{ needs.crds.outputs.context }} - envs: | - - linux: py310-test-jwst-cov-xdist + - macos: py311 + pytest-results-summary: true + - linux: py312 + - linux: py3-dev diff --git a/.github/workflows/contexts.yml b/.github/workflows/contexts.yml new file mode 100644 index 00000000..ec203905 --- /dev/null +++ b/.github/workflows/contexts.yml @@ -0,0 +1,39 @@ +name: contexts + +on: + workflow_call: + outputs: + jwst: + value: ${{ jobs.contexts.outputs.jwst }} + roman: + value: ${{ jobs.contexts.outputs.roman }} + workflow_dispatch: + +jobs: + contexts: + name: retrieve latest CRDS contexts + runs-on: ubuntu-latest + outputs: + jwst: ${{ steps.jwst_crds_context.outputs.pmap }} + roman: ${{ steps.roman_crds_context.outputs.pmap }} + steps: + - id: jwst_crds_context + env: + OBSERVATORY: jwst + CRDS_SERVER_URL: https://jwst-crds.stsci.edu + run: > + echo "pmap=$( + curl -s -X POST -d '{"jsonrpc": "1.0", "method": "get_default_context", "params": ["${{ env.OBSERVATORY }}", null], "id": 1}' ${{ env.CRDS_SERVER_URL }}/json/ --retry 8 --connect-timeout 10 | + python -c "import sys, json; print(json.load(sys.stdin)['result'])" + )" >> $GITHUB_OUTPUT + - run: if [[ ! -z "${{ steps.jwst_crds_context.outputs.pmap }}" ]]; then echo ${{ steps.jwst_crds_context.outputs.pmap }}; else exit 1; fi + - id: roman_crds_context + env: + OBSERVATORY: roman + CRDS_SERVER_URL: https://roman-crds.stsci.edu + run: > + echo "pmap=$( + curl -s -X POST -d '{"jsonrpc": "1.0", "method": "get_default_context", "params": ["${{ env.OBSERVATORY }}", null], "id": 1}' ${{ env.CRDS_SERVER_URL }}/json/ --retry 8 | + python -c "import sys, json; print(json.load(sys.stdin)['result'])" + )" >> $GITHUB_OUTPUT + - run: if [[ ! -z "${{ steps.roman_crds_context.outputs.pmap }}" ]]; then echo ${{ steps.roman_crds_context.outputs.pmap }}; else exit 1; fi diff --git a/.github/workflows/downstream.yml b/.github/workflows/downstream.yml new file mode 100644 index 00000000..434469f1 --- /dev/null +++ b/.github/workflows/downstream.yml @@ -0,0 +1,109 @@ +name: test downstream packages + +on: + workflow_dispatch: + schedule: + # Weekly Monday 9AM build + # * is a special character in YAML so you have to quote this string + - cron: '0 9 * * 1' + pull_request: + # We also want this workflow triggered if the `Downstream CI` label is + # added or present when PR is updated + types: + - synchronize + - labeled + push: + tags: + - '*' + +# Only cancel in-progress jobs or runs for the current workflow +# This cancels the already triggered workflows for a specific PR without canceling +# other instances of this workflow (other PRs, scheduled triggers, etc) when something +# within that PR re-triggers this CI +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true + +jobs: + environment: + runs-on: ubuntu-latest + steps: + - id: data_path + run: echo "path=${{ runner.temp }}/data" >> $GITHUB_OUTPUT + outputs: + data_path: ${{ steps.data_path.outputs.path }} + + crds_contexts: + uses: spacetelescope/crds/.github/workflows/contexts.yml@master + + jwst: + uses: OpenAstronomy/github-actions-workflows/.github/workflows/tox.yml@d68193b68216da64eafaa618f53c59f5d271c56e # v1.14.0 + needs: [ environment, crds_contexts ] + with: + setenv: | + CRDS_PATH: ${{ needs.environment.outputs.data_path }}/crds + CRDS_SERVER_URL: https://jwst-crds.stsci.edu + CRDS_CLIENT_RETRY_COUNT: 3 + CRDS_CLIENT_RETRY_DELAY_SECONDS: 20 + cache-path: ${{ needs.environment.outputs.data_path }}/crds + cache-key: crds-${{ needs.crds_contexts.outputs.jwst }} + envs: | + - linux: py311-test-jwst-cov-xdist + + romancal: + uses: OpenAstronomy/github-actions-workflows/.github/workflows/tox.yml@d68193b68216da64eafaa618f53c59f5d271c56e # v1.14.0 + needs: [ environment, crds_contexts ] + with: + setenv: | + CRDS_PATH: ${{ needs.environment.outputs.data_path }}/crds + CRDS_SERVER_URL: https://roman-crds.stsci.edu + CRDS_CLIENT_RETRY_COUNT: 3 + CRDS_CLIENT_RETRY_DELAY_SECONDS: 20 + cache-path: ${{ needs.environment.outputs.data_path }}/crds + cache-key: crds-${{ needs.crds_contexts.outputs.jwst }} + envs: | + - linux: py311-test-romancal-cov-xdist + + romanisim_data: + needs: [ environment ] + uses: spacetelescope/romanisim/.github/workflows/retrieve_cache.yml@main + with: + cache_path: ${{ needs.environment.outputs.data_path }} + + romanisim: + needs: [ romanisim_data ] + uses: OpenAstronomy/github-actions-workflows/.github/workflows/tox.yml@d68193b68216da64eafaa618f53c59f5d271c56e # v1.14.0 + with: + libraries: | + brew: + - eigen + - fftw + setenv: | + WEBBPSF_PATH: ${{ needs.romanisim_data.outputs.cache_path }}/webbpsf-data/ + GALSIM_CAT_PATH: ${{ needs.romanisim_data.outputs.cache_path }}/galsim_data/real_galaxy_catalog_23.5_example.fits + FFTW_DIR: /opt/homebrew/opt/fftw/lib/ + cache-path: ${{ needs.romanisim_data.outputs.cache_path }} + cache-key: ${{ needs.romanisim_data.outputs.cache_key }} + envs: | + - linux: py311-test-romanisim-cov-xdist + + astropy: + uses: OpenAstronomy/github-actions-workflows/.github/workflows/tox.yml@d68193b68216da64eafaa618f53c59f5d271c56e # v1.14.0 + if: (github.repository == 'spacetelescope/gwcs' && (github.event_name == 'schedule' || github.event_name == 'push' || github.event_name == 'workflow_dispatch' || contains(github.event.pull_request.labels.*.name, 'Downstream CI'))) + with: + submodules: false + # Any env name which does not start with `pyXY` will use this Python version. + default_python: '3.11' + envs: | + - linux: specutils + + third-party: + uses: OpenAstronomy/github-actions-workflows/.github/workflows/tox.yml@d68193b68216da64eafaa618f53c59f5d271c56e # v1.14.0 + if: (github.repository == 'spacetelescope/gwcs' && (github.event_name == 'schedule' || github.event_name == 'push' || github.event_name == 'workflow_dispatch' || contains(github.event.pull_request.labels.*.name, 'Downstream CI'))) + with: + submodules: false + # Any env name which does not start with `pyXY` will use this Python version. + default_python: '3.11' + envs: | + - linux: ndcube + - linux: dkist diff --git a/.github/workflows/romanisim_data.yml b/.github/workflows/romanisim_data.yml new file mode 100644 index 00000000..8d67697b --- /dev/null +++ b/.github/workflows/romanisim_data.yml @@ -0,0 +1,18 @@ +name: download and cache romanisim data + +on: + schedule: + - cron: "42 4 * * 3" + workflow_dispatch: + inputs: + webbpsf_minimal: + description: minimal WebbPSF dataset + type: boolean + required: false + default: true + +jobs: + download_romanisim_data: + uses: spacetelescope/romanisim/.github/workflows/data.yml@main + with: + webbpsf_minimal: ${{ github.event_name != 'workflow_dispatch' && true || inputs.webbpsf_minimal }} diff --git a/CHANGES.rst b/CHANGES.rst index b26657ad..4a62e01e 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,3 +1,34 @@ +0.22.0 (2024-12-19) +------------------- + +- Replace usages of ``copy_arrays`` with ``memmap`` [#503] + +- Fix an issue with units in ``wcs_from_points``. [#507] + +- Fix incorrect units being returned in the low level WCS API. [#512] + +- Synchronize ``region.py`` with the copies of it in JWST and Romancal. [#517] + +- Add support for compound bounding boxes and ignored bounding box entries. [#519] + + +- Add ``gwcs.examples`` module, based on the examples located in the testing ``conftest.py``. [#521] + +- Force ``bounding_box`` to always be returned as a ``F`` ordered box. [#522] + +- Move the bounding box attachment to the forward transform property. [#532] + +- Adjust ``world_to_array_index_values`` to round to integer coordinates as specified by APE 14. [#525] + +- Add warning filter to asdf extension to prevent the ``bounding_box`` order warning for gwcs objects originating from a file. [#526] + +- Fixed a bug where evaluating the inverse transform did not + respect the bounding box. [#498] + +- Improved reliability of inside/outside footprint computations when evaluating + inverse transform with bounding box. [#536] + + 0.21.0 (2024-03-10) ------------------- diff --git a/README.rst b/README.rst index 5fca535d..cb6bd6fc 100644 --- a/README.rst +++ b/README.rst @@ -43,7 +43,7 @@ To clone from github and install the master branch:: git clone https://github.com/spacetelescope/gwcs.git cd gwcs - python setup.py install + pip install --editable . Contributing Code, Documentation, or Feedback diff --git a/docs/gwcs/points_to_wcs.rst b/docs/gwcs/points_to_wcs.rst index d9d0c05c..e20b4de8 100644 --- a/docs/gwcs/points_to_wcs.rst +++ b/docs/gwcs/points_to_wcs.rst @@ -76,7 +76,10 @@ can be used to convert coordinates from pixel to world. >>> gwcs_obj(36.235,642.215) # doctest: +FLOAT_CMP (246.72158004206716, 43.46075091731673) -Or equivalently - >>> gwcs_obj.forward_transform(36.235,642.215) # doctest: +FLOAT_CMP +Or using the common WCS API + >>> gwcs_obj.pixel_to_world_values(36.235,642.215) # doctest: +FLOAT_CMP (246.72158004206716, 43.46075091731673) + >>> gwcs_obj.pixel_to_world(36.235,642.215) # doctest: +FLOAT_CMP + \ No newline at end of file diff --git a/docs/gwcs/using_wcs.rst b/docs/gwcs/using_wcs.rst index a8430c58..c3ff9d26 100644 --- a/docs/gwcs/using_wcs.rst +++ b/docs/gwcs/using_wcs.rst @@ -69,7 +69,43 @@ Calling the :meth:`~gwcs.WCS.footprint` returns the footprint on the sky. .. doctest-skip:: >>> wcsobj.footprint() + array([[ 5.5264392 , -72.05149213], + ... [ 5.54580408, -72.06422321], + ... [ 5.63010439, -72.05426843], + ... [ 5.610708 , -72.04173847]]) +.. warning:: + + GWCS and astropy default to different tuple ordering conventions for representing + multi-dimensional bounding boxes. + + * GWCS uses the ``"F"`` ordering convention, where the tuples are ordered + ``((x0min, x0max), (x1min, x1max), ..., (xnmin, xnmax))`` (x,y,z ordering). + * While astropy uses the ``"C"`` ordering convention, where tuples are ordered + ``((xnmin, xnmax), ..., (x1min, x1max), (x0min, x0max))`` (z, y, x ordering). + + This means that given the same tuple of tuples, say ``((a, b), (c, d))``, setting + the bounding box on the transform prior to creating the GWCS will result in a + different bounding box than if one sets the same tuple of tuples on the GWCS object + itself. Indeed, in this case the former will assume ``(c, d)`` is the bounding box + for ``x`` while the latter will assume ``(a, b)`` is the bounding box for ``x``. + + It is recommended that when working on GWCS objects that one sets the bounding + box on the GWCS object itself, rather than on the transform prior to creating + the GWCS object. + + Note if one wants to set the bounding box on the transform itself + rather than the GWCS object then it should be done with + `~astropy.modeling.bind_bounding_box` with the ``order`` argument properly set. + + +.. note :: + + The GWCS will always convert or assume the bounding box to the ``"F"`` ordering + convention when setting the bounding box on the GWCS object itself and will + perform this conversion on the first access to the bounding box through the GWCS + object. If conversion occurs on first access, GWCS will issue a warning to alert + the user that the bounding box has been converted. Manipulating Transforms ----------------------- @@ -82,7 +118,7 @@ Transforms between frames can be retrieved and evaluated separately. >>> dist = wcsobj.get_transform('detector', 'undistorted_frame') >>> dist(1, 2) # doctest: +FLOAT_CMP - (-292.4150238489997, -616.8680129899999) + (41.62325692108675, -12.68101006210054) Transforms in the pipeline can be replaced by new transforms. @@ -119,6 +155,8 @@ inverse transformation converts world coordinates back to image coordinates. .. doctest-skip:: + >>> import asdf + >>> from astropy.utils.data import get_pkg_data_filename >>> wcsobj = asdf.open(get_pkg_data_filename('imaging_wcs_wdist.asdf')).tree['wcs'] The most general method available for computing inverse coordinate @@ -140,7 +178,7 @@ most situations: (349.9999994163172, 200.00000017679295) >>> world = wcsobj([2, 350, -5000], [2, 200, 6000]) >>> print(wcsobj.invert(*world)) # convert multiple points at once - (array([ 2.00000000e+00, 3.49999999e+02, -5.00000000e+03]), array([1.99999972e+00, 2.00000002e+02, 6.00000000e+03]) + (array([ 1.99999752, 3.49999999e+02, -5.00000000e+03]), array([1.99999972e+00, 2.00000002e+02, 6.00000000e+03]) By default, parameter ``quiet`` is set to `True` in `WCS.numerical_inverse() ` and so it will return results "as is" without warning us about possible loss @@ -178,4 +216,3 @@ point far away from the image for which numerical inverse fails. [5.31656943e-06 2.72052603e-10] [6.81557583e-06 1.06560533e-06] [3.96365344e-04 6.41822468e-05]] - diff --git a/docs/index.rst b/docs/index.rst index 7a6ad21f..481f2f76 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -137,7 +137,7 @@ using `ASDF Standard `__ There are two ways to save the GWCS object to a files: -- `Save a WCS object as a pure ASDF file`_ +- `Save a WCS object as a pure ASDF file`_ - `Save a WCS object as a pure ASDF file`_ @@ -222,11 +222,11 @@ expressed with no associated transform ... unit=(u.deg, u.deg)) >>> wcsobj = wcs.WCS([(detector_frame, det2sky), ... (sky_frame, None) - ... ] + ... ]) >>> print(wcsobj) From Transform -------- ---------------- - detector linear_transform + detector detector_to_sky icrs None To convert a pixel (x, y) = (1, 2) to sky coordinates, call the WCS object as a function: @@ -235,16 +235,27 @@ To convert a pixel (x, y) = (1, 2) to sky coordinates, call the WCS object as a >>> sky = wcsobj(1, 2) >>> print(sky) - + (29.980402161089177, 44.98616499109102) The :meth:`~gwcs.wcs.WCS.invert` method evaluates the :meth:`~gwcs.wcs.WCS.backward_transform` if available, otherwise applies an iterative method to calculate the reverse coordinates. .. doctest-skip:: - >>> wcsobj.invert(sky) - (, ) + >>> wcsobj.invert(*sky) + (0.9999999996185807, 1.999999999186798) + +GWCS supports the common WCS interface which defines several methods +to work with high level Astropy objects: + +.. doctest-skip:: + + >>> sky_obj = wcsobj.pixel_to_world(1, 2) + >>> print(sky) + + >>> wcsobj.world_to_pixel(sky_obj) + (0.9999999996185807, 1.999999999186798) .. _save_as_asdf: diff --git a/gwcs/api.py b/gwcs/api.py index 4f2ce9fc..32f6f9c3 100644 --- a/gwcs/api.py +++ b/gwcs/api.py @@ -132,7 +132,13 @@ def world_to_pixel_values(self, *world_arrays): be returned in the ``(x, y)`` order, where for an image, ``x`` is the horizontal coordinate and ``y`` is the vertical coordinate. """ - world_arrays = self._add_units_input(world_arrays, self.backward_transform, self.output_frame) + try: + backward_transform = self.backward_transform + world_arrays = self._add_units_input(world_arrays, + backward_transform, + self.output_frame) + except NotImplementedError: + pass result = self.invert(*world_arrays, with_units=False) @@ -147,10 +153,14 @@ def world_to_array_index_values(self, *world_arrays): `~BaseLowLevelWCS.pixel_to_world_values`). The indices should be returned as rounded integers. """ - result = self.world_to_pixel_values(*world_arrays) - if self.pixel_n_dim != 1: - result = result[::-1] - return result + results = self.world_to_pixel_values(*world_arrays) + if self.pixel_n_dim == 1: + results = (results,) + else: + results = results[::-1] + + results = tuple(utils._toindex(result) for result in results) + return results[0] if self.pixel_n_dim == 1 else results @property def array_shape(self): @@ -313,7 +323,6 @@ def world_to_pixel(self, *world_objects): Convert world coordinates to pixel values. """ result = self.invert(*world_objects, with_units=True) - if self.input_frame.naxes > 1: first_res = result[0] if not utils.isnumerical(first_res): diff --git a/gwcs/converters/wcs.py b/gwcs/converters/wcs.py index a2f40b3f..313d6d10 100644 --- a/gwcs/converters/wcs.py +++ b/gwcs/converters/wcs.py @@ -1,6 +1,8 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst # -*- coding: utf-8 -*- +import warnings +from contextlib import suppress from asdf.extension import Converter @@ -14,10 +16,17 @@ class WCSConverter(Converter): types = ["gwcs.wcs.WCS"] def from_yaml_tree(self, node, tag, ctx): - from ..wcs import WCS + from ..wcs import WCS, GwcsBoundingBoxWarning gwcsobj = WCS(node['steps'], name=node['name']) if 'pixel_shape' in node: gwcsobj.pixel_shape = node['pixel_shape'] + + # Ignore the warning about the bounding box order for data read from a + # file. This is causing issues with files from MAST. + with suppress(AttributeError), warnings.catch_warnings(): + warnings.filterwarnings('ignore', category=GwcsBoundingBoxWarning) + _ = gwcsobj.bounding_box + return gwcsobj def to_yaml_tree(self, gwcsobj, tag, ctx): diff --git a/gwcs/coordinate_frames.py b/gwcs/coordinate_frames.py index 5131a4d4..5c36c253 100644 --- a/gwcs/coordinate_frames.py +++ b/gwcs/coordinate_frames.py @@ -355,8 +355,8 @@ def _world_axis_object_classes(self): @property def _world_axis_object_components(self): - return [('celestial', 0, 'spherical.lon'), - ('celestial', 1, 'spherical.lat')] + return [('celestial', 0, lambda sc: sc.spherical.lon.to_value(self.unit[0])), + ('celestial', 1, lambda sc: sc.spherical.lat.to_value(self.unit[1]))] def coordinates(self, *args): """ @@ -457,7 +457,7 @@ def _world_axis_object_classes(self): @property def _world_axis_object_components(self): - return [('spectral', 0, 'value')] + return [('spectral', 0, lambda sc: sc.to_value(self.unit[0]))] def coordinates(self, *args): # using SpectralCoord diff --git a/gwcs/examples.py b/gwcs/examples.py new file mode 100644 index 00000000..cd5713ae --- /dev/null +++ b/gwcs/examples.py @@ -0,0 +1,512 @@ +import numpy as np + +import astropy.units as u +from astropy.time import Time +from astropy import coordinates as coord +from astropy.modeling import models + +from . import coordinate_frames as cf +from . import spectroscopy as sp +from . import wcs + +# frames +DETECTOR_1D_FRAME = cf.CoordinateFrame(name='detector', axes_order=(0,), naxes=1, axes_type="detector") +DETECTOR_2D_FRAME = cf.Frame2D(name='detector', axes_order=(0, 1)) +ICRC_SKY_FRAME = cf.CelestialFrame(reference_frame=coord.ICRS(), + axes_order=(0, 1)) + +FREQ_FRAME = cf.SpectralFrame(name='freq', unit=u.Hz, axes_order=(0, )) +WAVE_FRAME = cf.SpectralFrame(name='wave', unit=u.m, axes_order=(2, ), + axes_names=('lambda', )) + +# transforms +MODEL_2D_SHIFT = models.Shift(1) & models.Shift(2) +MODEL_1D_SCALE = models.Scale(2) + + +def gwcs_2d_quantity_shift(): + frame = cf.CoordinateFrame(name="quantity", axes_order=(0, 1), naxes=2, axes_type=("SPATIAL", "SPATIAL"), unit=(u.km, u.km)) + pipe = [(DETECTOR_2D_FRAME, MODEL_2D_SHIFT), (frame, None)] + + return wcs.WCS(pipe) + + +def gwcs_2d_spatial_shift(): + """ + A simple one step spatial WCS, in ICRS with a 1 and 2 px shift. + """ + pipe = [(DETECTOR_2D_FRAME, MODEL_2D_SHIFT), (ICRC_SKY_FRAME, None)] + return wcs.WCS(pipe) + + +def gwcs_2d_spatial_reordered(): + """ + A simple one step spatial WCS, in ICRS with a 1 and 2 px shift. + """ + out_frame = cf.CelestialFrame(reference_frame=coord.ICRS(), + axes_order=(1, 0)) + return wcs.WCS(MODEL_2D_SHIFT | models.Mapping((1, 0)), input_frame=DETECTOR_2D_FRAME, output_frame=out_frame) + + +def gwcs_1d_freq(): + return wcs.WCS([(DETECTOR_1D_FRAME, MODEL_1D_SCALE), (FREQ_FRAME, None)]) + + +def gwcs_3d_spatial_wave(): + comp1 = cf.CompositeFrame([ICRC_SKY_FRAME, WAVE_FRAME]) + m = MODEL_2D_SHIFT & MODEL_1D_SCALE + + detector_frame = cf.CoordinateFrame(name="detector", naxes=3, + axes_order=(0, 1, 2), + axes_type=("pixel", "pixel", "pixel"), + axes_names=("x", "y", "z"), unit=(u.pix, u.pix, u.pix)) + + return wcs.WCS([(detector_frame, m), + (comp1, None)]) + + +def gwcs_2d_shift_scale(): + m1 = models.Shift(1) & models.Shift(2) + m2 = models.Scale(5) & models.Scale(10) + m3 = m1 | m2 + pipe = [(DETECTOR_2D_FRAME, m3), (ICRC_SKY_FRAME, None)] + return wcs.WCS(pipe) + +def gwcs_2d_bad_bounding_box_order(): + m1 = models.Shift(1) & models.Shift(2) + m2 = models.Scale(5) & models.Scale(10) + m3 = m1 | m2 + + # Purposefully set the bounding box in the wrong order + m3.bounding_box = ((1, 2), (3, 4)) + + pipe = [(DETECTOR_2D_FRAME, m3), (ICRC_SKY_FRAME, None)] + return wcs.WCS(pipe) + + +def gwcs_1d_freq_quantity(): + + detector_1d = cf.CoordinateFrame(name='detector', axes_order=(0,), naxes=1, unit=u.pix, axes_type="detector") + return wcs.WCS([(detector_1d, models.Multiply(1 * u.Hz / u.pix)), (FREQ_FRAME, None)]) + + +def gwcs_2d_shift_scale_quantity(): + m4 = models.Shift(1 * u.pix) & models.Shift(2 * u.pix) + m5 = models.Scale(5 * u.deg) + m6 = models.Scale(10 * u.deg) + m5.input_units_equivalencies = {'x': u.pixel_scale(1 * u.deg / u.pix)} + m6.input_units_equivalencies = {'x': u.pixel_scale(1 * u.deg / u.pix)} + m5.inverse = models.Scale(1. / 5 * u.pix) + m6.inverse = models.Scale(1. / 10 * u.pix) + m5.inverse.input_units_equivalencies = { + 'x': u.pixel_scale(1 * u.pix / u.deg) + } + m6.inverse.input_units_equivalencies = { + 'x': u.pixel_scale(1 * u.pix / u.deg) + } + m7 = m5 & m6 + m8 = m4 | m7 + pipe2 = [(DETECTOR_2D_FRAME, m8), (ICRC_SKY_FRAME, None)] + return wcs.WCS(pipe2) + + +def gwcs_3d_identity_units(): + """ + A simple 1-1 gwcs that converts from pixels to arcseconds + """ + identity = (models.Multiply(1 * u.arcsec / u.pixel) & + models.Multiply(1 * u.arcsec / u.pixel) & + models.Multiply(1 * u.nm / u.pixel)) + sky_frame = cf.CelestialFrame(axes_order=(0, 1), name='icrs', + reference_frame=coord.ICRS(), + axes_names=("longitude", "latitude"), + unit=(u.arcsec, u.arcsec)) + wave_frame = cf.SpectralFrame(axes_order=(2, ), unit=u.nm, axes_names=("wavelength",)) + + frame = cf.CompositeFrame([sky_frame, wave_frame]) + + detector_frame = cf.CoordinateFrame(name="detector", naxes=3, + axes_order=(0, 1, 2), + axes_type=("pixel", "pixel", "pixel"), + axes_names=("x", "y", "z"), unit=(u.pix, u.pix, u.pix)) + + return wcs.WCS(forward_transform=identity, output_frame=frame, input_frame=detector_frame) + + +def gwcs_4d_identity_units(): + """ + A simple 1-1 gwcs that converts from pixels to arcseconds + """ + identity = (models.Multiply(1*u.arcsec/u.pixel) & models.Multiply(1*u.arcsec/u.pixel) & + models.Multiply(1*u.nm/u.pixel) & models.Multiply(1*u.s/u.pixel)) + sky_frame = cf.CelestialFrame(axes_order=(0, 1), name='icrs', + reference_frame=coord.ICRS()) + wave_frame = cf.SpectralFrame(axes_order=(2, ), unit=u.nm) + time_frame = cf.TemporalFrame(axes_order=(3, ), unit=u.s, + reference_frame=Time("2000-01-01T00:00:00")) + + frame = cf.CompositeFrame([sky_frame, wave_frame, time_frame]) + + detector_frame = cf.CoordinateFrame(name="detector", naxes=4, + axes_order=(0, 1, 2, 3), + axes_type=("pixel", "pixel", "pixel", "pixel"), + axes_names=("x", "y", "z", "s"), unit=(u.pix, u.pix, u.pix, u.pix)) + + return wcs.WCS(forward_transform=identity, output_frame=frame, input_frame=detector_frame) + + +def gwcs_simple_imaging_units(): + shift_by_crpix = models.Shift(-2048*u.pix) & models.Shift(-1024*u.pix) + matrix = np.array([[1.290551569736E-05, 5.9525007864732E-06], + [5.0226382102765E-06 , -1.2644844123757E-05]]) + rotation = models.AffineTransformation2D(matrix * u.deg, + translation=[0, 0] * u.deg) + rotation.input_units_equivalencies = {"x": u.pixel_scale(1*u.deg/u.pix), + "y": u.pixel_scale(1*u.deg/u.pix)} + rotation.inverse = models.AffineTransformation2D(np.linalg.inv(matrix) * u.pix, + translation=[0, 0] * u.pix) + rotation.inverse.input_units_equivalencies = {"x": u.pixel_scale(1*u.pix/u.deg), + "y": u.pixel_scale(1*u.pix/u.deg)} + tan = models.Pix2Sky_TAN() + celestial_rotation = models.RotateNative2Celestial(5.63056810618*u.deg, + -72.05457184279*u.deg, + 180*u.deg) + det2sky = shift_by_crpix | rotation | tan | celestial_rotation + det2sky.name = "linear_transform" + + detector_frame = cf.Frame2D(name="detector", axes_names=("x", "y"), + unit=(u.pix, u.pix)) + sky_frame = cf.CelestialFrame(reference_frame=coord.ICRS(), name='icrs', + unit=(u.deg, u.deg)) + pipeline = [(detector_frame, det2sky), + (sky_frame, None) + ] + return wcs.WCS(pipeline) + + +def gwcs_simple_imaging(): + shift_by_crpix = models.Shift(-2048) & models.Shift(-1024) + matrix = np.array([[1.290551569736E-05, 5.9525007864732E-06], + [5.0226382102765E-06 , -1.2644844123757E-05]]) + rotation = models.AffineTransformation2D(matrix, + translation=[0, 0]) + rotation.inverse = models.AffineTransformation2D(np.linalg.inv(matrix), + translation=[0, 0]) + tan = models.Pix2Sky_TAN() + celestial_rotation = models.RotateNative2Celestial(5.63056810618, + -72.05457184279, + 180) + det2sky = shift_by_crpix | rotation | tan | celestial_rotation + det2sky.name = "linear_transform" + + detector_frame = cf.Frame2D(name="detector", axes_names=("x", "y")) + sky_frame = cf.CelestialFrame(reference_frame=coord.ICRS(), name='icrs') + pipeline = [(detector_frame, det2sky), + (sky_frame, None) + ] + return wcs.WCS(pipeline) + + +def gwcs_stokes_lookup(): + transform = models.Tabular1D([0, 1, 2, 3] * u.pix, [1, 2, 3, 4] * u.one, + method="nearest", fill_value=np.nan, bounds_error=False) + frame = cf.StokesFrame() + + detector_frame = cf.CoordinateFrame(name="detector", naxes=1, + axes_order=(0,), + axes_type=("pixel",), + axes_names=("x",), unit=(u.pix,)) + + return wcs.WCS(forward_transform=transform, output_frame=frame, input_frame=detector_frame) + + +def gwcs_3spectral_orders(): + comp1 = cf.CompositeFrame([ICRC_SKY_FRAME, WAVE_FRAME]) + detector_frame = cf.Frame2D(name="detector", axes_names=("x", "y"), + unit=(u.pix, u.pix)) + m = MODEL_2D_SHIFT & MODEL_1D_SCALE + + return wcs.WCS([(detector_frame, m), + (comp1, None)]) + + +def gwcs_with_frames_strings(): + transform = models.Shift(1) & models.Shift(1) & models.Polynomial2D(1) + pipe = [('detector', transform), + ('world', None) + ] + return wcs.WCS(pipe) + + +def sellmeier_glass(): + B_coef = [0.58339748, 0.46085267, 3.8915394] + C_coef = [0.00252643, 0.010078333, 1200.556] + return sp.SellmeierGlass(B_coef, C_coef) + + +def sellmeier_zemax(): + B_coef = [0.58339748, 0.46085267, 3.8915394] + C_coef = [0.00252643, 0.010078333, 1200.556] + D_coef = [-2.66e-05, 0.0, 0.0] + E_coef = [0., 0., 0.] + return sp.SellmeierZemax(65, 35, 0, 0, B_coef = B_coef, + C_coef=C_coef, D_coef=D_coef, + E_coef=E_coef) + + +def gwcs_3d_galactic_spectral(): + """ + This fixture has the axes ordered as lat, spectral, lon. + """ + # lat,wav,lon + crpix1, crpix2, crpix3 = 29, 39, 44 + crval1, crval2, crval3 = 10, 20, 25 + cdelt1, cdelt2, cdelt3 = -0.01, 0.5, 0.01 + + shift = models.Shift(-crpix3) & models.Shift(-crpix1) + scale = models.Multiply(cdelt3) & models.Multiply(cdelt1) + proj = models.Pix2Sky_CAR() + skyrot = models.RotateNative2Celestial(crval3, 90 + crval1, 180) + celestial = shift | scale | proj | skyrot + + wave_model = models.Shift(-crpix2) | models.Multiply(cdelt2) | models.Shift(crval2) + + transform = models.Mapping((2, 0, 1)) | celestial & wave_model | models.Mapping((1, 2, 0)) + + sky_frame = cf.CelestialFrame(axes_order=(2, 0), + reference_frame=coord.Galactic(), axes_names=("Longitude", "Latitude")) + wave_frame = cf.SpectralFrame(axes_order=(1, ), unit=u.Hz, axes_names=("Frequency",)) + + frame = cf.CompositeFrame([sky_frame, wave_frame]) + + detector_frame = cf.CoordinateFrame(name="detector", naxes=3, + axes_order=(0, 1, 2), + axes_type=("pixel", "pixel", "pixel"), + unit=(u.pix, u.pix, u.pix)) + + owcs = wcs.WCS(forward_transform=transform, output_frame=frame, input_frame=detector_frame) + owcs.bounding_box = ((-1, 35), (-2, 45), (5, 50)) + owcs.array_shape = (30, 20, 10) + owcs.pixel_shape = (10, 20, 30) + + return owcs + + +def gwcs_1d_spectral(): + """ + A simple 1D spectral WCS. + """ + wave_model = models.Shift(-5) | models.Multiply(3.7) | models.Shift(20) + wave_model.bounding_box = (7, 50) + wave_frame = cf.SpectralFrame(axes_order=(0, ), unit=u.Hz, axes_names=("Frequency",)) + + detector_frame = cf.CoordinateFrame( + name="detector", naxes=1, axes_order=(0, ), + axes_type=("pixel",), unit=(u.pix, ) + ) + + owcs = wcs.WCS(forward_transform=wave_model, output_frame=wave_frame, input_frame=detector_frame) + owcs.array_shape = (44, ) + owcs.pixel_shape = (44, ) + + return owcs + + +def gwcs_spec_cel_time_4d(): + """ + A complex 4D mixed celestial + spectral + time WCS. + """ + # spectroscopic frame: + wave_model = models.Shift(-5) | models.Multiply(3.7) | models.Shift(20) + wave_model.bounding_box = (7, 50) + wave_frame = cf.SpectralFrame(name='wave', unit=u.m, axes_order=(0,), axes_names=('lambda',)) + + # time frame: + time_model = models.Identity(1) # models.Linear1D(10, 0) + time_frame = cf.TemporalFrame(Time("2010-01-01T00:00"), name='time', unit=u.s, axes_order=(3,)) + + # Values from data/acs.hdr: + crpix = (12, 13) + crval = (5.63, -72.05) + cd = [[1.291E-05, 5.9532E-06], [5.02215E-06, -1.2645E-05]] + aff = models.AffineTransformation2D(matrix=cd, name='rotation') + offx = models.Shift(-crpix[0], name='x_translation') + offy = models.Shift(-crpix[1], name='y_translation') + wcslin = models.Mapping((1, 0)) | (offx & offy) | aff + tan = models.Pix2Sky_TAN(name='tangent_projection') + n2c = models.RotateNative2Celestial(*crval, 180, name='sky_rotation') + cel_model = wcslin | tan | n2c + icrs = cf.CelestialFrame(reference_frame=coord.ICRS(), name='sky', axes_order=(2, 1)) + + wcs_forward = wave_model & cel_model & time_model + + comp_frm = cf.CompositeFrame(frames=[wave_frame, icrs, time_frame], name='TEST 4D FRAME') + + detector_frame = cf.CoordinateFrame( + name="detector", naxes=4, axes_order=(0, 1, 2, 3), + axes_type=("pixel", "pixel", "pixel", "pixel"), + unit=(u.pix, u.pix, u.pix, u.pix) + ) + + w = wcs.WCS(forward_transform=wcs_forward, output_frame=comp_frm, input_frame=detector_frame) + + w.bounding_box = ((0, 63), (0, 127), (0, 255), (0, 9)) + w.array_shape = (10, 256, 128, 64) + w.pixel_shape = (64, 128, 256, 10) + return w + + +def gwcs_cube_with_separable_spectral(axes_order): + """ + GWCS cube with spectral axis separable from the celestial axes. + + Viable examples are (2, 0, 1) and (2, 1, 0). + (1, 0, 2) fails round-trip for -TAB axis 3 + """ + cube_size = (128, 64, 100) + + spectral_axes_order = (axes_order.index(2), ) + cel_axes_order = (axes_order.index(0), axes_order.index(1)) + + # Values from data/acs.hdr: + crpix = (64, 32) + crval = (5.63056810618, -72.0545718428) + cd = [[1.29058667557984E-05, 5.95320245884555E-06], + [5.02215195623825E-06, -1.2645010396976E-05]] + + aff = models.AffineTransformation2D(matrix=cd, name='rotation') + offx = models.Shift(-crpix[0], name='x_translation') + offy = models.Shift(-crpix[1], name='y_translation') + + wcslin = (offx & offy) | aff + tan = models.Pix2Sky_TAN(name='tangent_projection') + n2c = models.RotateNative2Celestial(*crval, 180, name='sky_rotation') + icrs = cf.CelestialFrame(reference_frame=coord.ICRS(), name='sky', + axes_order=cel_axes_order) + spec = cf.SpectralFrame( + name='wave', unit=[u.m,], axes_order=spectral_axes_order, + axes_names=('lambda',) + ) + comp_frm = cf.CompositeFrame(frames=[icrs, spec], name='TEST 3D FRAME WITH SPECTRAL AXIS') + wcs_forward = ((wcslin & models.Identity(1)) | + (tan & models.Identity(1)) | + (n2c & models.Identity(1)) | + models.Mapping(axes_order)) + + detector_frame = cf.CoordinateFrame(name="detector", naxes=3, + axes_order=(0, 1, 2), + axes_type=("pixel", "pixel", "pixel"), + unit=(u.pix, u.pix, u.pix)) + + w = wcs.WCS(forward_transform=wcs_forward, output_frame=comp_frm, + input_frame=detector_frame) + w.bounding_box = tuple((0, k - 1) for k in cube_size) + w.pixel_shape = cube_size + w.array_shape = w.pixel_shape[::-1] + + return w, axes_order + + +def gwcs_cube_with_separable_time(axes_order): + """ + A mixed celestial + time WCS. + + Viable examples are (2, 0, 1) and (2, 1, 0). + (0, 2, 1) fails round-trip for -TAB axis 2 + (1, 0, 2) fails round-trip for -TAB axis 3 + """ + cube_size = (64, 32, 128) + + time_axes_order = (axes_order.index(2), ) + cel_axes_order = (axes_order.index(0), axes_order.index(1)) + + detector_frame = cf.CoordinateFrame( + name="detector", naxes=3, axes_order=(0, 1, 2), + axes_type=("pixel", "pixel", "pixel"), unit=(u.pix, u.pix, u.pix) + ) + + # time frame: + time_model = models.Identity(1) # models.Linear1D(10, 0) + time_frame = cf.TemporalFrame(Time("2010-01-01T00:00"), name='time', + unit=u.s, axes_order=time_axes_order) + + # Values from data/acs.hdr: + crpix = (12, 13) + crval = (5.63, -72.05) + cd = [[1.291E-05, 5.9532E-06], [5.02215E-06, -1.2645E-05]] + aff = models.AffineTransformation2D(matrix=cd, name='rotation') + offx = models.Shift(-crpix[0], name='x_translation') + offy = models.Shift(-crpix[1], name='y_translation') + wcslin = models.Mapping((1, 0)) | (offx & offy) | aff + tan = models.Pix2Sky_TAN(name='tangent_projection') + n2c = models.RotateNative2Celestial(*crval, 180, name='sky_rotation') + cel_model = wcslin | tan | n2c + icrs = cf.CelestialFrame(reference_frame=coord.ICRS(), name='sky', + axes_order=cel_axes_order) + + wcs_forward = (cel_model & time_model) | models.Mapping(axes_order) + + comp_frm = cf.CompositeFrame(frames=[icrs, time_frame], + name='TEST 3D FRAME WITH TIME') + + w = wcs.WCS(forward_transform=wcs_forward, output_frame=comp_frm, + input_frame=detector_frame) + + w.bounding_box = tuple((0, k - 1) for k in cube_size) + w.pixel_shape = cube_size + w.array_shape = w.pixel_shape[::-1] + return w + + +def gwcs_7d_complex_mapping(): + """ + Useful features of this WCS (axes indices here are 0-based): + - includes two celestial axes: input (0, 1) maps to world (2 - RA, 1 - Dec) + - includes one separable frame with one axis: 4 -> 2 + - includes one frame with 3 input and 4 output axes (1 degenerate), + with separable world axes (3, 5) and (0, 6). + """ + offx = models.Shift(-64, name='x_translation') + offy = models.Shift(-32, name='y_translation') + cd = np.array([[1.2906, 0.59532], [0.50222, -1.2645]]) + aff = models.AffineTransformation2D(matrix=1e-5 * cd, name='rotation') + aff2 = models.AffineTransformation2D(matrix=cd, name='rotation2') + + wcslin = (offx & offy) | aff + tan = models.Pix2Sky_TAN(name='tangent_projection') + n2c = models.RotateNative2Celestial(5.630568, -72.0546, 180, name='skyrot') + icrs = cf.CelestialFrame(reference_frame=coord.ICRS(), name='sky', axes_order=(2, 1)) + spec = cf.SpectralFrame(name='wave', unit=[u.m], axes_order=(4,), axes_names=('lambda',)) + cmplx = cf.CoordinateFrame( + name="complex", + naxes=4, + axes_order=(3, 5, 0, 6), + axis_physical_types=(['em.wl', 'em.wl', 'time', 'time']), + axes_type=("SPATIAL", "SPATIAL", "TIME", "TIME"), + axes_names=("x", "y", "t", 'tau'), + unit=(u.m, u.m, u.second, u.second) + ) + + comp_frm = cf.CompositeFrame(frames=[icrs, spec, cmplx], name='TEST 7D') + wcs_forward = ((wcslin & models.Shift(-3.14) & models.Scale(2.7) & aff2) | + (tan & models.Identity(1) & models.Identity(1) & models.Identity(2)) | + (n2c & models.Identity(1) & models.Identity(1) & models.Identity(2)) | + models.Mapping((3, 1, 0, 4, 2, 5, 3))) + + + detector_frame = cf.CoordinateFrame( + name="detector", naxes=6, + axes_order=(0, 1, 2, 3, 4, 5), + axes_type=("pixel", "pixel", "pixel", "pixel", "pixel", "pixel"), + unit=(u.pix, u.pix, u.pix, u.pix, u.pix, u.pix) + ) + + # pipeline = [('detector', wcs_forward), (comp_frm, None)] + w = wcs.WCS(forward_transform=wcs_forward, output_frame=comp_frm, + input_frame=detector_frame) + w.bounding_box = ((0, 15), (0, 31), (0, 20), (0, 10), (0, 10), (0, 1)) + + w.array_shape = (2, 11, 11, 21, 32, 16) + w.pixel_shape = (16, 32, 21, 11, 11, 2) + + return w diff --git a/gwcs/region.py b/gwcs/region.py index 29665505..8faeaee6 100644 --- a/gwcs/region.py +++ b/gwcs/region.py @@ -14,9 +14,9 @@ __all__ = ['Region', 'Edge', 'Polygon'] +_INTERSECT_ATOL = 1e2 * np.finfo(float).eps -class Region(metaclass=abc.ABCMeta): - +class Region: """ Base class for regions. @@ -33,14 +33,14 @@ def __init__(self, rid, coordinate_frame): self._rid = rid @abc.abstractmethod - def __contains__(self, x, y): + def __contains__(self, px): """ Determines if a pixel is within a region. Parameters ---------- - x, y : float - x , y values of a pixel + px : tuple[float, float] + A pixel coordinate (x, y) Returns ------- @@ -49,6 +49,7 @@ def __contains__(self, x, y): Subclasses must define this method. """ + @abc.abstractmethod def scan(self, mask): """ Sets mask values to region id for all pixels within the region. @@ -68,7 +69,6 @@ def scan(self, mask): class Polygon(Region): - """ Represents a 2D polygon region with multiple vertices. @@ -86,11 +86,11 @@ class Polygon(Region): """ - def __init__(self, rid, vertices, coord_frame="detector"): + def __init__(self, rid, vertices, coord_system="Cartesian"): if len(vertices) < 4: - raise ValueError("Expected vertices to be " - "a list of minimum 4 tuples (x,y)") - super(Polygon, self).__init__(rid, coord_frame) + raise ValueError("Expected vertices to be a list of minimum 4 tuples (x,y)") + + super().__init__(rid, coord_system) # self._shiftx & self._shifty are introduced to shift the bottom-left # corner of the polygon's bounding box to (0,0) as a (hopefully @@ -114,8 +114,9 @@ def __init__(self, rid, vertices, coord_frame="detector"): self._shifty = int(round(self._shifty)) self._bbox = self._get_bounding_box() - self._scan_line_range = \ - list(range(self._bbox[1], self._bbox[3] + self._bbox[1] + 1)) + self._scan_line_range = list( + range(self._bbox[1], self._bbox[3] + self._bbox[1] + 1) + ) # constructs a Global Edge Table (GET) in bbox coordinates self._GET = self._construct_ordered_GET() @@ -130,7 +131,7 @@ def _construct_ordered_GET(self): """ Construct a Global Edge Table (GET) - The GET is an OrderedDict. Keys are scan line numbers, + The GET is an OrderedDict. Keys are scan line numbers, ordered from bbox.ymin to bbox.ymax, where bbox is the bounding box of the polygon. Values are lists of edges for which edge.ymin==scan_line_number. @@ -149,7 +150,7 @@ def _construct_ordered_GET(self): ymin_ind = (ymin == i).nonzero()[0] # a hack for incomplete filling .any() fails if 0 is in ymin_ind # if ymin_ind.any(): - yminindlen, = ymin_ind.shape + (yminindlen,) = ymin_ind.shape if yminindlen: GET[i] = [edges[ymin_ind[0]]] for j in ymin_ind[1:]: @@ -160,8 +161,10 @@ def get_edges(self): """ Create a list of Edge objects from vertices """ - return [Edge(name=f'E{i - 1}', start=self._vertices[i - 1], stop=self._vertices[i]) - for i in range(1, len(self._vertices))] + return [ + Edge(name=f"E{i - 1}", start=self._vertices[i - 1], stop=self._vertices[i]) + for i in range(1, len(self._vertices)) + ] def scan(self, data): """ @@ -186,11 +189,13 @@ def scan(self, data): 5. Set elements between pairs of X in the AET to the Edge's ID """ - # TODO: 1.This algorithm does not mark pixels in the top row and left most column. - # Pad the initial pixel description on top and left with 1 px to prevent this. - # 2. Currently it uses intersection of the scan line with edges. If this is - # too slow it should use the 1/m increment (replace 3 above) (or the increment - # should be removed from the GET entry). + # TODO: + # 1. This algorithm does not mark pixels in the top row and left most + # column. Pad the initial pixel description on top and left with 1 px + # to prevent this. + # 2. Currently it uses # intersection of the scan line with edges. If + # this is too slow it should use the 1/m increment (replace 3 above) + # (or the increment should be removed from the GET entry). # see comments in the __init__ function for the reason of introducing # polygon shifts (self._shiftx & self._shifty). Here we need to shift @@ -204,7 +209,6 @@ def scan(self, data): scline = self._scan_line_range[-1] while y <= scline: - if y < scline: AET = self.update_AET(y, AET) @@ -212,10 +216,16 @@ def scan(self, data): y += 1 continue - scan_line = Edge('scan_line', start=[self._bbox[0], y], - stop=[self._bbox[0] + self._bbox[2], y]) - x = [int(np.ceil(e.compute_AET_entry(scan_line)[1])) - for e in AET if e is not None] + scan_line = Edge( + "scan_line", + start=[self._bbox[0], y], + stop=[self._bbox[0] + self._bbox[2], y], + ) + x = [ + int(np.ceil(e.compute_AET_entry(scan_line)[1])) + for e in AET + if e is not None + ] xnew = np.sort(x) ysh = y + self._shifty @@ -226,7 +236,7 @@ def scan(self, data): for i, j in zip(xnew[::2], xnew[1::2]): xstart = max(0, i + self._shiftx) xend = min(j + self._shiftx, nx - 1) - data[ysh][xstart:xend + 1] = self._rid + data[ysh][xstart : xend + 1] = self._rid y += 1 @@ -254,13 +264,16 @@ def update_AET(self, y, AET): return AET def __contains__(self, px): - """even-odd algorithm or smth else better sould be used""" - return px[0] >= self._bbox[0] and px[0] <= self._bbox[0] + self._bbox[2] and \ - px[1] >= self._bbox[1] and px[1] <= self._bbox[1] + self._bbox[3] + """even-odd algorithm or something else better should be used""" + return ( + px[0] >= self._bbox[0] + and px[0] <= self._bbox[0] + self._bbox[2] + and px[1] >= self._bbox[1] + and px[1] <= self._bbox[1] + self._bbox[3] + ) class Edge: - """ Edge representation. @@ -271,7 +284,7 @@ class Edge: """ - def __init__(self, name=None, start=None, stop=None, next=None): + def __init__(self, name=None, start=None, stop=None, next=None): # noqa: A002 self._start = None if start is not None: self._start = np.asarray(start) @@ -329,8 +342,12 @@ def compute_GET_entry(self): if np.diff(earr[:, 1]).item() == 0: return None else: - entry = [self._ymax, self._yminx, - (np.diff(earr[:, 0]) / np.diff(earr[:, 1])).item(), None] + entry = [ + self._ymax, + self._yminx, + (np.diff(earr[:, 0]) / np.diff(earr[:, 1])).item(), + None, + ] return entry def compute_AET_entry(self, edge): @@ -347,19 +364,20 @@ def __repr__(self): fmt = "" if self._name is not None: fmt += self._name - next = self.next - while next is not None: + next_edge = self.next + while next_edge is not None: fmt += "-->" - fmt += next._name - next = next.next + fmt += next_edge._name + next_edge = next_edge.next + return fmt @property - def next(self): + def next(self): # noqa: A003 return self._next @next.setter - def next(self, edge): + def next(self, edge): # noqa: A003 if self._name is None: self._name = edge._name self._stop = edge._stop @@ -372,19 +390,31 @@ def intersection(self, edge): u = self._stop - self._start v = edge._stop - edge._start w = self._start - edge._start - D = np.cross(u, v) - if np.allclose(np.cross(u, v), 0, rtol=0, - atol=1e2 * np.finfo(float).eps): + D = _det(u, v) + + if abs(D) < _INTERSECT_ATOL: return np.array(self._start) - return np.cross(v, w) / D * u + self._start + return _det(v, w) / D * u + self._start def is_parallel(self, edge): u = self._stop - self._start v = edge._stop - edge._start - return np.allclose(np.cross(u, v), 0, rtol=0, - atol=1e2 * np.finfo(float).eps) + return np.allclose(_det(u, v), 0, rtol=0, atol=_INTERSECT_ATOL) + + +def _det(u, v): + """ + Find the determinant of the matrix formed by the vectors u and v + + Note: Originally this was computed using a numpy "2D" cross product, + however, this functionality has been deprecated and slated for + removal. + Note: This is marginally faster than using `np.linalg.det([u, v])` by ~10ms + during empirical testing + """ + return u[0] * v[1] - u[1] * v[0] def _round_vertex(v): diff --git a/gwcs/tests/conftest.py b/gwcs/tests/conftest.py index 0b0a8878..3bf243e5 100644 --- a/gwcs/tests/conftest.py +++ b/gwcs/tests/conftest.py @@ -3,349 +3,107 @@ """ import pytest -import numpy as np - -import astropy.units as u -from astropy.time import Time -from astropy import coordinates as coord -from astropy.modeling import models - -from .. import coordinate_frames as cf -from .. import spectroscopy as sp -from .. import wcs +from .. import examples from .. import geometry -# frames -detector_1d = cf.CoordinateFrame(name='detector', axes_order=(0,), naxes=1, axes_type="detector") -detector_2d = cf.Frame2D(name='detector', axes_order=(0, 1)) -icrs_sky_frame = cf.CelestialFrame(reference_frame=coord.ICRS(), - axes_order=(0, 1)) - -freq_frame = cf.SpectralFrame(name='freq', unit=u.Hz, axes_order=(0, )) -wave_frame = cf.SpectralFrame(name='wave', unit=u.m, axes_order=(2, ), - axes_names=('lambda', )) - -# transforms -model_2d_shift = models.Shift(1) & models.Shift(2) -model_1d_scale = models.Scale(2) - @pytest.fixture def gwcs_2d_quantity_shift(): - frame = cf.CoordinateFrame(name="quantity", axes_order=(0, 1), naxes=2, axes_type=("SPATIAL", "SPATIAL"), unit=(u.km, u.km)) - pipe = [(detector_2d, model_2d_shift), (frame, None)] - - return wcs.WCS(pipe) + return examples.gwcs_2d_quantity_shift() @pytest.fixture def gwcs_2d_spatial_shift(): - """ - A simple one step spatial WCS, in ICRS with a 1 and 2 px shift. - """ - pipe = [(detector_2d, model_2d_shift), (icrs_sky_frame, None)] - - return wcs.WCS(pipe) + return examples.gwcs_2d_spatial_shift() @pytest.fixture def gwcs_2d_spatial_reordered(): - """ - A simple one step spatial WCS, in ICRS with a 1 and 2 px shift. - """ - out_frame = cf.CelestialFrame(reference_frame=coord.ICRS(), - axes_order=(1, 0)) - return wcs.WCS(model_2d_shift | models.Mapping((1, 0)), input_frame=detector_2d, output_frame=out_frame) + return examples.gwcs_2d_spatial_reordered() @pytest.fixture def gwcs_1d_freq(): - return wcs.WCS([(detector_1d, model_1d_scale), (freq_frame, None)]) + return examples.gwcs_1d_freq() @pytest.fixture def gwcs_3d_spatial_wave(): - comp1 = cf.CompositeFrame([icrs_sky_frame, wave_frame]) - m = model_2d_shift & model_1d_scale - - detector_frame = cf.CoordinateFrame(name="detector", naxes=3, - axes_order=(0, 1, 2), - axes_type=("pixel", "pixel", "pixel"), - axes_names=("x", "y", "z"), unit=(u.pix, u.pix, u.pix)) - - return wcs.WCS([(detector_frame, m), - (comp1, None)]) - + return examples.gwcs_3d_spatial_wave() @pytest.fixture def gwcs_2d_shift_scale(): - m1 = models.Shift(1) & models.Shift(2) - m2 = models.Scale(5) & models.Scale(10) - m3 = m1 | m2 - pipe = [(detector_2d, m3), (icrs_sky_frame, None)] - return wcs.WCS(pipe) + return examples.gwcs_2d_shift_scale() @pytest.fixture def gwcs_1d_freq_quantity(): - - detector_1d = cf.CoordinateFrame(name='detector', axes_order=(0,), naxes=1, unit=u.pix, axes_type="detector") - return wcs.WCS([(detector_1d, models.Multiply(1 * u.Hz / u.pix)), (freq_frame, None)]) - + return examples.gwcs_1d_freq_quantity() @pytest.fixture def gwcs_2d_shift_scale_quantity(): - m4 = models.Shift(1 * u.pix) & models.Shift(2 * u.pix) - m5 = models.Scale(5 * u.deg) - m6 = models.Scale(10 * u.deg) - m5.input_units_equivalencies = {'x': u.pixel_scale(1 * u.deg / u.pix)} - m6.input_units_equivalencies = {'x': u.pixel_scale(1 * u.deg / u.pix)} - m5.inverse = models.Scale(1. / 5 * u.pix) - m6.inverse = models.Scale(1. / 10 * u.pix) - m5.inverse.input_units_equivalencies = { - 'x': u.pixel_scale(1 * u.pix / u.deg) - } - m6.inverse.input_units_equivalencies = { - 'x': u.pixel_scale(1 * u.pix / u.deg) - } - m7 = m5 & m6 - m8 = m4 | m7 - pipe2 = [(detector_2d, m8), (icrs_sky_frame, None)] - return wcs.WCS(pipe2) + return examples.gwcs_2d_shift_scale_quantity() @pytest.fixture def gwcs_3d_identity_units(): - """ - A simple 1-1 gwcs that converts from pixels to arcseconds - """ - identity = (models.Multiply(1 * u.arcsec / u.pixel) & - models.Multiply(1 * u.arcsec / u.pixel) & - models.Multiply(1 * u.nm / u.pixel)) - sky_frame = cf.CelestialFrame(axes_order=(0, 1), name='icrs', - reference_frame=coord.ICRS(), - axes_names=("longitude", "latitude")) - wave_frame = cf.SpectralFrame(axes_order=(2, ), unit=u.nm, axes_names=("wavelength",)) - - frame = cf.CompositeFrame([sky_frame, wave_frame]) - - detector_frame = cf.CoordinateFrame(name="detector", naxes=3, - axes_order=(0, 1, 2), - axes_type=("pixel", "pixel", "pixel"), - axes_names=("x", "y", "z"), unit=(u.pix, u.pix, u.pix)) - - return wcs.WCS(forward_transform=identity, output_frame=frame, input_frame=detector_frame) + return examples.gwcs_3d_identity_units() @pytest.fixture def gwcs_4d_identity_units(): - """ - A simple 1-1 gwcs that converts from pixels to arcseconds - """ - identity = (models.Multiply(1*u.arcsec/u.pixel) & models.Multiply(1*u.arcsec/u.pixel) & - models.Multiply(1*u.nm/u.pixel) & models.Multiply(1*u.s/u.pixel)) - sky_frame = cf.CelestialFrame(axes_order=(0, 1), name='icrs', - reference_frame=coord.ICRS()) - wave_frame = cf.SpectralFrame(axes_order=(2, ), unit=u.nm) - time_frame = cf.TemporalFrame(axes_order=(3, ), unit=u.s, - reference_frame=Time("2000-01-01T00:00:00")) - - frame = cf.CompositeFrame([sky_frame, wave_frame, time_frame]) - - detector_frame = cf.CoordinateFrame(name="detector", naxes=4, - axes_order=(0, 1, 2, 3), - axes_type=("pixel", "pixel", "pixel", "pixel"), - axes_names=("x", "y", "z", "s"), unit=(u.pix, u.pix, u.pix, u.pix)) - - return wcs.WCS(forward_transform=identity, output_frame=frame, input_frame=detector_frame) + return examples.gwcs_4d_identity_units() @pytest.fixture def gwcs_simple_imaging_units(): - shift_by_crpix = models.Shift(-2048*u.pix) & models.Shift(-1024*u.pix) - matrix = np.array([[1.290551569736E-05, 5.9525007864732E-06], - [5.0226382102765E-06 , -1.2644844123757E-05]]) - rotation = models.AffineTransformation2D(matrix * u.deg, - translation=[0, 0] * u.deg) - rotation.input_units_equivalencies = {"x": u.pixel_scale(1*u.deg/u.pix), - "y": u.pixel_scale(1*u.deg/u.pix)} - rotation.inverse = models.AffineTransformation2D(np.linalg.inv(matrix) * u.pix, - translation=[0, 0] * u.pix) - rotation.inverse.input_units_equivalencies = {"x": u.pixel_scale(1*u.pix/u.deg), - "y": u.pixel_scale(1*u.pix/u.deg)} - tan = models.Pix2Sky_TAN() - celestial_rotation = models.RotateNative2Celestial(5.63056810618*u.deg, - -72.05457184279*u.deg, - 180*u.deg) - det2sky = shift_by_crpix | rotation | tan | celestial_rotation - det2sky.name = "linear_transform" - - detector_frame = cf.Frame2D(name="detector", axes_names=("x", "y"), - unit=(u.pix, u.pix)) - sky_frame = cf.CelestialFrame(reference_frame=coord.ICRS(), name='icrs', - unit=(u.deg, u.deg)) - pipeline = [(detector_frame, det2sky), - (sky_frame, None) - ] - return wcs.WCS(pipeline) + return examples.gwcs_simple_imaging_units() @pytest.fixture -def gwcs_stokes_lookup(): - transform = models.Tabular1D([0, 1, 2, 3] * u.pix, [1, 2, 3, 4] * u.one, - method="nearest", fill_value=np.nan, bounds_error=False) - frame = cf.StokesFrame() +def gwcs_simple_imaging(): + return examples.gwcs_simple_imaging() - detector_frame = cf.CoordinateFrame(name="detector", naxes=1, - axes_order=(0,), - axes_type=("pixel",), - axes_names=("x",), unit=(u.pix,)) - return wcs.WCS(forward_transform=transform, output_frame=frame, input_frame=detector_frame) +@pytest.fixture +def gwcs_stokes_lookup(): + return examples.gwcs_stokes_lookup() @pytest.fixture def gwcs_3spectral_orders(): - comp1 = cf.CompositeFrame([icrs_sky_frame, wave_frame]) - detector_frame = cf.Frame2D(name="detector", axes_names=("x", "y"), - unit=(u.pix, u.pix)) - m = model_2d_shift & model_1d_scale - - return wcs.WCS([(detector_frame, m), - (comp1, None)]) + return examples.gwcs_3spectral_orders() @pytest.fixture def gwcs_with_frames_strings(): - transform = models.Shift(1) & models.Shift(1) & models.Polynomial2D(1) - pipe = [('detector', transform), - ('world', None) - ] - return wcs.WCS(pipe) + return examples.gwcs_with_frames_strings() @pytest.fixture def sellmeier_glass(): - B_coef = [0.58339748, 0.46085267, 3.8915394] - C_coef = [0.00252643, 0.010078333, 1200.556] - return sp.SellmeierGlass(B_coef, C_coef) + return examples.sellmeier_glass() @pytest.fixture def sellmeier_zemax(): - B_coef = [0.58339748, 0.46085267, 3.8915394] - C_coef = [0.00252643, 0.010078333, 1200.556] - D_coef = [-2.66e-05, 0.0, 0.0] - E_coef = [0., 0., 0.] - return sp.SellmeierZemax(65, 35, 0, 0, B_coef = B_coef, - C_coef=C_coef, D_coef=D_coef, - E_coef=E_coef) + return examples.sellmeier_zemax() @pytest.fixture(scope="function") def gwcs_3d_galactic_spectral(): - """ - This fixture has the axes ordered as lat, spectral, lon. - """ - # lat,wav,lon - crpix1, crpix2, crpix3 = 29, 39, 44 - crval1, crval2, crval3 = 10, 20, 25 - cdelt1, cdelt2, cdelt3 = -0.01, 0.5, 0.01 - - shift = models.Shift(-crpix3) & models.Shift(-crpix1) - scale = models.Multiply(cdelt3) & models.Multiply(cdelt1) - proj = models.Pix2Sky_CAR() - skyrot = models.RotateNative2Celestial(crval3, 90 + crval1, 180) - celestial = shift | scale | proj | skyrot - - wave_model = models.Shift(-crpix2) | models.Multiply(cdelt2) | models.Shift(crval2) - - transform = models.Mapping((2, 0, 1)) | celestial & wave_model | models.Mapping((1, 2, 0)) - transform.bounding_box = ((5, 50), (-2, 45), (-1, 35)) - - sky_frame = cf.CelestialFrame(axes_order=(2, 0), - reference_frame=coord.Galactic(), axes_names=("Longitude", "Latitude")) - wave_frame = cf.SpectralFrame(axes_order=(1, ), unit=u.Hz, axes_names=("Frequency",)) - - frame = cf.CompositeFrame([sky_frame, wave_frame]) - - detector_frame = cf.CoordinateFrame(name="detector", naxes=3, - axes_order=(0, 1, 2), - axes_type=("pixel", "pixel", "pixel"), - unit=(u.pix, u.pix, u.pix)) - - owcs = wcs.WCS(forward_transform=transform, output_frame=frame, input_frame=detector_frame) - owcs.array_shape = (30, 20, 10) - owcs.pixel_shape = (10, 20, 30) - - return owcs + return examples.gwcs_3d_galactic_spectral() @pytest.fixture(scope="function") def gwcs_1d_spectral(): - """ - A simple 1D spectral WCS. - """ - wave_model = models.Shift(-5) | models.Multiply(3.7) | models.Shift(20) - wave_model.bounding_box = (7, 50) - wave_frame = cf.SpectralFrame(axes_order=(0, ), unit=u.Hz, axes_names=("Frequency",)) - - detector_frame = cf.CoordinateFrame( - name="detector", naxes=1, axes_order=(0, ), - axes_type=("pixel",), unit=(u.pix, ) - ) - - owcs = wcs.WCS(forward_transform=wave_model, output_frame=wave_frame, input_frame=detector_frame) - owcs.array_shape = (44, ) - owcs.pixel_shape = (44, ) - - return owcs + return examples.gwcs_1d_spectral() @pytest.fixture(scope="function") def gwcs_spec_cel_time_4d(): - """ - A complex 4D mixed celestial + spectral + time WCS. - """ - # spectroscopic frame: - wave_model = models.Shift(-5) | models.Multiply(3.7) | models.Shift(20) - wave_model.bounding_box = (7, 50) - wave_frame = cf.SpectralFrame(name='wave', unit=u.m, axes_order=(0,), axes_names=('lambda',)) - - # time frame: - time_model = models.Identity(1) # models.Linear1D(10, 0) - time_frame = cf.TemporalFrame(Time("2010-01-01T00:00"), name='time', unit=u.s, axes_order=(3,)) - - # Values from data/acs.hdr: - crpix = (12, 13) - crval = (5.63, -72.05) - cd = [[1.291E-05, 5.9532E-06], [5.02215E-06, -1.2645E-05]] - aff = models.AffineTransformation2D(matrix=cd, name='rotation') - offx = models.Shift(-crpix[0], name='x_translation') - offy = models.Shift(-crpix[1], name='y_translation') - wcslin = models.Mapping((1, 0)) | (offx & offy) | aff - tan = models.Pix2Sky_TAN(name='tangent_projection') - n2c = models.RotateNative2Celestial(*crval, 180, name='sky_rotation') - cel_model = wcslin | tan | n2c - icrs = cf.CelestialFrame(reference_frame=coord.ICRS(), name='sky', axes_order=(2, 1)) - - wcs_forward = wave_model & cel_model & time_model - - comp_frm = cf.CompositeFrame(frames=[wave_frame, icrs, time_frame], name='TEST 4D FRAME') - - detector_frame = cf.CoordinateFrame( - name="detector", naxes=4, axes_order=(0, 1, 2, 3), - axes_type=("pixel", "pixel", "pixel", "pixel"), - unit=(u.pix, u.pix, u.pix, u.pix) - ) - - w = wcs.WCS(forward_transform=wcs_forward, output_frame=comp_frm, input_frame=detector_frame) - - w.bounding_box = ((0, 63), (0, 127), (0, 255), (0, 9)) - w.array_shape = (10, 256, 128, 64) - w.pixel_shape = (64, 128, 256, 10) - return w + return examples.gwcs_spec_cel_time_4d() @pytest.fixture( @@ -357,49 +115,8 @@ def gwcs_spec_cel_time_4d(): ] ) def gwcs_cube_with_separable_spectral(request): - cube_size = (128, 64, 100) - axes_order = request.param - spectral_axes_order = (axes_order.index(2), ) - cel_axes_order = (axes_order.index(0), axes_order.index(1)) - - # Values from data/acs.hdr: - crpix = (64, 32) - crval = (5.63056810618, -72.0545718428) - cd = [[1.29058667557984E-05, 5.95320245884555E-06], - [5.02215195623825E-06, -1.2645010396976E-05]] - - aff = models.AffineTransformation2D(matrix=cd, name='rotation') - offx = models.Shift(-crpix[0], name='x_translation') - offy = models.Shift(-crpix[1], name='y_translation') - - wcslin = (offx & offy) | aff - tan = models.Pix2Sky_TAN(name='tangent_projection') - n2c = models.RotateNative2Celestial(*crval, 180, name='sky_rotation') - icrs = cf.CelestialFrame(reference_frame=coord.ICRS(), name='sky', - axes_order=cel_axes_order) - spec = cf.SpectralFrame( - name='wave', unit=[u.m,], axes_order=spectral_axes_order, - axes_names=('lambda',) - ) - comp_frm = cf.CompositeFrame(frames=[icrs, spec], name='TEST 3D FRAME WITH SPECTRAL AXIS') - wcs_forward = ((wcslin & models.Identity(1)) | - (tan & models.Identity(1)) | - (n2c & models.Identity(1)) | - models.Mapping(axes_order)) - - detector_frame = cf.CoordinateFrame(name="detector", naxes=3, - axes_order=(0, 1, 2), - axes_type=("pixel", "pixel", "pixel"), - unit=(u.pix, u.pix, u.pix)) - - w = wcs.WCS(forward_transform=wcs_forward, output_frame=comp_frm, - input_frame=detector_frame) - w.bounding_box = tuple((0, k - 1) for k in cube_size) - w.pixel_shape = cube_size - w.array_shape = w.pixel_shape[::-1] - - return w, axes_order + return examples.gwcs_cube_with_separable_spectral(axes_order) @pytest.fixture( @@ -412,106 +129,13 @@ def gwcs_cube_with_separable_spectral(request): ] ) def gwcs_cube_with_separable_time(request): - """ - A mixed celestial + time WCS. - """ - cube_size = (64, 32, 128) - axes_order = request.param - time_axes_order = (axes_order.index(2), ) - cel_axes_order = (axes_order.index(0), axes_order.index(1)) - - detector_frame = cf.CoordinateFrame( - name="detector", naxes=3, axes_order=(0, 1, 2), - axes_type=("pixel", "pixel", "pixel"), unit=(u.pix, u.pix, u.pix) - ) - - # time frame: - time_model = models.Identity(1) # models.Linear1D(10, 0) - time_frame = cf.TemporalFrame(Time("2010-01-01T00:00"), name='time', - unit=u.s, axes_order=time_axes_order) - - # Values from data/acs.hdr: - crpix = (12, 13) - crval = (5.63, -72.05) - cd = [[1.291E-05, 5.9532E-06], [5.02215E-06, -1.2645E-05]] - aff = models.AffineTransformation2D(matrix=cd, name='rotation') - offx = models.Shift(-crpix[0], name='x_translation') - offy = models.Shift(-crpix[1], name='y_translation') - wcslin = models.Mapping((1, 0)) | (offx & offy) | aff - tan = models.Pix2Sky_TAN(name='tangent_projection') - n2c = models.RotateNative2Celestial(*crval, 180, name='sky_rotation') - cel_model = wcslin | tan | n2c - icrs = cf.CelestialFrame(reference_frame=coord.ICRS(), name='sky', - axes_order=cel_axes_order) - - wcs_forward = (cel_model & time_model) | models.Mapping(axes_order) - - comp_frm = cf.CompositeFrame(frames=[icrs, time_frame], - name='TEST 3D FRAME WITH TIME') - - w = wcs.WCS(forward_transform=wcs_forward, output_frame=comp_frm, - input_frame=detector_frame) - - w.bounding_box = tuple((0, k - 1) for k in cube_size) - w.pixel_shape = cube_size - w.array_shape = w.pixel_shape[::-1] - return w + return examples.gwcs_cube_with_separable_time(axes_order) @pytest.fixture(scope="function") def gwcs_7d_complex_mapping(): - """ - Useful features of this WCS (axes indices here are 0-based): - - includes two celestial axes: input (0, 1) maps to world (2 - RA, 1 - Dec) - - includes one separable frame with one axis: 4 -> 2 - - includes one frame with 3 input and 4 output axes (1 degenerate), - with separable world axes (3, 5) and (0, 6). - """ - offx = models.Shift(-64, name='x_translation') - offy = models.Shift(-32, name='y_translation') - cd = np.array([[1.2906, 0.59532], [0.50222, -1.2645]]) - aff = models.AffineTransformation2D(matrix=1e-5 * cd, name='rotation') - aff2 = models.AffineTransformation2D(matrix=cd, name='rotation2') - - wcslin = (offx & offy) | aff - tan = models.Pix2Sky_TAN(name='tangent_projection') - n2c = models.RotateNative2Celestial(5.630568, -72.0546, 180, name='skyrot') - icrs = cf.CelestialFrame(reference_frame=coord.ICRS(), name='sky', axes_order=(2, 1)) - spec = cf.SpectralFrame(name='wave', unit=[u.m], axes_order=(4,), axes_names=('lambda',)) - cmplx = cf.CoordinateFrame( - name="complex", - naxes=4, - axes_order=(3, 5, 0, 6), - axis_physical_types=(['em.wl', 'em.wl', 'time', 'time']), - axes_type=("SPATIAL", "SPATIAL", "TIME", "TIME"), - axes_names=("x", "y", "t", 'tau'), - unit=(u.m, u.m, u.second, u.second) - ) - - comp_frm = cf.CompositeFrame(frames=[icrs, spec, cmplx], name='TEST 7D') - wcs_forward = ((wcslin & models.Shift(-3.14) & models.Scale(2.7) & aff2) | - (tan & models.Identity(1) & models.Identity(1) & models.Identity(2)) | - (n2c & models.Identity(1) & models.Identity(1) & models.Identity(2)) | - models.Mapping((3, 1, 0, 4, 2, 5, 3))) - - - detector_frame = cf.CoordinateFrame( - name="detector", naxes=6, - axes_order=(0, 1, 2, 3, 4, 5), - axes_type=("pixel", "pixel", "pixel", "pixel", "pixel", "pixel"), - unit=(u.pix, u.pix, u.pix, u.pix, u.pix, u.pix) - ) - - pipeline = [('detector', wcs_forward), (comp_frm, None)] - w = wcs.WCS(forward_transform=wcs_forward, output_frame=comp_frm, - input_frame=detector_frame) - w.bounding_box = ((0, 15), (0, 31), (0, 20), (0, 10), (0, 10), (0, 1)) - - w.array_shape = (2, 11, 11, 21, 32, 16) - w.pixel_shape = (16, 32, 21, 11, 11, 2) - - return w + return examples.gwcs_7d_complex_mapping() @pytest.fixture diff --git a/gwcs/tests/data/miri_lrs_wcs.asdf b/gwcs/tests/data/miri_lrs_wcs.asdf index 48d0e44d..e5fe5bbd 100644 --- a/gwcs/tests/data/miri_lrs_wcs.asdf +++ b/gwcs/tests/data/miri_lrs_wcs.asdf @@ -3,97 +3,57 @@ %YAML 1.1 %TAG ! tag:stsci.edu:asdf/ --- !core/asdf-1.1.0 -asdf_library: !core/software-1.0.0 {author: Space Telescope Science Institute, homepage: 'http://github.com/spacetelescope/asdf', - name: asdf, version: 2.6.0} +asdf_library: !core/software-1.0.0 {author: The ASDF Developers, homepage: 'http://github.com/asdf-format/asdf', + name: asdf, version: 3.5.0} history: extensions: - !core/extension_metadata-1.0.0 - extension_class: asdf.extension.BuiltinExtension - software: !core/software-1.0.0 {name: asdf, version: 2.6.0} + extension_class: asdf.extension._manifest.ManifestExtension + extension_uri: asdf://asdf-format.org/core/extensions/core-1.5.0 + manifest_software: !core/software-1.0.0 {name: asdf_standard, version: 1.1.1} + software: !core/software-1.0.0 {name: asdf, version: 3.5.0} - !core/extension_metadata-1.0.0 - extension_class: gwcs.extension.GWCSExtension - software: !core/software-1.0.0 {name: gwcs, version: 0.13.1.dev16+g71f9b60} -wcs: ! + extension_class: asdf.extension._manifest.ManifestExtension + extension_uri: asdf://astropy.org/astropy/extensions/units-1.0.0 + software: !core/software-1.0.0 {name: asdf-astropy, version: 0.6.1} + - !core/extension_metadata-1.0.0 + extension_class: asdf.extension._manifest.ManifestExtension + extension_uri: asdf://asdf-format.org/astronomy/gwcs/extensions/gwcs-1.2.0 + manifest_software: !core/software-1.0.0 {name: asdf_wcs_schemas, version: 0.4.0} + software: !core/software-1.0.0 {name: gwcs, version: 0.22.0a1.dev14+gc46e932} + - !core/extension_metadata-1.0.0 + extension_class: asdf.extension._manifest.ManifestExtension + extension_uri: asdf://asdf-format.org/astronomy/coordinates/extensions/coordinates-1.0.0 + manifest_software: !core/software-1.0.0 {name: asdf_coordinates_schemas, version: 0.3.0} + software: !core/software-1.0.0 {name: asdf-astropy, version: 0.6.1} + - !core/extension_metadata-1.0.0 + extension_class: asdf.extension._manifest.ManifestExtension + extension_uri: asdf://asdf-format.org/transform/extensions/transform-1.5.0 + manifest_software: !core/software-1.0.0 {name: asdf_transform_schemas, version: 0.5.0} + software: !core/software-1.0.0 {name: asdf-astropy, version: 0.6.1} +wcs: ! name: '' + pixel_shape: null steps: - - ! + - ! frame: ! axes_names: [x, y] axes_order: [0, 1] axis_physical_types: ['custom:x', 'custom:y'] name: detector - unit: [!unit/unit-1.0.0 'pixel', !unit/unit-1.0.0 'pixel'] + unit: [!unit/unit-1.0.0 pixel, !unit/unit-1.0.0 pixel] transform: !transform/compose-1.2.0 - bounding_box: - - [6.5, 396.5] - - [302.5, 346.5] - bounds: - amplitude_7: [null, null] - angle_22: [null, null] - angle_3: [null, null] - c0_0_13: [null, null] - c0_0_14: [null, null] - c0_0_16: [null, null] - c0_0_17: [null, null] - c0_10: [null, null] - c0_11: [null, null] - c0_1_13: [null, null] - c0_1_14: [null, null] - c0_1_16: [null, null] - c0_1_17: [null, null] - c0_2_13: [null, null] - c0_2_14: [null, null] - c0_3_13: [null, null] - c0_3_14: [null, null] - c0_4_13: [null, null] - c0_4_14: [null, null] - c1_0_13: [null, null] - c1_0_14: [null, null] - c1_0_16: [null, null] - c1_0_17: [null, null] - c1_10: [null, null] - c1_11: [null, null] - c1_1_13: [null, null] - c1_1_14: [null, null] - c1_2_13: [null, null] - c1_2_14: [null, null] - c1_3_13: [null, null] - c1_3_14: [null, null] - c2_0_13: [null, null] - c2_0_14: [null, null] - c2_1_13: [null, null] - c2_1_14: [null, null] - c2_2_13: [null, null] - c2_2_14: [null, null] - c3_0_13: [null, null] - c3_0_14: [null, null] - c3_1_13: [null, null] - c3_1_14: [null, null] - c4_0_13: [null, null] - c4_0_14: [null, null] - offset_1: [null, null] - offset_2: [null, null] - offset_20: [null, null] - offset_21: [null, null] - offset_4: [null, null] - offset_5: [null, null] - offset_8: [null, null] - fixed: {amplitude_7: false, angle_22: false, angle_3: false, c0_0_13: false, - c0_0_14: false, c0_0_16: false, c0_0_17: false, c0_10: false, c0_11: false, - c0_1_13: false, c0_1_14: false, c0_1_16: false, c0_1_17: false, c0_2_13: false, - c0_2_14: false, c0_3_13: false, c0_3_14: false, c0_4_13: false, c0_4_14: false, - c1_0_13: false, c1_0_14: false, c1_0_16: false, c1_0_17: false, c1_10: false, - c1_11: false, c1_1_13: false, c1_1_14: false, c1_2_13: false, c1_2_14: false, - c1_3_13: false, c1_3_14: false, c2_0_13: false, c2_0_14: false, c2_1_13: false, - c2_1_14: false, c2_2_13: false, c2_2_14: false, c3_0_13: false, c3_0_14: false, - c3_1_13: false, c3_1_14: false, c4_0_13: false, c4_0_14: false, offset_1: false, - offset_2: false, offset_20: false, offset_21: false, offset_4: false, offset_5: false, - offset_8: false} + bounding_box: !transform/property/bounding_box-1.0.0 + ignore: [] + intervals: + x0: [302.5, 346.5] + x1: [6.5, 396.5] + order: F forward: - - !transform/remap_axes-1.2.0 - bounds: {} - fixed: {} + - !transform/remap_axes-1.3.0 + inputs: [x0, x1] mapping: [0, 1, 0, 1] + outputs: [x0, x1, x2, x3] - !transform/concatenate-1.2.0 forward: - !transform/compose-1.2.0 @@ -103,51 +63,54 @@ wcs: ! - !transform/concatenate-1.2.0 forward: - !transform/shift-1.2.0 - bounds: - offset: &id001 [null, null] - fixed: {offset: false} + inputs: [x] offset: -325.13 + outputs: [y] - !transform/shift-1.2.0 - bounds: - offset: *id001 - fixed: {offset: false} + inputs: [x] offset: -299.7 - - !transform/rotate2d-1.2.0 + outputs: [y] + inputs: [x0, x1] + outputs: [y0, y1] + - !transform/rotate2d-1.3.0 angle: 0.24155757363306068 - bounds: - angle: &id002 [null, null] - fixed: {angle: false} + inputs: [x, y] + outputs: [x, y] + inputs: [x0, x1] + outputs: [x, y] - !transform/compose-1.2.0 forward: - !transform/concatenate-1.2.0 forward: - !transform/shift-1.2.0 - bounds: - offset: *id001 - fixed: {offset: false} + inputs: [x] offset: 325.13 + outputs: [y] - !transform/shift-1.2.0 - bounds: - offset: *id001 - fixed: {offset: false} + inputs: [x] offset: 299.7 + outputs: [y] + inputs: [x0, x1] + outputs: [y0, y1] - !transform/compose-1.2.0 forward: - !transform/concatenate-1.2.0 forward: - !transform/identity-1.2.0 - bounds: {} - fixed: {} - - !transform/constant-1.2.0 - bounds: - amplitude: [null, null] - fixed: {amplitude: false} - inverse: !transform/constant-1.2.0 - bounds: - amplitude: [null, null] - fixed: {amplitude: false} + inputs: [x0] + outputs: [x0] + - !transform/constant-1.4.0 + dimensions: 1 + inputs: [x] + inverse: !transform/constant-1.4.0 + dimensions: 1 + inputs: [x] + outputs: [y] value: 299.7 + outputs: [y] value: 299.7 + inputs: [x0, x] + outputs: [x0, y] - !transform/compose-1.2.0 forward: - !transform/compose-1.2.0 @@ -165,92 +128,75 @@ wcs: ! - !transform/concatenate-1.2.0 forward: - !transform/shift-1.2.0 - bounds: - offset: *id001 - fixed: {offset: false} + inputs: [x] offset: -4.0 + outputs: [y] - !transform/identity-1.2.0 - bounds: {} - fixed: {} + inputs: [x0] + outputs: [x0] + inputs: [x, x0] + outputs: [y, x0] - !transform/concatenate-1.2.0 forward: - !transform/polynomial-1.2.0 - bounds: - c0: [null, null] - c1: [null, null] coefficients: !core/ndarray-1.0.0 source: 0 datatype: float64 byteorder: little shape: [2] domain: [-1, 1] - fixed: {c0: false, c1: false} + inputs: [x] inverse: !transform/polynomial-1.2.0 - bounds: - c0: [null, null] - c1: [null, null] coefficients: !core/ndarray-1.0.0 source: 1 datatype: float64 byteorder: little shape: [2] domain: [-1, 1] - fixed: {c0: false, c1: false} + inputs: [x] + outputs: [y] window: [-1, 1] name: M_column_correction + outputs: [y] window: [-1, 1] - !transform/polynomial-1.2.0 - bounds: - c0: [null, null] - c1: [null, null] coefficients: !core/ndarray-1.0.0 source: 2 datatype: float64 byteorder: little shape: [2] domain: [-1, 1] - fixed: {c0: false, c1: false} + inputs: [x] inverse: !transform/polynomial-1.2.0 - bounds: - c0: [null, null] - c1: [null, null] coefficients: !core/ndarray-1.0.0 source: 3 datatype: float64 byteorder: little shape: [2] domain: [-1, 1] - fixed: {c0: false, c1: false} + inputs: [x] + outputs: [y] window: [-1, 1] name: M_row_correction + outputs: [y] window: [-1, 1] - - !transform/remap_axes-1.2.0 - bounds: {} - fixed: {} + inputs: [x0, x1] + outputs: [y0, y1] + inputs: [x, x0] + outputs: [y0, y1] + - !transform/remap_axes-1.3.0 + inputs: [x0, x1] inverse: !transform/identity-1.2.0 - bounds: {} - fixed: {} + inputs: [x0, x1] n_dims: 2 + outputs: [x0, x1] mapping: [0, 1, 0, 1] + outputs: [x0, x1, x2, x3] + inputs: [x, x0] + outputs: [x0, x1, x2, x3] - !transform/concatenate-1.2.0 forward: - !transform/polynomial-1.2.0 - bounds: - c0_0: [null, null] - c0_1: [null, null] - c0_2: [null, null] - c0_3: [null, null] - c0_4: [null, null] - c1_0: [null, null] - c1_1: [null, null] - c1_2: [null, null] - c1_3: [null, null] - c2_0: [null, null] - c2_1: [null, null] - c2_2: [null, null] - c3_0: [null, null] - c3_1: [null, null] - c4_0: [null, null] coefficients: !core/ndarray-1.0.0 source: 4 datatype: float64 @@ -259,27 +205,8 @@ wcs: ! domain: - [-1, 1] - [-1, 1] - fixed: {c0_0: false, c0_1: false, c0_2: false, c0_3: false, - c0_4: false, c1_0: false, c1_1: false, c1_2: false, - c1_3: false, c2_0: false, c2_1: false, c2_2: false, - c3_0: false, c3_1: false, c4_0: false} + inputs: [x, y] inverse: !transform/polynomial-1.2.0 - bounds: - c0_0: [null, null] - c0_1: [null, null] - c0_2: [null, null] - c0_3: [null, null] - c0_4: [null, null] - c1_0: [null, null] - c1_1: [null, null] - c1_2: [null, null] - c1_3: [null, null] - c2_0: [null, null] - c2_1: [null, null] - c2_2: [null, null] - c3_0: [null, null] - c3_1: [null, null] - c4_0: [null, null] coefficients: !core/ndarray-1.0.0 source: 5 datatype: float64 @@ -288,34 +215,17 @@ wcs: ! domain: - [-1, 1] - [-1, 1] - fixed: {c0_0: false, c0_1: false, c0_2: false, c0_3: false, - c0_4: false, c1_0: false, c1_1: false, c1_2: false, - c1_3: false, c2_0: false, c2_1: false, c2_2: false, - c3_0: false, c3_1: false, c4_0: false} + inputs: [x, y] + outputs: [z] window: - [-1, 1] - [-1, 1] name: B_correction + outputs: [z] window: - [-1, 1] - [-1, 1] - !transform/polynomial-1.2.0 - bounds: - c0_0: [null, null] - c0_1: [null, null] - c0_2: [null, null] - c0_3: [null, null] - c0_4: [null, null] - c1_0: [null, null] - c1_1: [null, null] - c1_2: [null, null] - c1_3: [null, null] - c2_0: [null, null] - c2_1: [null, null] - c2_2: [null, null] - c3_0: [null, null] - c3_1: [null, null] - c4_0: [null, null] coefficients: !core/ndarray-1.0.0 source: 6 datatype: float64 @@ -324,27 +234,8 @@ wcs: ! domain: - [-1, 1] - [-1, 1] - fixed: {c0_0: false, c0_1: false, c0_2: false, c0_3: false, - c0_4: false, c1_0: false, c1_1: false, c1_2: false, - c1_3: false, c2_0: false, c2_1: false, c2_2: false, - c3_0: false, c3_1: false, c4_0: false} + inputs: [x, y] inverse: !transform/polynomial-1.2.0 - bounds: - c0_0: [null, null] - c0_1: [null, null] - c0_2: [null, null] - c0_3: [null, null] - c0_4: [null, null] - c1_0: [null, null] - c1_1: [null, null] - c1_2: [null, null] - c1_3: [null, null] - c2_0: [null, null] - c2_1: [null, null] - c2_2: [null, null] - c3_0: [null, null] - c3_1: [null, null] - c4_0: [null, null] coefficients: !core/ndarray-1.0.0 source: 7 datatype: float64 @@ -353,32 +244,33 @@ wcs: ! domain: - [-1, 1] - [-1, 1] - fixed: {c0_0: false, c0_1: false, c0_2: false, c0_3: false, - c0_4: false, c1_0: false, c1_1: false, c1_2: false, - c1_3: false, c2_0: false, c2_1: false, c2_2: false, - c3_0: false, c3_1: false, c4_0: false} + inputs: [x, y] + outputs: [z] window: - [-1, 1] - [-1, 1] name: A_correction + outputs: [z] window: - [-1, 1] - [-1, 1] - - !transform/remap_axes-1.2.0 - bounds: {} - fixed: {} - inverse: !transform/remap_axes-1.2.0 - bounds: {} - fixed: {} + inputs: [x0, y0, x1, y1] + outputs: [z0, z1] + inputs: [x, x0] + outputs: [z0, z1] + - !transform/remap_axes-1.3.0 + inputs: [x0, x1] + inverse: !transform/remap_axes-1.3.0 + inputs: [x0, x1] mapping: [0, 1, 0, 1] + outputs: [x0, x1, x2, x3] mapping: [0, 1, 0, 1] + outputs: [x0, x1, x2, x3] + inputs: [x, x0] + outputs: [x0, x1, x2, x3] - !transform/concatenate-1.2.0 forward: - !transform/polynomial-1.2.0 - bounds: - c0_0: [null, null] - c0_1: [null, null] - c1_0: [null, null] coefficients: !core/ndarray-1.0.0 source: 8 datatype: float64 @@ -387,16 +279,13 @@ wcs: ! domain: - [-1, 1] - [-1, 1] - fixed: {c0_0: false, c0_1: false, c1_0: false} + inputs: [x, y] name: TI_row_correction + outputs: [z] window: - [-1, 1] - [-1, 1] - !transform/polynomial-1.2.0 - bounds: - c0_0: [null, null] - c0_1: [null, null] - c1_0: [null, null] coefficients: !core/ndarray-1.0.0 source: 9 datatype: float64 @@ -405,23 +294,38 @@ wcs: ! domain: - [-1, 1] - [-1, 1] - fixed: {c0_0: false, c0_1: false, c1_0: false} + inputs: [x, y] name: TI_column_correction + outputs: [z] window: - [-1, 1] - [-1, 1] + inputs: [x0, y0, x1, y1] + outputs: [z0, z1] + inputs: [x, x0] + outputs: [z0, z1] - !transform/identity-1.2.0 - bounds: {} - fixed: {} - inverse: !transform/remap_axes-1.2.0 - bounds: {} - fixed: {} + inputs: [x0, x1] + inverse: !transform/remap_axes-1.3.0 + inputs: [x0, x1] mapping: [0, 1, 0, 1] + outputs: [x0, x1, x2, x3] n_dims: 2 - - !transform/remap_axes-1.2.0 - bounds: {} - fixed: {} + outputs: [x0, x1] + inputs: [x, x0] + outputs: [x0, x1] + - !transform/remap_axes-1.3.0 + inputs: [x0, x1] mapping: [1, 0] + outputs: [x0, x1] + inputs: [x, x0] + outputs: [x0, x1] + inputs: [x0, x] + outputs: [x0, x1] + inputs: [x0, x1] + outputs: [x0, x1] + inputs: [x0, x1] + outputs: [x0, x1] - !transform/compose-1.2.0 forward: - !transform/compose-1.2.0 @@ -431,42 +335,53 @@ wcs: ! - !transform/concatenate-1.2.0 forward: - !transform/shift-1.2.0 - bounds: - offset: *id001 - fixed: {offset: false} + inputs: [x] offset: -325.13 + outputs: [y] - !transform/shift-1.2.0 - bounds: - offset: *id001 - fixed: {offset: false} + inputs: [x] offset: -299.7 - - !transform/rotate2d-1.2.0 + outputs: [y] + inputs: [x0, x1] + outputs: [y0, y1] + - !transform/rotate2d-1.3.0 angle: 0.24155757363306068 - bounds: - angle: *id002 - fixed: {angle: false} - - !transform/remap_axes-1.2.0 - bounds: {} - fixed: {} + inputs: [x, y] + outputs: [x, y] + inputs: [x0, x1] + outputs: [x, y] + - !transform/remap_axes-1.3.0 + inputs: [x0, x1] mapping: [1] + outputs: [x0] + inputs: [x0, x1] + outputs: [x0] - !transform/tabular-1.2.0 - bounding_box: [-291.60259000740217, 95.69787960086536] - bounds: {} + bounding_box: !transform/property/bounding_box-1.0.0 + ignore: [] + intervals: + x: [-291.60259000740217, 95.69787960086536] + order: C bounds_error: false fill_value: .nan - fixed: {} + inputs: [x] inverse: !transform/tabular-1.2.0 - bounding_box: [3.596040329703134, 14.387876373781149] - bounds: {} + bounding_box: !transform/property/bounding_box-1.0.0 + ignore: [] + intervals: + x: [3.596040329703134, 14.387876373781149] + order: C bounds_error: false fill_value: .nan - fixed: {} + inputs: [x] lookup_table: !core/ndarray-1.0.0 source: 12 datatype: float64 byteorder: little shape: [699] + method: linear name: waverefinv + outputs: [y] points: - !core/ndarray-1.0.0 source: 13 @@ -478,72 +393,21 @@ wcs: ! datatype: float64 byteorder: little shape: [388] + method: linear name: waveref + outputs: [y] points: - !core/ndarray-1.0.0 source: 11 datatype: float64 byteorder: little shape: [388] + inputs: [x0, x1] + outputs: [y] + inputs: [x00, x10, x01, x11] + outputs: [x0, x1, y] + inputs: [x0, x1] inverse: !transform/compose-1.2.0 - bounds: - angle_14: [null, null] - angle_17: [null, null] - c0_0_2: [null, null] - c0_0_3: [null, null] - c0_0_5: [null, null] - c0_0_6: [null, null] - c0_1_2: [null, null] - c0_1_3: [null, null] - c0_1_5: [null, null] - c0_1_6: [null, null] - c0_2_5: [null, null] - c0_2_6: [null, null] - c0_3_5: [null, null] - c0_3_6: [null, null] - c0_4_5: [null, null] - c0_4_6: [null, null] - c0_8: [null, null] - c0_9: [null, null] - c1_0_2: [null, null] - c1_0_3: [null, null] - c1_0_5: [null, null] - c1_0_6: [null, null] - c1_1_5: [null, null] - c1_1_6: [null, null] - c1_2_5: [null, null] - c1_2_6: [null, null] - c1_3_5: [null, null] - c1_3_6: [null, null] - c1_8: [null, null] - c1_9: [null, null] - c2_0_5: [null, null] - c2_0_6: [null, null] - c2_1_5: [null, null] - c2_1_6: [null, null] - c2_2_5: [null, null] - c2_2_6: [null, null] - c3_0_5: [null, null] - c3_0_6: [null, null] - c3_1_5: [null, null] - c3_1_6: [null, null] - c4_0_5: [null, null] - c4_0_6: [null, null] - offset_10: [null, null] - offset_12: [null, null] - offset_13: [null, null] - offset_18: [null, null] - offset_19: [null, null] - fixed: {angle_14: false, angle_17: false, c0_0_2: false, c0_0_3: false, c0_0_5: false, - c0_0_6: false, c0_1_2: false, c0_1_3: false, c0_1_5: false, c0_1_6: false, - c0_2_5: false, c0_2_6: false, c0_3_5: false, c0_3_6: false, c0_4_5: false, - c0_4_6: false, c0_8: false, c0_9: false, c1_0_2: false, c1_0_3: false, c1_0_5: false, - c1_0_6: false, c1_1_5: false, c1_1_6: false, c1_2_5: false, c1_2_6: false, - c1_3_5: false, c1_3_6: false, c1_8: false, c1_9: false, c2_0_5: false, c2_0_6: false, - c2_1_5: false, c2_1_6: false, c2_2_5: false, c2_2_6: false, c3_0_5: false, - c3_0_6: false, c3_1_5: false, c3_1_6: false, c4_0_5: false, c4_0_6: false, - offset_10: false, offset_12: false, offset_13: false, offset_18: false, - offset_19: false} forward: - !transform/concatenate-1.2.0 forward: @@ -553,25 +417,21 @@ wcs: ! forward: - !transform/compose-1.2.0 forward: - - !transform/remap_axes-1.2.0 - bounds: {} - fixed: {} + - !transform/remap_axes-1.3.0 + inputs: [x0, x1] mapping: [1, 0] + outputs: [x0, x1] - !transform/compose-1.2.0 forward: - - !transform/remap_axes-1.2.0 - bounds: {} - fixed: {} + - !transform/remap_axes-1.3.0 + inputs: [x0, x1] mapping: [0, 1, 0, 1] + outputs: [x0, x1, x2, x3] - !transform/compose-1.2.0 forward: - !transform/concatenate-1.2.0 forward: - !transform/polynomial-1.2.0 - bounds: - c0_0: [null, null] - c0_1: [null, null] - c1_0: [null, null] coefficients: !core/ndarray-1.0.0 source: 14 datatype: float64 @@ -580,16 +440,13 @@ wcs: ! domain: - [-1, 1] - [-1, 1] - fixed: {c0_0: false, c0_1: false, c1_0: false} + inputs: [x, y] name: T_row_correction + outputs: [z] window: - [-1, 1] - [-1, 1] - !transform/polynomial-1.2.0 - bounds: - c0_0: [null, null] - c0_1: [null, null] - c1_0: [null, null] coefficients: !core/ndarray-1.0.0 source: 15 datatype: float64 @@ -598,38 +455,25 @@ wcs: ! domain: - [-1, 1] - [-1, 1] - fixed: {c0_0: false, c0_1: false, c1_0: false} + inputs: [x, y] name: T_column_correction + outputs: [z] window: - [-1, 1] - [-1, 1] + inputs: [x0, y0, x1, y1] + outputs: [z0, z1] - !transform/compose-1.2.0 forward: - - !transform/remap_axes-1.2.0 - bounds: {} - fixed: {} + - !transform/remap_axes-1.3.0 + inputs: [x0, x1] mapping: [0, 1, 0, 1] + outputs: [x0, x1, x2, x3] - !transform/compose-1.2.0 forward: - !transform/concatenate-1.2.0 forward: - !transform/polynomial-1.2.0 - bounds: - c0_0: [null, null] - c0_1: [null, null] - c0_2: [null, null] - c0_3: [null, null] - c0_4: [null, null] - c1_0: [null, null] - c1_1: [null, null] - c1_2: [null, null] - c1_3: [null, null] - c2_0: [null, null] - c2_1: [null, null] - c2_2: [null, null] - c3_0: [null, null] - c3_1: [null, null] - c4_0: [null, null] coefficients: !core/ndarray-1.0.0 source: 16 datatype: float64 @@ -638,30 +482,12 @@ wcs: ! domain: - [-1, 1] - [-1, 1] - fixed: {c0_0: false, c0_1: false, c0_2: false, c0_3: false, - c0_4: false, c1_0: false, c1_1: false, c1_2: false, - c1_3: false, c2_0: false, c2_1: false, c2_2: false, - c3_0: false, c3_1: false, c4_0: false} + inputs: [x, y] + outputs: [z] window: - [-1, 1] - [-1, 1] - !transform/polynomial-1.2.0 - bounds: - c0_0: [null, null] - c0_1: [null, null] - c0_2: [null, null] - c0_3: [null, null] - c0_4: [null, null] - c1_0: [null, null] - c1_1: [null, null] - c1_2: [null, null] - c1_3: [null, null] - c2_0: [null, null] - c2_1: [null, null] - c2_2: [null, null] - c3_0: [null, null] - c3_1: [null, null] - c4_0: [null, null] coefficients: !core/ndarray-1.0.0 source: 17 datatype: float64 @@ -670,119 +496,148 @@ wcs: ! domain: - [-1, 1] - [-1, 1] - fixed: {c0_0: false, c0_1: false, c0_2: false, c0_3: false, - c0_4: false, c1_0: false, c1_1: false, c1_2: false, - c1_3: false, c2_0: false, c2_1: false, c2_2: false, - c3_0: false, c3_1: false, c4_0: false} + inputs: [x, y] + outputs: [z] window: - [-1, 1] - [-1, 1] + inputs: [x0, y0, x1, y1] + outputs: [z0, z1] - !transform/compose-1.2.0 forward: - !transform/identity-1.2.0 - bounds: {} - fixed: {} + inputs: [x0, x1] n_dims: 2 + outputs: [x0, x1] - !transform/compose-1.2.0 forward: - !transform/concatenate-1.2.0 forward: - !transform/polynomial-1.2.0 - bounds: - c0: [null, null] - c1: [null, null] coefficients: !core/ndarray-1.0.0 source: 18 datatype: float64 byteorder: little shape: [2] domain: [-1, 1] - fixed: {c0: false, c1: false} + inputs: [x] + outputs: [y] window: [-1, 1] - !transform/polynomial-1.2.0 - bounds: - c0: [null, null] - c1: [null, null] coefficients: !core/ndarray-1.0.0 source: 19 datatype: float64 byteorder: little shape: [2] domain: [-1, 1] - fixed: {c0: false, c1: false} + inputs: [x] + outputs: [y] window: [-1, 1] + inputs: [x0, x1] + outputs: [y0, y1] - !transform/concatenate-1.2.0 forward: - !transform/shift-1.2.0 - bounds: - offset: &id003 [null, null] - fixed: {offset: false} + inputs: [x] offset: 4.0 + outputs: [y] - !transform/identity-1.2.0 - bounds: {} - fixed: {} + inputs: [x0] + outputs: [x0] + inputs: [x, x0] + outputs: [y, x0] + inputs: [x0, x1] + outputs: [y, x0] + inputs: [x0, x1] + outputs: [y, x0] + inputs: [x0, y0, x1, y1] + outputs: [y, x0] + inputs: [x0, x1] + outputs: [y, x0] + inputs: [x0, y0, x1, y1] + outputs: [y, x0] + inputs: [x0, x1] + outputs: [y, x0] + inputs: [x0, x1] + outputs: [y, x0] - !transform/compose-1.2.0 forward: - !transform/concatenate-1.2.0 forward: - !transform/shift-1.2.0 - bounds: - offset: *id003 - fixed: {offset: false} + inputs: [x] offset: -325.13 + outputs: [y] - !transform/shift-1.2.0 - bounds: - offset: *id003 - fixed: {offset: false} + inputs: [x] offset: -299.7 - - !transform/rotate2d-1.2.0 + outputs: [y] + inputs: [x0, x1] + outputs: [y0, y1] + - !transform/rotate2d-1.3.0 angle: 0.24155757363306068 - bounds: - angle: &id004 [null, null] - fixed: {angle: false} - - !transform/remap_axes-1.2.0 - bounds: {} - fixed: {} + inputs: [x, y] + outputs: [x, y] + inputs: [x0, x1] + outputs: [x, y] + inputs: [x0, x1] + outputs: [x, y] + - !transform/remap_axes-1.3.0 + inputs: [x0, x1] mapping: [0] n_inputs: 2 + outputs: [x0] + inputs: [x0, x1] + outputs: [x0] - !transform/tabular-1.2.0 - bounding_box: [3.596040329703134, 14.387876373781149] - bounds: {} + bounding_box: !transform/property/bounding_box-1.0.0 + ignore: [] + intervals: + x: [3.596040329703134, 14.387876373781149] + order: C bounds_error: false fill_value: .nan - fixed: {} + inputs: [x] lookup_table: !core/ndarray-1.0.0 source: 20 datatype: float64 byteorder: little shape: [699] + method: linear name: waverefinv + outputs: [y] points: - !core/ndarray-1.0.0 source: 21 datatype: float64 byteorder: little shape: [699] + inputs: [x0, x1, x] + outputs: [x0, y] - !transform/compose-1.2.0 forward: - - !transform/rotate2d-1.2.0 + - !transform/rotate2d-1.3.0 angle: -0.24155757363306068 - bounds: - angle: *id004 - fixed: {angle: false} + inputs: [x, y] + outputs: [x, y] - !transform/concatenate-1.2.0 forward: - !transform/shift-1.2.0 - bounds: - offset: *id003 - fixed: {offset: false} + inputs: [x] offset: 325.13 + outputs: [y] - !transform/shift-1.2.0 - bounds: - offset: *id003 - fixed: {offset: false} + inputs: [x] offset: 299.7 - - ! + outputs: [y] + inputs: [x0, x1] + outputs: [y0, y1] + inputs: [x, y] + outputs: [y0, y1] + inputs: [x0, x1, x] + outputs: [y0, y1] + outputs: [x0, x1, y] + - ! frame: ! frames: - ! @@ -790,47 +645,45 @@ wcs: ! axes_order: [0, 1] axis_physical_types: ['custom:v2', 'custom:v3'] name: v2v3_spatial - unit: [!unit/unit-1.0.0 'arcsec', !unit/unit-1.0.0 'arcsec'] + unit: [!unit/unit-1.0.0 arcsec, !unit/unit-1.0.0 arcsec] - ! axes_names: [lambda] axes_order: [2] axis_physical_types: [em.wl] name: spec - unit: [!unit/unit-1.0.0 'um'] + unit: [!unit/unit-1.0.0 um] name: v2v3 transform: !transform/concatenate-1.2.0 - bounds: - angles_2: [null, null] - factor_0: [null, null] - factor_1: [null, null] - fixed: {angles_2: false, factor_0: false, factor_1: false} forward: - !transform/compose-1.2.0 forward: - !transform/concatenate-1.2.0 forward: - !transform/scale-1.2.0 - bounds: - factor: &id005 [null, null] factor: 0.0002777777777777778 - fixed: {factor: false} + inputs: [x] + outputs: [y] - !transform/scale-1.2.0 - bounds: - factor: *id005 factor: 0.0002777777777777778 - fixed: {factor: false} + inputs: [x] + outputs: [y] + inputs: [x0, x1] + outputs: [y0, y1] - !transform/rotate_sequence_3d-1.0.0 angles: [-0.12597594444444443, 0.10374517305555556, 0.0, -72.0545718, -5.630568] axes_order: zyxyz - bounds: - angles: [null, null] - fixed: {angles: false} + inputs: [lon, lat] name: v23tosky + outputs: [lon, lat] rotation_type: spherical + inputs: [x0, x1] + outputs: [lon, lat] - !transform/identity-1.2.0 - bounds: {} - fixed: {} - - ! + inputs: [x0] + outputs: [x0] + inputs: [x00, x10, x01] + outputs: [lon, lat, x0] + - ! frame: ! frames: - ! @@ -840,14 +693,15 @@ wcs: ! name: icrs reference_frame: ! frame_attributes: {} - unit: [!unit/unit-1.0.0 'deg', !unit/unit-1.0.0 'deg'] + unit: [!unit/unit-1.0.0 deg, !unit/unit-1.0.0 deg] - ! axes_names: [lambda] axes_order: [2] axis_physical_types: [em.wl] name: spec - unit: [!unit/unit-1.0.0 'um'] + unit: [!unit/unit-1.0.0 um] name: world + transform: null ... BLK0۴ke-.u33333)?BLK0Eտ5Nh{@D@BLK0۴ke-.u33333)?BLK0Eտ5Nh{@D@BLK0 Z惂<LCjX?_@ r+9:mM)3><ǹ?"KmL-?@n.>TF>6 \UnҪx->Fa>V9>`f>BLK0*W˲VN{C&IiujBCƞ>@ Ո8x:?Z4~\W?rfƾd;P>6X=#>`hg=eLD-BLK0ȓr)ǸW`ӷi?ϝ_2*0IAn d>([RѨ@H`f0'9YE (p]>o-L'? t>3MX/Egy>\` ^>BLK09UHw o+u F_|0W4?Wʽxн{D=_=mR?w8j{?h  X,,R[ތ>Iİ؃>%Mu,i"qF!Z9 @@ -930,26 +784,26 @@ A>8)@ #+@y)*+@J5_v2+@A<:+@L~B+@XJ+@dӅR+@_p tZ+@0|^b+@jIj+@ғG3r+@$z+@t+@E +@ûݑ+@Θ(ș+@u+@R0+@Z/+@+ 8r+@ \+@?G+@!1+@o-G+@@9^+@E;O+@P+@\V+@hڰ,@Ut^,@&,@ifp,@ȗFZ ,@#nE(,@j/0,@;u8,@ Ǻ@,@җ}G,@tO,@QW,@P. _,@! g,@ o,@Ɣnw,@%Y,@e1C,@6=] .,@I:,@T(,@`,@zl/خ,@Kx¶,@7,@h,@#ASDF BLOCK INDEX %YAML 1.1 --- -- 35600 -- 35670 -- 35740 -- 35810 -- 35880 -- 36134 -- 36388 -- 36642 -- 36896 -- 36982 -- 37068 -- 40226 -- 43384 -- 49030 -- 54676 -- 54762 -- 54848 -- 55102 -- 55356 -- 55426 -- 55496 -- 61142 +- 27843 +- 27913 +- 27983 +- 28053 +- 28123 +- 28377 +- 28631 +- 28885 +- 29139 +- 29225 +- 29311 +- 32469 +- 35627 +- 41273 +- 46919 +- 47005 +- 47091 +- 47345 +- 47599 +- 47669 +- 47739 +- 53385 ... diff --git a/gwcs/tests/data/miriwcs.asdf b/gwcs/tests/data/miriwcs.asdf index e70bc620..2cc3a875 100644 --- a/gwcs/tests/data/miriwcs.asdf +++ b/gwcs/tests/data/miriwcs.asdf @@ -1,237 +1,470 @@ #ASDF 1.0.0 -#ASDF_STANDARD 1.4.0 +#ASDF_STANDARD 1.5.0 %YAML 1.1 %TAG ! tag:stsci.edu:asdf/ --- !core/asdf-1.1.0 -asdf_library: !core/software-1.0.0 {author: Space Telescope Science Institute, homepage: 'http://github.com/spacetelescope/asdf', - name: asdf, version: 2.5.2} +asdf_library: !core/software-1.0.0 {author: The ASDF Developers, homepage: 'http://github.com/asdf-format/asdf', + name: asdf, version: 3.5.0} history: extensions: - !core/extension_metadata-1.0.0 - extension_class: astropy.io.misc.asdf.extension.AstropyAsdfExtension - software: {name: astropy, version: '4.0'} + extension_class: asdf.extension._manifest.ManifestExtension + extension_uri: asdf://asdf-format.org/core/extensions/core-1.5.0 + manifest_software: !core/software-1.0.0 {name: asdf_standard, version: 1.1.1} + software: !core/software-1.0.0 {name: asdf, version: 3.5.0} - !core/extension_metadata-1.0.0 - extension_class: gwcs.extension.GWCSExtension - software: {name: gwcs, version: 0.12.0} + extension_class: asdf.extension._manifest.ManifestExtension + extension_uri: asdf://astropy.org/astropy/extensions/units-1.0.0 + software: !core/software-1.0.0 {name: asdf-astropy, version: 0.6.1} - !core/extension_metadata-1.0.0 - extension_class: asdf.extension.BuiltinExtension - software: {name: asdf, version: 2.5.2} -wcs: ! + extension_class: asdf.extension._manifest.ManifestExtension + extension_uri: asdf://asdf-format.org/transform/extensions/transform-1.5.0 + manifest_software: !core/software-1.0.0 {name: asdf_transform_schemas, version: 0.5.0} + software: !core/software-1.0.0 {name: asdf-astropy, version: 0.6.1} + - !core/extension_metadata-1.0.0 + extension_class: asdf.extension._manifest.ManifestExtension + extension_uri: asdf://asdf-format.org/astronomy/gwcs/extensions/gwcs-1.2.0 + manifest_software: !core/software-1.0.0 {name: asdf_wcs_schemas, version: 0.4.0} + software: !core/software-1.0.0 {name: gwcs, version: 0.22.0a1.dev14+gc46e932} + - !core/extension_metadata-1.0.0 + extension_class: asdf.extension._manifest.ManifestExtension + extension_uri: asdf://asdf-format.org/astronomy/coordinates/extensions/coordinates-1.0.0 + manifest_software: !core/software-1.0.0 {name: asdf_coordinates_schemas, version: 0.3.0} + software: !core/software-1.0.0 {name: asdf-astropy, version: 0.6.1} +wcs: ! name: '' + pixel_shape: null steps: - - ! + - ! frame: ! axes_names: [x, y] axes_order: [0, 1] axis_physical_types: ['custom:x', 'custom:y'] name: detector - unit: [!unit/unit-1.0.0 'pixel', !unit/unit-1.0.0 'pixel'] - transform: !transform/compose-1.1.0 - bounding_box: - - [-0.5, 1023.5] - - [3.5, 1027.5] + unit: [!unit/unit-1.0.0 pixel, !unit/unit-1.0.0 pixel] + transform: !transform/compose-1.2.0 + bounding_box: !transform/property/bounding_box-1.0.0 + ignore: [] + intervals: + x0: [3.5, 1027.5] + x1: [-0.5, 1023.5] + order: F forward: - - !transform/concatenate-1.1.0 + - !transform/concatenate-1.2.0 forward: - - !transform/shift-1.2.0 {offset: 0.15} - - !transform/shift-1.2.0 {offset: -0.59} - - !transform/compose-1.1.0 + - !transform/shift-1.2.0 + inputs: [x] + offset: 0.15 + outputs: [y] + - !transform/shift-1.2.0 + inputs: [x] + offset: -0.59 + outputs: [y] + inputs: [x0, x1] + outputs: [y0, y1] + - !transform/compose-1.2.0 forward: - - !transform/compose-1.1.0 + - !transform/compose-1.2.0 forward: - - !transform/compose-1.1.0 + - !transform/compose-1.2.0 forward: - - !transform/compose-1.1.0 + - !transform/compose-1.2.0 forward: - - !transform/compose-1.1.0 + - !transform/compose-1.2.0 forward: - - !transform/compose-1.1.0 + - !transform/compose-1.2.0 forward: - - !transform/compose-1.1.0 + - !transform/compose-1.2.0 forward: - - !transform/concatenate-1.1.0 + - !transform/concatenate-1.2.0 forward: - - !transform/shift-1.2.0 {offset: -4.0} - - !transform/identity-1.1.0 {} - - !transform/concatenate-1.1.0 + - !transform/shift-1.2.0 + inputs: [x] + offset: -4.0 + outputs: [y] + - !transform/identity-1.2.0 + inputs: [x0] + outputs: [x0] + inputs: [x, x0] + outputs: [y, x0] + - !transform/concatenate-1.2.0 forward: - - !transform/polynomial-1.1.0 + - !transform/polynomial-1.2.0 coefficients: !core/ndarray-1.0.0 source: 0 datatype: float64 byteorder: little shape: [2] - inverse: !transform/polynomial-1.1.0 + domain: &id001 [-1, 1] + inputs: [x] + inverse: !transform/polynomial-1.2.0 coefficients: !core/ndarray-1.0.0 source: 1 datatype: float64 byteorder: little shape: [2] + domain: *id001 + inputs: [x] + outputs: [y] + window: *id001 name: M_column_correction - - !transform/polynomial-1.1.0 + outputs: [y] + window: *id001 + - !transform/polynomial-1.2.0 coefficients: !core/ndarray-1.0.0 source: 2 datatype: float64 byteorder: little shape: [2] - inverse: !transform/polynomial-1.1.0 + domain: *id001 + inputs: [x] + inverse: !transform/polynomial-1.2.0 coefficients: !core/ndarray-1.0.0 source: 3 datatype: float64 byteorder: little shape: [2] + domain: *id001 + inputs: [x] + outputs: [y] + window: *id001 name: M_row_correction - - !transform/remap_axes-1.1.0 - inverse: !transform/identity-1.1.0 {n_dims: 2} + outputs: [y] + window: *id001 + inputs: [x0, x1] + outputs: [y0, y1] + inputs: [x, x0] + outputs: [y0, y1] + - !transform/remap_axes-1.3.0 + inputs: [x0, x1] + inverse: !transform/identity-1.2.0 + inputs: [x0, x1] + n_dims: 2 + outputs: [x0, x1] mapping: [0, 1, 0, 1] - - !transform/concatenate-1.1.0 + outputs: [x0, x1, x2, x3] + inputs: [x, x0] + outputs: [x0, x1, x2, x3] + - !transform/concatenate-1.2.0 forward: - - !transform/polynomial-1.1.0 + - !transform/polynomial-1.2.0 coefficients: !core/ndarray-1.0.0 source: 4 datatype: float64 byteorder: little shape: [5, 5] - inverse: !transform/polynomial-1.1.0 + domain: + - *id001 + - *id001 + inputs: [x, y] + inverse: !transform/polynomial-1.2.0 coefficients: !core/ndarray-1.0.0 source: 5 datatype: float64 byteorder: little shape: [5, 5] + domain: + - *id001 + - *id001 + inputs: [x, y] + outputs: [z] + window: + - *id001 + - *id001 name: B_correction - - !transform/polynomial-1.1.0 + outputs: [z] + window: + - *id001 + - *id001 + - !transform/polynomial-1.2.0 coefficients: !core/ndarray-1.0.0 source: 6 datatype: float64 byteorder: little shape: [5, 5] - inverse: !transform/polynomial-1.1.0 + domain: + - *id001 + - *id001 + inputs: [x, y] + inverse: !transform/polynomial-1.2.0 coefficients: !core/ndarray-1.0.0 source: 7 datatype: float64 byteorder: little shape: [5, 5] + domain: + - *id001 + - *id001 + inputs: [x, y] + outputs: [z] + window: + - *id001 + - *id001 name: A_correction - - !transform/remap_axes-1.1.0 - inverse: !transform/remap_axes-1.1.0 + outputs: [z] + window: + - *id001 + - *id001 + inputs: [x0, y0, x1, y1] + outputs: [z0, z1] + inputs: [x, x0] + outputs: [z0, z1] + - !transform/remap_axes-1.3.0 + inputs: [x0, x1] + inverse: !transform/remap_axes-1.3.0 + inputs: [x0, x1] mapping: [0, 1, 0, 1] + outputs: [x0, x1, x2, x3] mapping: [0, 1, 0, 1] - - !transform/concatenate-1.1.0 + outputs: [x0, x1, x2, x3] + inputs: [x, x0] + outputs: [x0, x1, x2, x3] + - !transform/concatenate-1.2.0 forward: - - !transform/polynomial-1.1.0 + - !transform/polynomial-1.2.0 coefficients: !core/ndarray-1.0.0 source: 8 datatype: float64 byteorder: little shape: [2, 2] + domain: + - *id001 + - *id001 + inputs: [x, y] name: TI_row_correction - - !transform/polynomial-1.1.0 + outputs: [z] + window: + - *id001 + - *id001 + - !transform/polynomial-1.2.0 coefficients: !core/ndarray-1.0.0 source: 9 datatype: float64 byteorder: little shape: [2, 2] + domain: + - *id001 + - *id001 + inputs: [x, y] name: TI_column_correction - - !transform/identity-1.1.0 - inverse: !transform/remap_axes-1.1.0 + outputs: [z] + window: + - *id001 + - *id001 + inputs: [x0, y0, x1, y1] + outputs: [z0, z1] + inputs: [x, x0] + outputs: [z0, z1] + - !transform/identity-1.2.0 + inputs: [x0, x1] + inverse: !transform/remap_axes-1.3.0 + inputs: [x0, x1] mapping: [0, 1, 0, 1] + outputs: [x0, x1, x2, x3] n_dims: 2 - - !transform/remap_axes-1.1.0 + outputs: [x0, x1] + inputs: [x, x0] + outputs: [x0, x1] + - !transform/remap_axes-1.3.0 + inputs: [x0, x1] mapping: [1, 0] - inverse: !transform/compose-1.1.0 + outputs: [x0, x1] + inputs: [x, x0] + outputs: [x0, x1] + inputs: [x0, x1] + inverse: !transform/compose-1.2.0 forward: - - !transform/compose-1.1.0 + - !transform/compose-1.2.0 forward: - - !transform/remap_axes-1.1.0 + - !transform/remap_axes-1.3.0 + inputs: [x0, x1] mapping: [1, 0] - - !transform/compose-1.1.0 + outputs: [x0, x1] + - !transform/compose-1.2.0 forward: - - !transform/remap_axes-1.1.0 + - !transform/remap_axes-1.3.0 + inputs: [x0, x1] mapping: [0, 1, 0, 1] - - !transform/compose-1.1.0 + outputs: [x0, x1, x2, x3] + - !transform/compose-1.2.0 forward: - - !transform/concatenate-1.1.0 + - !transform/concatenate-1.2.0 forward: - - !transform/polynomial-1.1.0 + - !transform/polynomial-1.2.0 coefficients: !core/ndarray-1.0.0 source: 10 datatype: float64 byteorder: little shape: [2, 2] + domain: + - *id001 + - *id001 + inputs: [x, y] name: T_row_correction - - !transform/polynomial-1.1.0 + outputs: [z] + window: + - *id001 + - *id001 + - !transform/polynomial-1.2.0 coefficients: !core/ndarray-1.0.0 source: 11 datatype: float64 byteorder: little shape: [2, 2] + domain: + - *id001 + - *id001 + inputs: [x, y] name: T_column_correction - - !transform/compose-1.1.0 + outputs: [z] + window: + - *id001 + - *id001 + inputs: [x0, y0, x1, y1] + outputs: [z0, z1] + - !transform/compose-1.2.0 forward: - - !transform/remap_axes-1.1.0 + - !transform/remap_axes-1.3.0 + inputs: [x0, x1] mapping: [0, 1, 0, 1] - - !transform/compose-1.1.0 + outputs: [x0, x1, x2, x3] + - !transform/compose-1.2.0 forward: - - !transform/concatenate-1.1.0 + - !transform/concatenate-1.2.0 forward: - - !transform/polynomial-1.1.0 + - !transform/polynomial-1.2.0 coefficients: !core/ndarray-1.0.0 source: 12 datatype: float64 byteorder: little shape: [5, 5] - - !transform/polynomial-1.1.0 + domain: + - *id001 + - *id001 + inputs: [x, y] + outputs: [z] + window: + - *id001 + - *id001 + - !transform/polynomial-1.2.0 coefficients: !core/ndarray-1.0.0 source: 13 datatype: float64 byteorder: little shape: [5, 5] - - !transform/compose-1.1.0 + domain: + - *id001 + - *id001 + inputs: [x, y] + outputs: [z] + window: + - *id001 + - *id001 + inputs: [x0, y0, x1, y1] + outputs: [z0, z1] + - !transform/compose-1.2.0 forward: - - !transform/identity-1.1.0 {n_dims: 2} - - !transform/compose-1.1.0 + - !transform/identity-1.2.0 + inputs: [x0, x1] + n_dims: 2 + outputs: [x0, x1] + - !transform/compose-1.2.0 forward: - - !transform/concatenate-1.1.0 + - !transform/concatenate-1.2.0 forward: - - !transform/polynomial-1.1.0 + - !transform/polynomial-1.2.0 coefficients: !core/ndarray-1.0.0 source: 14 datatype: float64 byteorder: little shape: [2] - - !transform/polynomial-1.1.0 + domain: *id001 + inputs: [x] + outputs: [y] + window: *id001 + - !transform/polynomial-1.2.0 coefficients: !core/ndarray-1.0.0 source: 15 datatype: float64 byteorder: little shape: [2] - - !transform/concatenate-1.1.0 + domain: *id001 + inputs: [x] + outputs: [y] + window: *id001 + inputs: [x0, x1] + outputs: [y0, y1] + - !transform/concatenate-1.2.0 forward: - - !transform/shift-1.2.0 {offset: 4.0} - - !transform/identity-1.1.0 {} - - !transform/concatenate-1.1.0 + - !transform/shift-1.2.0 + inputs: [x] + offset: 4.0 + outputs: [y] + - !transform/identity-1.2.0 + inputs: [x0] + outputs: [x0] + inputs: [x, x0] + outputs: [y, x0] + inputs: [x0, x1] + outputs: [y, x0] + inputs: [x0, x1] + outputs: [y, x0] + inputs: [x0, y0, x1, y1] + outputs: [y, x0] + inputs: [x0, x1] + outputs: [y, x0] + inputs: [x0, y0, x1, y1] + outputs: [y, x0] + inputs: [x0, x1] + outputs: [y, x0] + inputs: [x0, x1] + outputs: [y, x0] + - !transform/concatenate-1.2.0 forward: - - !transform/shift-1.2.0 {offset: -0.15} - - !transform/shift-1.2.0 {offset: 0.59} - - ! + - !transform/shift-1.2.0 + inputs: [x] + offset: -0.15 + outputs: [y] + - !transform/shift-1.2.0 + inputs: [x] + offset: 0.59 + outputs: [y] + inputs: [x0, x1] + outputs: [y0, y1] + inputs: [x0, x1] + outputs: [y0, y1] + outputs: [x0, x1] + - ! frame: ! axes_names: [x, y] axes_order: [0, 1] axis_physical_types: ['custom:x', 'custom:y'] name: v2v3 - unit: [!unit/unit-1.0.0 'arcsec', !unit/unit-1.0.0 'arcsec'] - transform: !transform/compose-1.1.0 + unit: [!unit/unit-1.0.0 arcsec, !unit/unit-1.0.0 arcsec] + transform: !transform/compose-1.2.0 forward: - - !transform/concatenate-1.1.0 + - !transform/concatenate-1.2.0 forward: - - !transform/scale-1.2.0 {factor: 0.0002777777777777778} - - !transform/scale-1.2.0 {factor: 0.0002777777777777778} + - !transform/scale-1.2.0 + factor: 0.0002777777777777778 + inputs: [x] + outputs: [y] + - !transform/scale-1.2.0 + factor: 0.0002777777777777778 + inputs: [x] + outputs: [y] + inputs: [x0, x1] + outputs: [y0, y1] - !transform/rotate_sequence_3d-1.0.0 angles: [-0.12597594444444443, 0.10374517305555556, 0.0, -72.0545718, -5.630568] axes_order: zyxyz + inputs: [lon, lat] name: v23tosky + outputs: [lon, lat] rotation_type: spherical - - ! + inputs: [x0, x1] + outputs: [lon, lat] + - ! frame: ! axes_names: [lon, lat] axes_order: [0, 1] @@ -239,7 +472,8 @@ wcs: ! name: world reference_frame: ! frame_attributes: {} - unit: [!unit/unit-1.0.0 'deg', !unit/unit-1.0.0 'deg'] + unit: [!unit/unit-1.0.0 deg, !unit/unit-1.0.0 deg] + transform: null ... BLK0۴ke-.u33333)?BLK0Eտ5Nh{@D@BLK0۴ke-.u33333)?BLK0Eտ5Nh{@D@BLK0 Z惂<LCjX?_@ r+9:mM)3><ǹ?"KmL-?@n.>TF>6 \UnҪx->Fa>V9>`f>BLK0*W˲VN{C&IiujBCƞ>@ Ո8x:?Z4~\W?rfƾd;P>6X=#>`hg=eLD-BLK0ȓr)ǸW`ӷi?ϝ_2*0IAn d>([RѨ@H`f0'9YE (p]>o-L'? t>3MX/Egy>\` ^>BLK09UHw o+u F_|0W4?Wʽxн{D=_=mR?w8j{?h  X,,R[ތ>Iİ؃>%Mu,i"qF!Z9 @@ -249,20 +483,20 @@ X,,R ؐ=BLK0*W˲VN{C&IiujBCƞ>@ Ո8x:?Z4~\W?rfƾd;P>6X=#>`hg=eLD-BLK0Eտ5Nh{@D@BLK0Eտ5Nh{@D@#ASDF BLOCK INDEX %YAML 1.1 --- -- 9730 -- 9800 -- 9870 -- 9940 -- 10010 -- 10264 -- 10518 -- 10772 -- 11026 -- 11112 -- 11198 -- 11284 -- 11370 -- 11624 -- 11878 -- 11948 +- 17762 +- 17832 +- 17902 +- 17972 +- 18042 +- 18296 +- 18550 +- 18804 +- 19058 +- 19144 +- 19230 +- 19316 +- 19402 +- 19656 +- 19910 +- 19980 ... diff --git a/gwcs/tests/data/nircamwcs.asdf b/gwcs/tests/data/nircamwcs.asdf index f0e2c14b..c9d739db 100644 Binary files a/gwcs/tests/data/nircamwcs.asdf and b/gwcs/tests/data/nircamwcs.asdf differ diff --git a/gwcs/tests/test_api.py b/gwcs/tests/test_api.py index fd6f916c..cdb49c76 100644 --- a/gwcs/tests/test_api.py +++ b/gwcs/tests/test_api.py @@ -63,7 +63,6 @@ def wcs_ndim_types_units(request): @fixture_all_wcses def test_lowlevel_types(wcsobj): - pytest.importorskip("typeguard") try: # Skip this on older versions of astropy where it dosen't exist. from astropy.wcs.wcsapi.tests.utils import validate_low_level_wcs_types @@ -165,8 +164,10 @@ def test_array_index_to_world_values(gwcs_2d_spatial_shift, x, y): def test_world_axis_object_components_2d(gwcs_2d_spatial_shift): waoc = gwcs_2d_spatial_shift.world_axis_object_components - assert waoc == [('celestial', 0, 'spherical.lon'), - ('celestial', 1, 'spherical.lat')] + assert waoc[0][:2] == ('celestial', 0) + assert callable(waoc[0][2]) + assert waoc[1][:2] == ('celestial', 1) + assert callable(waoc[1][2]) def test_world_axis_object_components_2d_generic(gwcs_2d_quantity_shift): @@ -177,15 +178,19 @@ def test_world_axis_object_components_2d_generic(gwcs_2d_quantity_shift): def test_world_axis_object_components_1d(gwcs_1d_freq): waoc = gwcs_1d_freq.world_axis_object_components - assert waoc == [('spectral', 0, 'value')] + assert [c[:2] for c in waoc] == [('spectral', 0)] + assert callable(waoc[0][2]) def test_world_axis_object_components_4d(gwcs_4d_identity_units): waoc = gwcs_4d_identity_units.world_axis_object_components - assert waoc[0:3] == [('celestial', 0, 'spherical.lon'), - ('celestial', 1, 'spherical.lat'), - ('spectral', 0, 'value')] - assert waoc[3][0:2] == ('temporal', 0) + first_two = [c[:2] for c in waoc] + last_one = [c[2] for c in waoc] + assert first_two == [('celestial', 0), + ('celestial', 1), + ('spectral', 0), + ('temporal', 0)] + assert all([callable(l) for l in last_one]) def test_world_axis_object_classes_2d(gwcs_2d_spatial_shift): @@ -230,12 +235,12 @@ def test_world_axis_object_classes_4d(gwcs_4d_identity_units): def _compare_frame_output(wc1, wc2): if isinstance(wc1, coord.SkyCoord): assert isinstance(wc1.frame, type(wc2.frame)) - assert u.allclose(wc1.spherical.lon, wc2.spherical.lon) - assert u.allclose(wc1.spherical.lat, wc2.spherical.lat) - assert u.allclose(wc1.spherical.distance, wc2.spherical.distance) + assert u.allclose(wc1.spherical.lon, wc2.spherical.lon, equal_nan=True) + assert u.allclose(wc1.spherical.lat, wc2.spherical.lat, equal_nan=True) + assert u.allclose(wc1.spherical.distance, wc2.spherical.distance, equal_nan=True) elif isinstance(wc1, u.Quantity): - assert u.allclose(wc1, wc2) + assert u.allclose(wc1, wc2, equal_nan=True) elif isinstance(wc1, time.Time): assert u.allclose((wc1 - wc2).to(u.s), 0*u.s) @@ -252,12 +257,6 @@ def _compare_frame_output(wc1, wc2): @fixture_all_wcses def test_high_level_wrapper(wcsobj, request): - if request.node.callspec.params['wcsobj'] in ('gwcs_4d_identity_units', 'gwcs_stokes_lookup'): - pytest.importorskip("astropy", minversion="4.0dev0") - - # Remove the bounding box because the type test is a little broken with the - # bounding box. - del wcsobj._pipeline[0].transform.bounding_box hlvl = HighLevelWCSWrapper(wcsobj) @@ -280,8 +279,6 @@ def test_high_level_wrapper(wcsobj, request): def test_stokes_wrapper(gwcs_stokes_lookup): - pytest.importorskip("astropy", minversion="4.0dev0") - hlvl = HighLevelWCSWrapper(gwcs_stokes_lookup) pixel_input = [0, 1, 2, 3] @@ -463,10 +460,11 @@ def test_world_to_pixel(gwcs_2d_spatial_shift, sky_ra_dec): assert_allclose(wcsobj.world_to_pixel(sky), wcsobj.invert(ra, dec, with_units=False)) -def test_world_to_array_index(gwcs_2d_spatial_shift, sky_ra_dec): - wcsobj = gwcs_2d_spatial_shift +def test_world_to_array_index(gwcs_simple_imaging, sky_ra_dec): + wcsobj = gwcs_simple_imaging sky, ra, dec = sky_ra_dec - assert_allclose(wcsobj.world_to_array_index(sky), wcsobj.invert(ra, dec, with_units=False)[::-1]) + + assert_allclose(wcsobj.world_to_array_index(sky), wcsobj.invert(ra * u.deg, dec * u.deg, with_units=False)[::-1]) def test_world_to_pixel_values(gwcs_2d_spatial_shift, sky_ra_dec): @@ -476,12 +474,12 @@ def test_world_to_pixel_values(gwcs_2d_spatial_shift, sky_ra_dec): assert_allclose(wcsobj.world_to_pixel_values(sky), wcsobj.invert(ra, dec, with_units=False)) -def test_world_to_array_index_values(gwcs_2d_spatial_shift, sky_ra_dec): - wcsobj = gwcs_2d_spatial_shift +def test_world_to_array_index_values(gwcs_simple_imaging, sky_ra_dec): + wcsobj = gwcs_simple_imaging sky, ra, dec = sky_ra_dec assert_allclose(wcsobj.world_to_array_index_values(sky), - wcsobj.invert(ra, dec, with_units=False)[::-1]) + wcsobj.invert(ra * u.deg, dec * u.deg, with_units=False)[::-1]) def test_ndim_str_frames(gwcs_with_frames_strings): @@ -521,3 +519,19 @@ def test_coordinate_frame_api(): pixel2 = wcs.invert(world) assert u.allclose(pixel2, 0*u.pix) + + +def test_world_axis_object_components_units(gwcs_3d_identity_units): + from astropy.wcs.wcsapi.high_level_api import high_level_objects_to_values + + wcs = gwcs_3d_identity_units + world = wcs.pixel_to_world(1, 1, 1) + + values = high_level_objects_to_values(*world, low_level_wcs=wcs) + + expected_values = [world[0].spherical.lon.to_value(wcs.output_frame.unit[0]), + world[0].spherical.lon.to_value(wcs.output_frame.unit[1]), + world[1].to_value(wcs.output_frame.unit[2])] + + assert not any([isinstance(o, u.Quantity) for o in values]) + np.testing.assert_allclose(values, expected_values) diff --git a/gwcs/tests/test_api_slicing.py b/gwcs/tests/test_api_slicing.py index d4866242..626a0e2a 100644 --- a/gwcs/tests/test_api_slicing.py +++ b/gwcs/tests/test_api_slicing.py @@ -44,9 +44,13 @@ def test_ellipsis(gwcs_3d_galactic_spectral): assert_equal(wcs.axis_correlation_matrix, [[True, False, True], [False, True, False], [True, False, True]]) - assert wcs.world_axis_object_components == [('celestial', 1, 'spherical.lat'), - ('spectral', 0, 'value'), - ('celestial', 0, 'spherical.lon')] + first_two = [c[:2] for c in wcs.world_axis_object_components] + last_one = [c[2] for c in wcs.world_axis_object_components] + assert first_two == [('celestial', 1), + ('spectral', 0), + ('celestial', 0)] + + assert all([callable(l) for l in last_one]) assert wcs.world_axis_object_classes['celestial'][0] is SkyCoord assert wcs.world_axis_object_classes['celestial'][1] == () @@ -106,8 +110,12 @@ def test_spectral_slice(gwcs_3d_galactic_spectral): assert_equal(wcs.axis_correlation_matrix, [[True, True], [True, True]]) - assert wcs.world_axis_object_components == [('celestial', 1, 'spherical.lat'), - ('celestial', 0, 'spherical.lon')] + first_two = [c[:2] for c in wcs.world_axis_object_components] + last_one = [c[2] for c in wcs.world_axis_object_components] + assert first_two == [('celestial', 1), + ('celestial', 0)] + + assert all([callable(l) for l in last_one]) assert wcs.world_axis_object_classes['celestial'][0] is SkyCoord assert wcs.world_axis_object_classes['celestial'][1] == () @@ -166,9 +174,13 @@ def test_spectral_range(gwcs_3d_galactic_spectral): assert_equal(wcs.axis_correlation_matrix, [[True, False, True], [False, True, False], [True, False, True]]) - assert wcs.world_axis_object_components == [('celestial', 1, 'spherical.lat'), - ('spectral', 0, 'value'), - ('celestial', 0, 'spherical.lon')] + first_two = [c[:2] for c in wcs.world_axis_object_components] + last_one = [c[2] for c in wcs.world_axis_object_components] + assert first_two == [('celestial', 1), + ('spectral', 0), + ('celestial', 0)] + + assert all([callable(l) for l in last_one]) assert wcs.world_axis_object_classes['celestial'][0] is SkyCoord assert wcs.world_axis_object_classes['celestial'][1] == () @@ -230,9 +242,13 @@ def test_celestial_slice(gwcs_3d_galactic_spectral): assert_equal(wcs.axis_correlation_matrix, [[False, True], [True, False], [False, True]]) - assert wcs.world_axis_object_components == [('celestial', 1, 'spherical.lat'), - ('spectral', 0, 'value'), - ('celestial', 0, 'spherical.lon')] + first_two = [c[:2] for c in wcs.world_axis_object_components] + last_one = [c[2] for c in wcs.world_axis_object_components] + assert first_two == [('celestial', 1), + ('spectral', 0), + ('celestial', 0)] + + assert all([callable(l) for l in last_one]) assert wcs.world_axis_object_classes['celestial'][0] is SkyCoord assert wcs.world_axis_object_classes['celestial'][1] == () @@ -246,8 +262,8 @@ def test_celestial_slice(gwcs_3d_galactic_spectral): assert_allclose(wcs.pixel_to_world_values(39, 44), (10.24, 20, 25)) assert_allclose(wcs.array_index_to_world_values(44, 39), (10.24, 20, 25)) - assert_allclose(wcs.world_to_pixel_values(12.4, 20, 25), (39., 44.)) - assert_equal(wcs.world_to_array_index_values(12.4, 20, 25), (44, 39)) + assert_allclose(wcs.world_to_pixel_values(10.24, 20, 25), (39., 44.)) + assert_equal(wcs.world_to_array_index_values(10.24, 20, 25), (44, 39)) assert_equal(wcs.pixel_bounds, [(-2, 45), (5, 50)]) @@ -295,9 +311,13 @@ def test_celestial_range(gwcs_3d_galactic_spectral): assert_equal(wcs.axis_correlation_matrix, [[True, False, True], [False, True, False], [True, False, True]]) - assert wcs.world_axis_object_components == [('celestial', 1, 'spherical.lat'), - ('spectral', 0, 'value'), - ('celestial', 0, 'spherical.lon')] + first_two = [c[:2] for c in wcs.world_axis_object_components] + last_one = [c[2] for c in wcs.world_axis_object_components] + assert first_two == [('celestial', 1), + ('spectral', 0), + ('celestial', 0)] + + assert all([callable(l) for l in last_one]) assert wcs.world_axis_object_classes['celestial'][0] is SkyCoord assert wcs.world_axis_object_classes['celestial'][1] == () @@ -363,9 +383,13 @@ def test_no_array_shape(gwcs_3d_galactic_spectral): assert_equal(wcs.axis_correlation_matrix, [[True, False, True], [False, True, False], [True, False, True]]) - assert wcs.world_axis_object_components == [('celestial', 1, 'spherical.lat'), - ('spectral', 0, 'value'), - ('celestial', 0, 'spherical.lon')] + first_two = [c[:2] for c in wcs.world_axis_object_components] + last_one = [c[2] for c in wcs.world_axis_object_components] + assert first_two == [('celestial', 1), + ('spectral', 0), + ('celestial', 0)] + + assert all([callable(l) for l in last_one]) assert wcs.world_axis_object_classes['celestial'][0] is SkyCoord assert wcs.world_axis_object_classes['celestial'][1] == () @@ -431,9 +455,13 @@ def test_ellipsis_none_types(gwcs_3d_galactic_spectral): assert_equal(wcs.axis_correlation_matrix, [[True, False, True], [False, True, False], [True, False, True]]) - assert wcs.world_axis_object_components == [('celestial', 1, 'spherical.lat'), - ('spectral', 0, 'value'), - ('celestial', 0, 'spherical.lon')] + first_two = [c[:2] for c in wcs.world_axis_object_components] + last_one = [c[2] for c in wcs.world_axis_object_components] + assert first_two == [('celestial', 1), + ('spectral', 0), + ('celestial', 0)] + + assert all([callable(l) for l in last_one]) assert wcs.world_axis_object_classes['celestial'][0] is SkyCoord assert wcs.world_axis_object_classes['celestial'][1] == () diff --git a/gwcs/tests/test_bounding_box.py b/gwcs/tests/test_bounding_box.py new file mode 100644 index 00000000..1cc46404 --- /dev/null +++ b/gwcs/tests/test_bounding_box.py @@ -0,0 +1,87 @@ +import numpy as np +from numpy.testing import assert_array_equal, assert_allclose + +import pytest + + +x = [-1, 2, 4, 13] +y = [np.nan, np.nan, 4, np.nan] +y1 = [np.nan, np.nan, 4, np.nan] + + +@pytest.mark.parametrize((("input", "output")), [((2, 4), (2, 4)), + ((100, 200), (np.nan, np.nan)), + ((x, x),(y, y)) + ]) +def test_2d_spatial(gwcs_2d_spatial_shift, input, output): + w = gwcs_2d_spatial_shift + w.bounding_box = ((-.5, 21), (4, 12)) + + assert_array_equal(w.invert(*w(*input)), output) + assert_array_equal(w.world_to_pixel_values(*w.pixel_to_world_values(*input)), output) + assert_array_equal(w.world_to_pixel(w.pixel_to_world(*input)), output) + + +@pytest.mark.parametrize((("input", "output")), [((2, 4), (2, 4)), + ((100, 200), (np.nan, np.nan)), + ((x, x), (y, y)) + ]) +def test_2d_spatial_coordinate(gwcs_2d_quantity_shift, input, output): + w = gwcs_2d_quantity_shift + w.bounding_box = ((-.5, 21), (4, 12)) + + assert_array_equal(w.invert(*w(*input)), output) + assert_array_equal(w.world_to_pixel_values(*w.pixel_to_world_values(*input)), output) + assert_array_equal(w.world_to_pixel(*w.pixel_to_world(*input)), output) + + +@pytest.mark.parametrize((("input", "output")), [((2, 4), (2, 4)), + ((100, 200), (np.nan, np.nan)), + ((x, x), (y, y)) + ]) +def test_2d_spatial_coordinate_reordered(gwcs_2d_spatial_reordered, input, output): + w = gwcs_2d_spatial_reordered + w.bounding_box = ((-.5, 21), (4, 12)) + + assert_array_equal(w.invert(*w(*input)), output) + assert_array_equal(w.world_to_pixel_values(*w.pixel_to_world_values(*input)), output) + assert_array_equal(w.world_to_pixel(w.pixel_to_world(*input)), output) + + +@pytest.mark.parametrize((("input", "output")), [(2, 2), + ((10, 200), (10, np.nan)), + (x, (np.nan, 2, 4, 13)) + ]) +def test_1d_freq(gwcs_1d_freq, input, output): + w = gwcs_1d_freq + w.bounding_box = (-.5, 21) + print(f"input {input}, {output}") + assert_array_equal(w.invert(w(input)), output) + assert_array_equal(w.world_to_pixel_values(w.pixel_to_world_values(input)), output) + assert_array_equal(w.world_to_pixel(w.pixel_to_world(input)), output) + + +@pytest.mark.parametrize((("input", "output")), [((2, 4, 5), (2, 4, 5)), + ((100, 200, 5), (np.nan, np.nan, np.nan)), + ((x, x, x), (y1, y1, y1)) + ]) +def test_3d_spatial_wave(gwcs_3d_spatial_wave, input, output): + w = gwcs_3d_spatial_wave + w.bounding_box = ((-.5, 21), (4, 12), (3, 21)) + + assert_array_equal(w.invert(*w(*input)), output) + assert_array_equal(w.world_to_pixel_values(*w.pixel_to_world_values(*input)), output) + assert_array_equal(w.world_to_pixel(*w.pixel_to_world(*input)), output) + + +@pytest.mark.parametrize((("input", "output")), [((1, 2, 3, 4), (1., 2., 3., 4.)), + ((100, 3, 3, 3), (np.nan, 3, 3, 3)), + ((x, x, x, x), [[np.nan, 2., 4., 13.], + [np.nan, 2., 4., 13.], + [np.nan, 2., 4., 13.], + [np.nan, 2., 4., np.nan]]) + ]) +def test_gwcs_spec_cel_time_4d(gwcs_spec_cel_time_4d, input, output): + w = gwcs_spec_cel_time_4d + + assert_allclose(w.invert(*w(*input, with_bounding_box=False)), output, atol=1e-8) diff --git a/gwcs/tests/test_coordinate_systems.py b/gwcs/tests/test_coordinate_systems.py index 967657f8..ac2c7303 100644 --- a/gwcs/tests/test_coordinate_systems.py +++ b/gwcs/tests/test_coordinate_systems.py @@ -190,7 +190,6 @@ def test_temporal_relative(): assert a[1] == Time("2018-01-01T00:00:00") + 20 * u.s -@pytest.mark.skipif(astropy_version<"4", reason="Requires astropy 4.0 or higher") def test_temporal_absolute(): t = cf.TemporalFrame(reference_frame=Time([], format='isot')) assert t.coordinates("2018-01-01T00:00:00") == Time("2018-01-01T00:00:00") @@ -240,7 +239,6 @@ def test_coordinate_to_quantity_spectral(inp): (Time("2011-01-01T00:00:10"),), (10 * u.s,) ]) -@pytest.mark.skipif(astropy_version<"4", reason="Requires astropy 4.0 or higher.") def test_coordinate_to_quantity_temporal(inp): temp = cf.TemporalFrame(reference_frame=Time("2011-01-01T00:00:00"), unit=u.s) @@ -325,7 +323,6 @@ def test_coordinate_to_quantity_frame_2d(): assert_quantity_allclose(output, exp) -@pytest.mark.skipif(astropy_version<"4", reason="Requires astropy 4.0 or higher.") def test_coordinate_to_quantity_error(): frame = cf.Frame2D(unit=(u.one, u.arcsec)) with pytest.raises(ValueError): diff --git a/gwcs/tests/test_utils.py b/gwcs/tests/test_utils.py index e69ec536..748b1472 100644 --- a/gwcs/tests/test_utils.py +++ b/gwcs/tests/test_utils.py @@ -6,6 +6,8 @@ from astropy import units as u from astropy import coordinates as coord from astropy.modeling import models +from astropy import table + from astropy.tests.helper import assert_quantity_allclose import pytest from numpy.testing import assert_allclose @@ -104,6 +106,10 @@ def test_isnumerical(): assert gwutils.isnumerical(np.array(0, dtype='>f8')) assert gwutils.isnumerical(np.array(0, dtype='>i4')) + # check a table column + t = table.Table(data=[[1,2,3], [4,5,6]], names=['x', 'y']) + assert not gwutils.isnumerical(t['x']) + def test_get_values(): args = 2 * u.cm diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py index 77627176..23fc82c1 100644 --- a/gwcs/tests/test_wcs.py +++ b/gwcs/tests/test_wcs.py @@ -7,13 +7,15 @@ import numpy as np from numpy.testing import assert_allclose, assert_equal -from astropy.modeling import models +from astropy.modeling import models, bind_compound_bounding_box +from astropy.modeling.bounding_box import ModelBoundingBox from astropy import coordinates as coord from astropy.io import fits from astropy import units as u from astropy import wcs as astwcs from astropy.wcs import wcsapi from astropy.time import Time +from astropy.utils.introspection import minversion import asdf from .. import wcs @@ -23,6 +25,7 @@ from ..utils import CoordinateFrameError from .utils import _gwcs_from_hst_fits_wcs from . import data +from ..examples import gwcs_2d_bad_bounding_box_order data_path = os.path.split(os.path.abspath(data.__file__))[0] @@ -50,6 +53,14 @@ y = np.linspace(0, 1, ny) xv, yv = np.meshgrid(x, y) + +def asdf_open_memory_mapping_kwarg(memmap: bool) -> dict: + if minversion("asdf", "3.1.0"): + return {"memmap": memmap} + else : + return {"copy_arrays": not memmap} + + # Test initializing a WCS @@ -382,6 +393,61 @@ def test_grid_from_bounding_box_step(): with pytest.raises(ValueError): grid_from_bounding_box(bb, step=(1, 2, 1)) +def test_grid_from_model_bounding_box(): + bbox = ((-1, 1), (0, 1)) + # Truth grid + grid_truth = grid_from_bounding_box(bbox) + + # Create a bounding box + model = models.Const2D() & models.Const1D() + model.inputs = ("x", "y", "slit_name") + model.bounding_box = ModelBoundingBox( + { + "x": bbox[0], + "y": bbox[1], + }, + model=model, + ignored=["slit_name"], + order="F", + ) + grid = grid_from_bounding_box(model.bounding_box) + + assert np.all(grid == grid_truth) + + +def test_grid_from_compound_bounding_box(): + bbox = ((-1, 1), (0, 1)) + # Truth grid + grid_truth = grid_from_bounding_box(bbox) + + # Create a compound bounding box + model = models.Const2D() & models.Const1D() + model.inputs = ("x", "y", "slit_name") + bind_compound_bounding_box( + model, + { + (200,) : { + "x": bbox[0], + "y": bbox[1], + }, + (300,) :{ + "x": (-2, 2), + "y": (0, 2), + } + }, + [("slit_name",)], + order="F", + ) + grid = grid_from_bounding_box(model.bounding_box, selector=(200,)) + + assert np.all(grid == grid_truth) + + # Capture errors + with pytest.raises(ValueError, match=r"Cannot use selector with a non-CompoundBoundingBox"): + grid_from_bounding_box(model.bounding_box[(300,)], selector=(300,)) + with pytest.raises(ValueError, match=r"selector must be set when bounding_box is a CompoundBoundingBox"): + grid_from_bounding_box(model.bounding_box) + def test_wcs_from_points(): np.random.seed(0) @@ -414,6 +480,14 @@ def test_wcs_from_points(): assert_allclose(newra, ra, atol=10**-6) assert_allclose(newdec, dec, atol=10**-6) + newra, newdec = w.pixel_to_world_values(x, y) + assert_allclose(newra, ra, atol=10**-6) + assert_allclose(newdec, dec, atol=10**-6) + + newsky = w.pixel_to_world(x, y) + assert_allclose(newsky.data.lon.deg, ra, atol=10**-6) + assert_allclose(newsky.data.lat.deg, dec, atol=10**-6) + def test_grid_from_bounding_box_2(): bb = ((-0.5, 5.5), (-0.5, 4.5)) @@ -648,7 +722,7 @@ def test_to_fits_sip(): xflat = np.ravel(x[1:-1, 1:-1]) yflat = np.ravel(y[1:-1, 1:-1]) fn = os.path.join(data_path, 'miriwcs.asdf') - with asdf.open(fn, lazy_load=False, copy_arrays=True, ignore_missing_extensions=True) as af: + with asdf.open(fn, lazy_load=False, ignore_missing_extensions=True, **asdf_open_memory_mapping_kwarg(memmap=False)) as af: miriwcs = af.tree['wcs'] bounding_box = ((0, 1024), (0, 1024)) mirisip = miriwcs.to_fits_sip(bounding_box, max_inv_pix_error=0.1, verbose=True) @@ -1009,7 +1083,7 @@ def test_to_fits_tab_time_cube(gwcs_cube_with_separable_time): def test_to_fits_tab_miri_image(): # gWCS: fn = os.path.join(data_path, 'miriwcs.asdf') - with asdf.open(fn, copy_arrays=True, lazy_load=False, ignore_missing_extensions=True) as af: + with asdf.open(fn, lazy_load=False, ignore_missing_extensions=True, **asdf_open_memory_mapping_kwarg(memmap=False)) as af: w = af.tree['wcs'] # FITS WCS -TAB: @@ -1033,7 +1107,7 @@ def test_to_fits_tab_miri_image(): def test_to_fits_tab_miri_lrs(): fn = os.path.join(data_path, 'miri_lrs_wcs.asdf') - with asdf.open(fn, copy_arrays=True, lazy_load=False, ignore_missing_extensions=True) as af: + with asdf.open(fn, lazy_load=False, ignore_missing_extensions=True, **asdf_open_memory_mapping_kwarg(memmap=False)) as af: w = af.tree['wcs'] # FITS WCS -TAB: @@ -1089,8 +1163,8 @@ def test_in_image(): assert np.isscalar(w2.in_image(2, 6)) assert not np.isscalar(w2.in_image([2], [6])) - assert w2.in_image(4, 6) - assert not w2.in_image(5, 0) + assert (w2.in_image(4, 6)) + assert not (w2.in_image(5, 0)) assert np.array_equal( w2.in_image( [[9, 10, 11, 15], [8, 9, 67, 98], [2, 2, np.nan, 102]], @@ -1102,7 +1176,7 @@ def test_in_image(): def test_iter_inv(): fn = os.path.join(data_path, 'nircamwcs.asdf') - with asdf.open(fn, lazy_load=False, copy_arrays=False, ignore_missing_extensions=True) as af: + with asdf.open(fn, lazy_load=False, ignore_missing_extensions=True, **asdf_open_memory_mapping_kwarg(memmap=True)) as af: w = af.tree['wcs'] # remove analytic/user-supplied inverse: w.pipeline[0].transform.inverse = None @@ -1125,11 +1199,12 @@ def test_iter_inv(): *w(x, y), adaptive=True, detect_divergence=True, + tolerance=1e-4, maxiter=50, quiet=False ) assert np.allclose((x, y), (xp, yp)) - with asdf.open(fn, lazy_load=False, copy_arrays=False, ignore_missing_extensions=True) as af: + with asdf.open(fn, lazy_load=False, ignore_missing_extensions=True, **asdf_open_memory_mapping_kwarg(memmap=True)) as af: w = af.tree['wcs'] # test single point @@ -1144,6 +1219,7 @@ def test_iter_inv(): xp, yp = w.numerical_inverse( *w(x, y), adaptive=True, + tolerance=1e-5, maxiter=50, detect_divergence=False, quiet=False ) @@ -1178,6 +1254,7 @@ def test_iter_inv(): xp, yp = w.numerical_inverse( *w(x, y, with_bounding_box=False), adaptive=False, + tolerance=1e-5, maxiter=50, detect_divergence=True, quiet=False, with_bounding_box=False @@ -1191,6 +1268,7 @@ def test_iter_inv(): xp, yp = w.numerical_inverse( *w(x, y, with_bounding_box=False), adaptive=False, + tolerance=1e-5, maxiter=50, detect_divergence=True, quiet=False, with_bounding_box=False @@ -1264,7 +1342,7 @@ def test_sip_roundtrip(): assert hdr[k] == hdr_back[k] for k in ['cd1_1', 'cd1_2', 'cd2_1', 'cd2_2']: - assert np.allclose(hdr[k], hdr_back[k], atol=0, rtol=1e-9) + assert np.allclose(hdr[k], hdr_back[k], atol=1e-14, rtol=1e-9) for t in ('a', 'b'): order = hdr[f'{t}_order'] @@ -1275,7 +1353,7 @@ def test_sip_roundtrip(): assert np.allclose( hdr[k], hdr_back[k], - atol=0.0, + atol=1e-15, rtol=1.0e-8 * 10**(i + j) ) @@ -1321,3 +1399,78 @@ def test_spatial_spectral_stokes(): def test_wcs_str(): w = wcs.WCS(output_frame="icrs") assert 'icrs' in str(w) + + +def test_bounding_box_is_returned_F(): + bbox_tuple = ((1, 2), (3, 4)) + + detector_2d_frame = cf.Frame2D(name='detector', axes_order=(0, 1)) + model_2d_shift = models.Shift(1) & models.Shift(2) + + model_2d_shift_bbox = model_2d_shift.copy() + model_2d_shift_bbox.bounding_box = bbox_tuple + + frame = cf.CoordinateFrame(name="quantity", axes_order=(0, 1), naxes=2, axes_type=("SPATIAL", "SPATIAL"), unit=(u.km, u.km)) + + # Demonstrate that model_2d_shift does not have a bounding box + with pytest.raises(NotImplementedError): + model_2d_shift.bounding_box + + # Demonstrate that model_2d_shift_bbox does have a bounding box + assert model_2d_shift_bbox.bounding_box == bbox_tuple + + # Demonstrate the model_2d_shift_bbox has order "C" + assert model_2d_shift_bbox.bounding_box.order == "C" + + # Create a WCS and then set a bounding box on it + pipeline_bbox_after = [(detector_2d_frame, model_2d_shift), (frame, None)] + gwcs_object_after = wcs.WCS(pipeline_bbox_after) + gwcs_object_after.bounding_box = bbox_tuple + + assert gwcs_object_after.bounding_box == bbox_tuple + assert gwcs_object_after.bounding_box.order == "F" + + # Create a WCS on transform with a bounding box + pipeline_bbox_before = [(detector_2d_frame, model_2d_shift_bbox), (frame, None)] + gwcs_object_before = wcs.WCS(pipeline_bbox_before) + + # Check that first access in this case will raise a warning + with pytest.warns(wcs.GwcsBoundingBoxWarning): + gwcs_object_before.bounding_box + + # Check order is returned as F + assert gwcs_object_before.bounding_box.order == "F" + + # The bounding box tuple will now be ordered differently than the original + # Tuple due to the order change + assert gwcs_object_before.bounding_box != bbox_tuple + assert gwcs_object_before.bounding_box.bounding_box(order="C") == bbox_tuple + + # Show the the bounding box is different between the two WCS objects + assert gwcs_object_after.bounding_box != gwcs_object_before.bounding_box + + +def test_no_bounding_box_if_read_from_file(tmp_path): + bad_wcs = gwcs_2d_bad_bounding_box_order() + + # Check the waring is issued for the bounding box of this WCS object + with pytest.warns(wcs.GwcsBoundingBoxWarning): + bad_wcs.bounding_box + + # Check that the warning is not issued again the second time + with warnings.catch_warnings(): + warnings.simplefilter("error") + bad_wcs.bounding_box + + # Write a bad wcs bounding box to an asdf file + asdf_file = tmp_path / "bad_wcs.asdf" + af = asdf.AsdfFile({"wcs": gwcs_2d_bad_bounding_box_order()}) # re-create the bad wcs object + af.write_to(asdf_file) + + with asdf.open(asdf_file) as af: + wcs_from_file = af["wcs"] + + # Check that no warning is issued for the bounding box of this WCS object + with warnings.catch_warnings(): + warnings.simplefilter("error") + wcs_from_file.bounding_box diff --git a/gwcs/utils.py b/gwcs/utils.py index 104558cf..dcae0558 100644 --- a/gwcs/utils.py +++ b/gwcs/utils.py @@ -12,6 +12,7 @@ from astropy import coordinates as coords from astropy import units as u from astropy.time import Time, TimeDelta +from astropy import table from astropy.wcs import Celprm @@ -470,14 +471,12 @@ def isnumerical(val): Determine if a value is numerical (number or np.array of numbers). """ isnum = True - if isinstance(val, coords.SkyCoord): - isnum = False - elif isinstance(val, u.Quantity): - isnum = False - elif isinstance(val, (Time, TimeDelta)): + astropy_types=(coords.SkyCoord, u.Quantity, Time, TimeDelta, table.Column, table.Row) + if isinstance(val, astropy_types): isnum = False elif (isinstance(val, np.ndarray) and not np.issubdtype(val.dtype, np.floating) - and not np.issubdtype(val.dtype, np.integer)): + and not np.issubdtype(val.dtype, np.integer) + ): isnum = False return isnum diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 0432e541..e0cbf191 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -6,7 +6,6 @@ import astropy.io.fits as fits import numpy as np import numpy.linalg as npla -from astropy import units as u from astropy.modeling import fix_inputs, projections from astropy.modeling.bounding_box import CompoundBoundingBox from astropy.modeling.bounding_box import ModelBoundingBox as Bbox @@ -14,7 +13,11 @@ from astropy.modeling.models import (Const1D, Identity, Mapping, Polynomial2D, RotateCelestial2Native, Shift, Sky2Pix_TAN) +from astropy.modeling.parameters import _tofloat from astropy.wcs.utils import celestial_frame_to_wcs, proj_plane_pixel_scales +from astropy.wcs.wcsapi.high_level_api import high_level_objects_to_values, values_to_high_level_objects + +from astropy import units as u from scipy import linalg, optimize from . import coordinate_frames as cf @@ -73,6 +76,12 @@ def __init__(self, *args, best_solution=None, accuracy=None, niter=None, self.slow_conv = slow_conv +class GwcsBoundingBoxWarning(UserWarning): + """ + A warning class to report issues with bounding boxes in GWCS. + """ + + class _WorldAxisInfo(): def __init__(self, axis, frame, world_axis_order, cunit, ctype, input_axes): """ @@ -253,15 +262,20 @@ def set_transform(self, from_frame, to_frame, transform): def forward_transform(self): """ Return the total forward transform - from input to output coordinate frame. - """ - - if self._pipeline: - #return functools.reduce(lambda x, y: x | y, [step[1] for step in self._pipeline[: -1]]) - return functools.reduce(lambda x, y: x | y, [step.transform for step in self._pipeline[:-1]]) - else: + if not self._pipeline: return None + transform = functools.reduce(lambda x, y: x | y, [step.transform for step in self._pipeline[:-1]]) + + if self.bounding_box is not None: + # Currently compound models do not attempt to combine individual model + # bounding boxes. Get the forward transform and assign the bounding_box to it + # before evaluating it. The order Model.bounding_box is reversed. + transform.bounding_box = self.bounding_box + + return transform + @property def backward_transform(self): """ @@ -349,12 +363,6 @@ def __call__(self, *args, **kwargs): if 'fill_value' not in kwargs: kwargs['fill_value'] = np.nan - if self.bounding_box is not None: - # Currently compound models do not attempt to combine individual model - # bounding boxes. Get the forward transform and assign the bounding_box to it - # before evaluating it. The order Model.bounding_box is reversed. - transform.bounding_box = self.bounding_box - result = transform(*args, **kwargs) if with_units: @@ -392,34 +400,12 @@ def in_image(self, *args, **kwargs): and `False` if input is outside the footprint. """ - kwargs['with_bounding_box'] = True - kwargs['fill_value'] = np.nan - coords = self.invert(*args, **kwargs) result = np.isfinite(coords) if self.input_frame.naxes > 1: result = np.all(result, axis=0) - if self.bounding_box is None or not np.any(result): - return result - - if self.input_frame.naxes == 1: - x1, x2 = self.bounding_box.bounding_box() - - if len(np.shape(args[0])) > 0: - result[result] = (coords[result] >= x1) & (coords[result] <= x2) - elif result: - result = (coords >= x1) and (coords <= x2) - - else: - if len(np.shape(args[0])) > 0: - for c, (x1, x2) in zip(coords, self.bounding_box): - result[result] = (c[result] >= x1) & (c[result] <= x2) - - elif result: - result = all([(c >= x1) and (c <= x2) for c, (x1, x2) in zip(coords, self.bounding_box)]) - return result def invert(self, *args, **kwargs): @@ -466,29 +452,35 @@ def invert(self, *args, **kwargs): """ with_units = kwargs.pop('with_units', False) + try: + btrans = self.backward_transform + except NotImplementedError: + btrans = None if not utils.isnumerical(args[0]): + # convert astropy objects to numbers and arrays args = self.output_frame.coordinate_to_quantity(*args) if self.output_frame.naxes == 1: args = [args] - try: - if not self.backward_transform.uses_quantity: - args = utils.get_values(self.output_frame.unit, *args) - except (NotImplementedError, KeyError): - args = utils.get_values(self.output_frame.unit, *args) - if 'with_bounding_box' not in kwargs: - kwargs['with_bounding_box'] = True + # if the transform does not use units, getthe numerical values + if btrans is not None and not btrans.uses_quantity: + args = utils.get_values(self.output_frame.unit, *args) - if 'fill_value' not in kwargs: - kwargs['fill_value'] = np.nan + with_bounding_box = kwargs.pop('with_bounding_box', True) + fill_value = kwargs.pop('fill_value', np.nan) + akwargs = {k: v for k, v in kwargs.items() if k not in _ITER_INV_KWARGS} + if with_bounding_box and self.bounding_box is not None: + args = self.outside_footprint(args) - try: - # remove iterative inverse-specific keyword arguments: - akwargs = {k: v for k, v in kwargs.items() if k not in _ITER_INV_KWARGS} - result = self.backward_transform(*args, **akwargs) - except (NotImplementedError, KeyError): + if btrans is not None: + result = btrans(*args, **akwargs) + else: result = self.numerical_inverse(*args, **kwargs, with_units=with_units) + # deal with values outside the bounding box + if with_bounding_box and self.bounding_box is not None: + result = self.out_of_bounds(result, fill_value=fill_value) + if with_units and self.input_frame: if self.input_frame.naxes == 1: return self.input_frame.coordinates(result) @@ -497,7 +489,69 @@ def invert(self, *args, **kwargs): else: return result - def numerical_inverse(self, *args, tolerance=1e-5, maxiter=50, adaptive=True, + def outside_footprint(self, world_arrays): + world_arrays = list(world_arrays) + + axes_types = set(self.output_frame.axes_type) + axes_phys_types = self.world_axis_physical_types + footprint = self.footprint() + not_numerical = False + if not utils.isnumerical(world_arrays[0]): + not_numerical = True + world_arrays = high_level_objects_to_values(*world_arrays, low_level_wcs=self) + for axtyp in axes_types: + ind = np.asarray((np.asarray(self.output_frame.axes_type) == axtyp)) + + for idim, (coord, phys) in enumerate(zip(world_arrays, axes_phys_types)): + coord = _tofloat(coord) + if np.asarray(ind).sum() > 1: + axis_range = footprint[:, idim] + else: + axis_range = footprint + range = [axis_range.min(), axis_range.max()] + + if (axtyp == 'SPATIAL' and str(phys).endswith((".ra", ".lon")) + and range[1] - range[0] > 180): + # most likely this coordinate is wrapped at 360 + d = np.mean(range) + range = [ + axis_range[axis_range < d].max(), + axis_range[axis_range > d].min() + ] + outside = (coord >= range[0]) & (coord < range[1]) + else: + outside = (coord < range[0]) | (coord > range[1]) + if np.any(outside): + if np.isscalar(coord): + coord = np.nan + else: + coord[outside] = np.nan + world_arrays[idim] = coord + if not_numerical: + world_arrays = values_to_high_level_objects(*world_arrays, low_level_wcs=self) + return world_arrays + + + def out_of_bounds(self, pixel_arrays, fill_value=np.nan): + if np.isscalar(pixel_arrays) or self.input_frame.naxes == 1: + pixel_arrays = [pixel_arrays] + + pixel_arrays = list(pixel_arrays) + bbox = self.bounding_box + for idim, pix in enumerate(pixel_arrays): + outside = (pix < bbox[idim][0]) | (pix > bbox[idim][1]) + if np.any(outside): + if np.isscalar(pix): + pixel_arrays[idim] = np.nan + else: + pix = pixel_arrays[idim].astype(float, copy=True) + pix[outside] = np.nan + pixel_arrays[idim] = pix + if self.input_frame.naxes == 1: + pixel_arrays = pixel_arrays[0] + return pixel_arrays + + def numerical_inverse(self, *args, tolerance=1e-5, maxiter=30, adaptive=True, detect_divergence=True, quiet=True, with_bounding_box=True, fill_value=np.nan, with_units=False, **kwargs): """ @@ -679,7 +733,7 @@ def numerical_inverse(self, *args, tolerance=1e-5, maxiter=50, adaptive=True, >>> import numpy as np >>> filename = get_pkg_data_filename('data/nircamwcs.asdf', package='gwcs.tests') - >>> with asdf.open(filename, copy_arrays=True, lazy_load=False, ignore_missing_extensions=True) as af: + >>> with asdf.open(filename, lazy_load=False, ignore_missing_extensions=True) as af: ... w = af.tree['wcs'] >>> ra, dec = w([1,2,3], [1,1,1]) @@ -1282,11 +1336,32 @@ def bounding_box(self): frames = self.available_frames transform_0 = self.get_transform(frames[0], frames[1]) + + if transform_0 is None: + return None + try: bb = transform_0.bounding_box except NotImplementedError: return None + if ( + # Check that the bounding_box was set on the instance (not a default) + transform_0._user_bounding_box is not None + # Check the order of that bounding_box is C + and bb.order == "C" + # Check that the bounding_box is not a single value + and (isinstance(bb, CompoundBoundingBox) or len(bb) > 1) + ): + warnings.warn( + "The bounding_box was set in C order on the transform prior to being used in the gwcs!\n" + "Check that you indended that ordering for the bounding_box, and consider setting it in F order.\n" + "The bounding_box will remain meaning the same but will be converted to F order for consistency in the GWCS.", + GwcsBoundingBoxWarning + ) + self.bounding_box = bb.bounding_box(order="F") + bb = self.bounding_box + return bb @bounding_box.setter @@ -1386,17 +1461,22 @@ def _order_clockwise(v): if bounding_box is None: if self.bounding_box is None: raise TypeError("Need a valid bounding_box to compute the footprint.") - bb = self.bounding_box + bb = self.bounding_box.bounding_box(order='F') else: bb = bounding_box all_spatial = all([t.lower() == "spatial" for t in self.output_frame.axes_type]) - - if all_spatial: + if self.output_frame.naxes == 1: + if isinstance(bb[0], u.Quantity): + bb = np.asarray([b.value for b in bb]) * bb[0].unit + vertices = (bb,) + elif all_spatial: vertices = _order_clockwise(bb) else: vertices = np.array(list(itertools.product(*bb))).T + # workaround an issue with bbox with quantity, interval needs to be a cquantity, not a list of quantities + # strip units if center: vertices = utils._toindex(vertices) @@ -1410,14 +1490,16 @@ def _order_clockwise(v): axtyp_ind = np.array([t.lower() for t in self.output_frame.axes_type]) == axis_type if not axtyp_ind.any(): raise ValueError('This WCS does not have axis of type "{}".'.format(axis_type)) - result = np.asarray([(r.min(), r.max()) for r in result[axtyp_ind]]) + if len(axtyp_ind) > 1: + result = np.asarray([(r.min(), r.max()) for r in result[axtyp_ind]]) if axis_type == "spatial": result = _order_clockwise(result) else: result.sort() result = np.squeeze(result) - + if self.output_frame.naxes == 1: + return np.array([result]).T return result.T def fix_inputs(self, fixed): @@ -2031,16 +2113,12 @@ def find_frame(axis_number): # index of the axis in this frame's fidx = frame.axes_order.index(axno) - if hasattr(frame.unit[fidx], 'get_format_name'): - cunit = frame.unit[fidx].get_format_name(u.format.Fits).upper() - else: - cunit = '' axis_info = _WorldAxisInfo( axis=axno, frame=frame, world_axis_order=self.output_frame.axes_order.index(axno), - cunit=cunit, + cunit=frame.unit[fidx].to_string('fits', fraction=True).upper(), ctype=cf.get_ctype_from_ucd(self.world_axis_physical_types[axno]), input_axes=mapping[axno] ) diff --git a/gwcs/wcstools.py b/gwcs/wcstools.py index 4de18578..179f7773 100644 --- a/gwcs/wcstools.py +++ b/gwcs/wcstools.py @@ -5,6 +5,7 @@ from astropy.modeling.core import Model from astropy.modeling import projections from astropy.modeling import models, fitting +from astropy.modeling.bounding_box import CompoundBoundingBox, ModelBoundingBox from astropy import coordinates as coord from astropy import units as u @@ -139,7 +140,7 @@ def _frame2D_transform(fiducial, **kwargs): } -def grid_from_bounding_box(bounding_box, step=1, center=True): +def grid_from_bounding_box(bounding_box, step=1, center=True, selector=None): """ Create a grid of input points from the WCS bounding_box. @@ -151,11 +152,14 @@ def grid_from_bounding_box(bounding_box, step=1, center=True): Parameters ---------- - bounding_box : tuple + bounding_box : tuple | ~astropy.modeling.bounding_box.ModelBoundingBox | ~astropy.modeling.bounding_box.CompoundBoundingBox The bounding_box of a WCS object, `~gwcs.wcs.WCS.bounding_box`. step : scalar or tuple Step size for grid in each dimension. Scalar applies to all dimensions. center : bool + selector : tuple | None + If selector is set then it must be a selector tuple and bounding_box must + be a CompoundBoundingBox. The bounding_box is in order of X, Y [, Z] and the output will be in the same order. @@ -187,6 +191,26 @@ def grid_from_bounding_box(bounding_box, step=1, center=True): """ def _bbox_to_pixel(bbox): return (np.floor(bbox[0] + 0.5), np.ceil(bbox[1] - 0.5)) + + if selector is not None and not isinstance(bounding_box, CompoundBoundingBox): + raise ValueError("Cannot use selector with a non-CompoundBoundingBox") + + if isinstance(bounding_box, CompoundBoundingBox): + if selector is None: + raise ValueError("selector must be set when bounding_box is a CompoundBoundingBox") + + bounding_box = bounding_box[selector] + + if isinstance(bounding_box, ModelBoundingBox): + input_names = bounding_box.model.inputs + + # Get tuple of tuples of the bounding box values + bounding_box = tuple( + tuple(bounding_box[name]) + for name in input_names + if name not in bounding_box.ignored_inputs + ) + # 1D case if np.isscalar(bounding_box[0]): nd = 1 @@ -247,8 +271,9 @@ def wcs_from_points(xy, world_coords, proj_point='center', be passed in. projection : `~astropy.modeling.projections.Projection` A projection type. One of the projections in - `~astropy.modeling.projections.projcodes`. Defaults to TAN projection - (`astropy.modeling.projections.Sky2Pix_TAN`). + `~astropy.modeling.projections.projcodes`. + The direction is from sky to detector. + Defaults to TAN projection (`astropy.modeling.projections.Sky2Pix_TAN`). poly_degree : int Degree of polynomial model to be fit to data. Defaults to 4. polynomial_type : str @@ -302,7 +327,7 @@ def wcs_from_points(xy, world_coords, proj_point='center', "Only one of {} is supported.".format(polynomial_type, supported_poly_types.keys())) - skyrot = models.RotateCelestial2Native(crval[0], crval[1], 180*u.deg) + skyrot = models.RotateCelestial2Native(crval[0].deg, crval[1].deg, 180) trans = (skyrot | projection) projection_x, projection_y = trans(lon, lat) poly = supported_poly_types[polynomial_type](poly_degree) diff --git a/pyproject.toml b/pyproject.toml index 97aa010b..129d10fa 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,17 +1,17 @@ [project] name = "gwcs" description = "Generalized World Coordinate System" -requires-python = ">=3.9" +requires-python = ">=3.10" authors = [ { name = "gwcs developers", email = "help@stsci.edu" }, ] dependencies = [ - "asdf >= 2.8.1", - "astropy >= 5.3", - "numpy", - "scipy", + "asdf >= 3.3.0", + "astropy >= 6.0", + "numpy>=1.24", + "scipy>=1.14.1", "asdf_wcs_schemas >= 0.4.0", - "asdf-astropy >= 0.2.0", + "asdf-astropy >= 0.5.0", ] dynamic = [ "version", @@ -46,20 +46,21 @@ docs = [ ] test = [ "ci-watson>=0.3.0", - "pytest>=4.6.0", - "pytest-astropy", + "pytest>=7.0.0", + "pytest-astropy>=0.11.0", ] [build-system] requires = [ "setuptools>=61.2", "setuptools_scm[toml]>=3.4", - "wheel", ] build-backend = "setuptools.build_meta" [tool.setuptools.package-data] -'gwcs.tests.data' = ["*"] +"gwcs.tests.data" = [ + "*", +] [tool.setuptools.packages.find] namespaces = false @@ -74,17 +75,21 @@ upload-dir = "docs/_build/html" show-response = 1 [tool.pytest.ini_options] -minversion = "4.6" +minversion = 4.6 +doctest_plus = true +doctest_rst = true +text_file_format = "rst" +addopts = [ + "--color=yes", + "--doctest-rst", +] norecursedirs = [ "build", "docs/_build", ".tox", ] -doctest_plus = "enabled" -addopts = "--doctest-rst" filterwarnings = [ "ignore:Models in math_functions:astropy.utils.exceptions.AstropyUserWarning", - # importing astropy code causes this warning to be issued, so ignore it "ignore:numpy.ndarray size changed:RuntimeWarning", ] diff --git a/requirements-dev.txt b/requirements-dev.txt index 2ad2d68e..8ca4d21a 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -1,15 +1,13 @@ git+https://github.com/asdf-format/asdf git+https://github.com/asdf-format/asdf-standard -# Use weekly astropy dev build ---extra-index-url https://pypi.anaconda.org/astropy/simple astropy --pre - git+https://github.com/astropy/asdf-astropy git+https://github.com/asdf-format/asdf-transform-schemas git+https://github.com/asdf-format/asdf-coordinates-schemas git+https://github.com/asdf-format/asdf-wcs-schemas -# Use Bi-weekly numpy/scipy dev builds ---extra-index-url https://pypi.anaconda.org/scientific-python-nightly-wheels/simple +astropy>=0.0.dev0 +pyerfa>=0.0.dev0 + numpy>=0.0.dev0 scipy>=0.0.dev0 diff --git a/tox.ini b/tox.ini index e24fdcc5..17021f3a 100644 --- a/tox.ini +++ b/tox.ini @@ -1,9 +1,10 @@ [tox] envlist = - check-{style,security,build} - test{,-dev}{,-pyargs,-cov} - test-numpy{125,122} - build-{docs,dist} + check-{style,security} + test{,-dev}{,-oldestdeps,-cov}{-jwst,-romancal,-romanisim,-specutils,-dkist,-ndcube} + docs +requres = + tox-uv # tox environments are constructed with so-called 'factors' (or terms) # separated by hyphens, e.g. test-devdeps-cov. Lines below starting with factor: @@ -13,6 +14,13 @@ envlist = # # tox -l -v # +[main] +jwst_repo = https://github.com/spacetelescope/jwst.git +romancal_repo = https://github.com/spacetelescope/romancal.git +romanisim_repo = https://github.com/spacetelescope/romanisim.git +specutils_repo = https://github.com/astropy/specutils.git +dkist_repo = https://github.com/DKISTDC/dkist.git +ndcube_repo = https://github.com/sunpy/ndcube.git [testenv:check-style] description = check code style, e.g. with flake8 @@ -30,73 +38,78 @@ deps = commands = bandit -r -ll -c .bandit.yaml gwcs -[testenv:check-build] -description = check build sdist/wheel and a strict twine check for metadata -skip_install = true -deps = - twine>=3.3 - build +[testenv:docs] +description = invoke sphinx-build to build the HTML docs +extras = docs commands = - python -m build . - twine check --strict dist/* + sphinx-build -W docs docs/_build [testenv] description = run tests jwst: of JWST pipeline romancal: of Romancal pipeline + romanisim: of Romanisim image simulation + specutils: of Specutils + dkist: of DKIST + ndcube: of NDCube dev: with the latest developer version of key dependencies - pyargs: with --pyargs on installed package - warnings: treating warnings as errors + oldestdeps: with the oldest supported version of key dependencies cov: with coverage xdist: using parallel processing -passenv = +allowlist_externals = + jwst,romancal,romanisim,specutils,dkist,ndcube: git +pass_env = HOME GITHUB_* TOXENV CI CODECOV_* DISPLAY + jwst,romancal: CRDS_* + romanisim,romancal: WEBBPSF_PATH + romanisim: GALSIM_CAT_PATH + romanisim: FFTW_DIR + romanisim: LIBRARY_PATH set_env = - dev: PIP_EXTRA_INDEX_URL = https://pypi.anaconda.org/astropy/simple https://pypi.anaconda.org/liberfa/simple https://pypi.anaconda.org/scientific-python-nightly-wheels/simple - -args_are_paths = false -change_dir = pyargs: {env:HOME} + dev: UV_INDEX = https://pypi.anaconda.org/astropy/simple https://pypi.anaconda.org/liberfa/simple https://pypi.anaconda.org/scientific-python-nightly-wheels/simple + dev: UV_INDEX_STRATEGY = unsafe-any-match +change_dir = + jwst,romancal,romanisim,specutils,dkist,ndcube: {env_tmp_dir} extras = test alldeps: all +uv_resolution = + oldestdeps: lowest-direct deps = xdist: pytest-xdist cov: pytest-cov - jwst: jwst[test] @ git+https://github.com/spacetelescope/jwst.git - romancal: romancal[test] @ git+https://github.com/spacetelescope/romancal.git - numpy123: numpy==1.23.* - numpy125: numpy==1.25.* -pass_env = - jwst,romancal: CRDS_* + jwst: jwst[test] @ git+{[main]jwst_repo} + romancal: romancal[test] @ git+{[main]romancal_repo} + romanisim: romanisim[test] @ git+{[main]romanisim_repo} + specutils: specutils[test] @ git+{[main]specutils_repo} + dkist: dkist[tests] @ git+{[main]dkist_repo} + ndcube: ndcube[dev] @ git+{[main]ndcube_repo} + dev: -r requirements-dev.txt commands_pre = - dev: pip install -r requirements-dev.txt -U --upgrade-strategy eager - pip freeze + {list_dependencies_command} + jwst: git clone -n --depth=1 --filter=blob:none {[main]jwst_repo} target_repo + romancal: git clone -n --depth=1 --filter=blob:none {[main]romancal_repo} target_repo + romanisim: git clone -n --depth=1 --filter=blob:none {[main]romanisim_repo} target_repo + specutils: git clone -n --depth=1 --filter=blob:none {[main]specutils_repo} target_repo + dkist: git clone -n --depth=1 --filter=blob:none {[main]dkist_repo} target_repo + ndcube: git clone -n --depth=1 --filter=blob:none {[main]ndcube_repo} target_repo + + jwst,romancal,romanisim,specutils,dkist,ndcube: git --git-dir={env_tmp_dir}/target_repo/.git checkout HEAD pyproject.toml commands = pytest \ - warnings: -W error \ + jwst,romancal,romanisim,specutils,dkist,ndcube: --config-file={env_tmp_dir}/pyproject.toml --pyargs \ + jwst: jwst --ignore-glob=timeconversion --ignore-glob=associations --ignore-glob=scripts --show-capture=no \ + romancal: romancal \ + romanisim: romanisim \ + specutils: specutils \ + dkist: dkist --benchmark-skip \ + ndcube: ndcube \ xdist: -n auto \ - pyargs: {toxinidir}/docs --pyargs gwcs \ - jwst: --pyargs jwst --ignore-glob=timeconversion --ignore-glob=associations --ignore-glob=scripts\ - romancal: --pyargs romancal \ cov: --cov=. --cov-config=pyproject.toml --cov-report=term-missing --cov-report=xml \ {posargs} - -[testenv:build-docs] -description = invoke sphinx-build to build the HTML docs -extras = docs -commands = - sphinx-build -W docs docs/_build - -[testenv:build-dist] -description = build wheel and sdist -skip_install = true -deps = - build -commands = - python -m build .