diff --git a/.gitignore b/.gitignore index 0f077b45..395a25e2 100644 --- a/.gitignore +++ b/.gitignore @@ -15,6 +15,7 @@ __pycache__ htmlcov .coverage MANIFEST +.pytest_cache # Sphinx docs/api @@ -45,4 +46,4 @@ distribute-*.tar.gz # Mac OSX .DS_Store -.idea \ No newline at end of file +.idea diff --git a/.rtd-environment.yml b/.rtd-environment.yml new file mode 100644 index 00000000..cd08184f --- /dev/null +++ b/.rtd-environment.yml @@ -0,0 +1,12 @@ + +dependencies: + - python>=3 + - numpy + - pip: + - astropy + - sphinx_astropy + - sphinx-automodapi + - stsci_rtd_theme + - asdf + - sphinx_rtd_theme + diff --git a/.travis.yml b/.travis.yml index 89ad6d19..6733f6f3 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,4 +1,4 @@ -language: python +language: c # Setting sudo to false opts in to Travis-CI container-based builds. sudo: false @@ -12,34 +12,36 @@ addons: - texlive-latex-extra - dvipng -python: - - 3.5 - - 3.6 - env: global: # The following versions are the 'default' for tests, unless # overidden underneath. They are defined here in order to save having # to repeat them for all configurations. + - PYTHON_VERSION=3.6 - NUMPY_VERSION=stable - - ASTROPY_VERSION=development + #- ASTROPY_VERSION=development - MAIN_CMD='python setup.py' - CONDA_DEPENDENCIES='scipy' - ASDF_GIT='git+https://github.com/spacetelescope/asdf.git#egg=asdf' - - PIP_DEPENDENCIES="pytest-astropy $ASDF_GIT" + - ASTROPY_GIT='git+https://github.com/astropy/astropy.git@v3.1.x' + - PIP_DEPENDENCIES="$ASTROPY_GIT pytest-astropy asdf" - ASTROPY_USE_SYSTEM_PYTEST=1 + - PYTEST_VERSION=3.6 matrix: - SETUP_CMD='egg_info' - SETUP_CMD='test' + - PYTHON_VERSION=3.5 SETUP_CMD='test' + - PYTHON_VERSION=3.6 SETUP_CMD='test' + - PYTHON_VERSION=3.7 SETUP_CMD='test' matrix: # Don't wait for allowed failures fast_finish: true - + include: - + # Do a coverage test. - python: 3.5 env: SETUP_CMD='test --coverage --open-files --remote-data' @@ -49,19 +51,39 @@ matrix: - python: 3.6 env: SETUP_CMD='build_docs -w' - # Numpy - python: 3.6 - env: NUMPY_VERSION=1.12 SETUP_CMD="test" - - python: 3.6 - env: NUMPY_VERSION=1.13 SETUP_CMD="test" + env: NUMPY_VERSION=1.14 SETUP_CMD="test" # Do a PEP8 test with flake8 + - python: 3.6 + env: MAIN_CMD='flake8 count --select=F, E101, E111, E112, E113, E401, E402, E711, E722 --max-line-length=110' + SETUP_CMD='gwcs" + - python: 3.5 - env: MAIN_CMD='flake8 gwcs --count' SETUP_CMD='' + env: MAIN_CMD='flake8 count --max-line-length=110' SETUP_CMD='gwcs' + + - python: 3.6 + env: NUMPY_VERSION=dev SETUP_CMD='test' + + - python: 3.6 + env: ASTROPY_VERSION=development SETUP_CMD='test' + + - python: 3.6 + env: PIP_DEPENDENCIES="$ASTROPY_GIT $ASDF_GIT" SETUP_CMD='test' allow_failures: - python: 3.5 - env: MAIN_CMD='flake8 gwcs --count' SETUP_CMD='' + env: MAIN_CMD='flake8 count --max-line-length=110' SETUP_CMD='gwcs' + + - python: 3.6 + env: NUMPY_VERSION=dev SETUP_CMD='test' + + - python: 3.6 + env: ASTROPY_VERSION=development SETUP_CMD='test' + + - python: 3.6 + env: PIP_DEPENDENCIES="$ASTROPY_GIT $ASDF_GIT" SETUP_CMD='test' + install: - git clone git://github.com/astropy/ci-helpers.git diff --git a/CHANGES.rst b/CHANGES.rst index 5ca360af..c66f9a56 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,3 +1,20 @@ +0.10.0 (12/20/2018) +------------------- + +New Features +^^^^^^^^^^^^ + +- Initializing a ``WCS`` object with a ``pipeline`` list now keeps + the complete ``CoordinateFrame`` objects in the ``WCS.pipeline``. + The effect is that a ``WCS`` object can now be initialized with + a ``pipeline`` from a different ``WCS`` object. [#174] + +- Implement support for astropy APE 14 + (https://doi.org/10.5281/zenodo.1188875). [#146] + +- Added a ``wcs_from_[points`` function which creates a WCS object + two matching sets of points ``(x,y)`` and ``(ra, dec)``. [#42] + 0.9.0 (2018-05-23) ------------------ @@ -16,8 +33,10 @@ New Features - Add a ``StokesFrame`` which converts from 'I', 'Q', 'U', 'V' to 0-3. [#133] -- Support serialising the base ``CoordinateFrame`` class to asdf, by making - a specific tag and schema for ``Frame2D`` [#150] +- Support serializing the base ``CoordinateFrame`` class to asdf, by making + a specific tag and schema for ``Frame2D``. [#150] + +- Generalized the footrpint calculation to all output axes. [#167] API Changes @@ -26,6 +45,8 @@ API Changes - The argument ``output="numerical_plus"`` was replaced by a bool argument ``with_units``. [#156] +- Added a new flag ``axis_type`` to the footprint method. It controls what + type of footprint to calculate. [#167] Bug Fixes ^^^^^^^^^ diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 00000000..b3df6cff --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,5 @@ +Please open a new issue or new pull request for bugs, feedback, or new features you would like to see. If there is an issue you would like to work on, please leave a comment and we will be happy to assist. New contributions and contributors are very welcome! + +New to github or open source projects? If you are unsure about where to start or haven't used github before, please feel free to contact the package maintainers. + +Feedback and feature requests? Is there something missing you would like to see? Please open an issue or send an email to the maintainers. This package follows the Spacetelescope [Code of Conduct](CODE_OF_CONDUCT.md) strives to provide a welcoming community to all of our users and contributors. diff --git a/README.rst b/README.rst index 53e37b44..06fb9dd0 100644 --- a/README.rst +++ b/README.rst @@ -1,12 +1,34 @@ -GWCS -==== +.. + GWCS - Generalized World Coordinate System + ========================================== -Generalized World Coordinate System (GWCS) is an `Astropy`_ affiliated package providing tools for managing the World Coordinate System of astromnomical data. +.. raw:: html +

GWCS - Generalized World Coordinate System

+

+ Documentation Status + Build Status + Coverage Status + license + stsci + astropy +

+ + +Generalized World Coordinate System (GWCS) is an `Astropy`_ affiliated package providing tools for managing the World Coordinate System of astronomical data. + +GWCS takes a general approach to the problem of expressing transformations between pixel and world coordinates. It supports a data model which includes the entire transformation pipeline from input coordinates (detector by default) to world coordinates. It is tightly integrated with `Astropy`_. + +- Transforms are instances of ``astropy.Model``. They can be chained, joined or combined with arithmetic operators using the flexible framework of compound models in `astropy.modeling`_. +- Celestial coordinates are instances of ``astropy.SkyCoord`` and are transformed to other standard celestial frames using `astropy.coordinates`_. +- Time coordinates are represented by ``astropy.Time`` and can be further manipulated using the tools in `astropy.time`_ +- Spectral coordinates are ``astropy.Quantity`` objects and can be converted to other units using the tools in `astropy.units`_. + +For complete features and usage examples see the `documentation`_ site. Note ---- -Beginning with the next version (0.9), GWCS requires Python 3.5 and above. +Beginning with version 0.9 GWCS requires Python 3.5 and above. Installation @@ -22,15 +44,17 @@ To clone from github and install the master branch:: cd gwcs python setup.py install -Status ------- + +Contributing Code, Documentation, or Feedback +--------------------------------------------- -.. image:: https://travis-ci.org/spacetelescope/gwcs.svg - :target: https://travis-ci.org/spacetelescope/gwcs - :alt: GWCS's Travis CI Status -.. image:: https://coveralls.io/repos/github/spacetelescope/gwcs/badge.svg?branch=master - :target: https://coveralls.io/github/spacetelescope/gwcs?branch=master - :alt: Test Coverage Status +We welcome feedback and contributions to the project. Contributions of +code, documentation, or general feedback are all appreciated. Please +follow the `contributing guidelines `__ to submit an +issue or a pull request. + +We strive to provide a welcoming community to all of our users by +abiding to the `Code of Conduct `__. Citing GWCS @@ -41,9 +65,10 @@ Citing GWCS If you use GWCS, please cite the package via its Zenodo record. -License -------- - -GWCS is licensed under a 3-clause BSD style license. - .. _Astropy: http://www.astropy.org/ + +.. _astropy.time: http://docs.astropy.org/en/stable/time/ +.. _astropy.modeling: http://docs.astropy.org/en/stable/modeling/ +.. _astropy.units: http://docs.astropy.org/en/stable/units/ +.. _astropy.coordinates: http://docs.astropy.org/en/stable/coordinates/ +.. _documentation: http://gwcs.readthedocs.org/en/latest/ diff --git a/ah_bootstrap.py b/ah_bootstrap.py index e55254ee..ece5aa98 100644 --- a/ah_bootstrap.py +++ b/ah_bootstrap.py @@ -38,7 +38,6 @@ import contextlib import errno -import imp import io import locale import os @@ -46,55 +45,35 @@ import subprocess as sp import sys +__minimum_python_version__ = (3, 5) + +if sys.version_info < __minimum_python_version__: + print("ERROR: Python {} or later is required by astropy-helpers".format( + __minimum_python_version__)) + sys.exit(1) + try: from ConfigParser import ConfigParser, RawConfigParser except ImportError: from configparser import ConfigParser, RawConfigParser -if sys.version_info[0] < 3: - _str_types = (str, unicode) - _text_type = unicode - PY3 = False -else: - _str_types = (str, bytes) - _text_type = str - PY3 = True +_str_types = (str, bytes) # What follows are several import statements meant to deal with install-time # issues with either missing or misbehaving pacakges (including making sure # setuptools itself is installed): +# Check that setuptools 1.0 or later is present +from distutils.version import LooseVersion -# Some pre-setuptools checks to ensure that either distribute or setuptools >= -# 0.7 is used (over pre-distribute setuptools) if it is available on the path; -# otherwise the latest setuptools will be downloaded and bootstrapped with -# ``ez_setup.py``. This used to be included in a separate file called -# setuptools_bootstrap.py; but it was combined into ah_bootstrap.py try: - import pkg_resources - _setuptools_req = pkg_resources.Requirement.parse('setuptools>=0.7') - # This may raise a DistributionNotFound in which case no version of - # setuptools or distribute is properly installed - _setuptools = pkg_resources.get_distribution('setuptools') - if _setuptools not in _setuptools_req: - # Older version of setuptools; check if we have distribute; again if - # this results in DistributionNotFound we want to give up - _distribute = pkg_resources.get_distribution('distribute') - if _setuptools != _distribute: - # It's possible on some pathological systems to have an old version - # of setuptools and distribute on sys.path simultaneously; make - # sure distribute is the one that's used - sys.path.insert(1, _distribute.location) - _distribute.activate() - imp.reload(pkg_resources) -except: - # There are several types of exceptions that can occur here; if all else - # fails bootstrap and use the bootstrapped version - from ez_setup import use_setuptools - use_setuptools() - + import setuptools + assert LooseVersion(setuptools.__version__) >= LooseVersion('1.0') +except (ImportError, AssertionError): + print("ERROR: setuptools 1.0 or later is required by astropy-helpers") + sys.exit(1) # typing as a dependency for 1.6.1+ Sphinx causes issues when imported after # initializing submodule with ah_boostrap.py @@ -150,11 +129,7 @@ # TODO: Maybe enable checking for a specific version of astropy_helpers? DIST_NAME = 'astropy-helpers' PACKAGE_NAME = 'astropy_helpers' - -if PY3: - UPPER_VERSION_EXCLUSIVE = None -else: - UPPER_VERSION_EXCLUSIVE = '3' +UPPER_VERSION_EXCLUSIVE = None # Defaults for other options DOWNLOAD_IF_NEEDED = True @@ -186,7 +161,7 @@ def __init__(self, path=None, index_url=None, use_git=None, offline=None, if not (isinstance(path, _str_types) or path is False): raise TypeError('path must be a string or False') - if PY3 and not isinstance(path, _text_type): + if not isinstance(path, str): fs_encoding = sys.getfilesystemencoding() path = path.decode(fs_encoding) # path to unicode @@ -829,9 +804,9 @@ def run_cmd(cmd): stdio_encoding = 'latin1' # Unlikely to fail at this point but even then let's be flexible - if not isinstance(stdout, _text_type): + if not isinstance(stdout, str): stdout = stdout.decode(stdio_encoding, 'replace') - if not isinstance(stderr, _text_type): + if not isinstance(stderr, str): stderr = stderr.decode(stdio_encoding, 'replace') return (p.returncode, stdout, stderr) diff --git a/astropy_helpers b/astropy_helpers index 5e30a46d..9f82aac6 160000 --- a/astropy_helpers +++ b/astropy_helpers @@ -1 +1 @@ -Subproject commit 5e30a46dd0146a241f5c7ca8de0cc0144455b2e1 +Subproject commit 9f82aac6c2141b425e2d639560f7260189d90b54 diff --git a/gwcs/conftest.py b/conftest.py similarity index 81% rename from gwcs/conftest.py rename to conftest.py index e9dc3d04..29ffc17e 100644 --- a/gwcs/conftest.py +++ b/conftest.py @@ -4,11 +4,12 @@ # by importing them here in conftest.py they are discoverable by py.test # no matter how it is invoked within the source tree. -from astropy.tests.pytest_plugins import * +from astropy.tests.helper import enable_deprecations_as_exceptions +from astropy.tests.plugins.display import PYTEST_HEADER_MODULES, TESTED_VERSIONS # Uncomment the following line to treat all DeprecationWarnings as # exceptions -enable_deprecations_as_exceptions() +#enable_deprecations_as_exceptions() try: PYTEST_HEADER_MODULES['Astropy'] = 'astropy' diff --git a/convert_schemas.py b/convert_schemas.py new file mode 100644 index 00000000..feb5d327 --- /dev/null +++ b/convert_schemas.py @@ -0,0 +1,411 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +# -*- coding: utf-8 -*- + +from collections import OrderedDict +import io +import json +import os +import sys +import textwrap + +import yaml + + +def write_if_different(filename, data): + """ Write ``data`` to ``filename``, if the content of the file is different. + + Parameters + ---------- + filename : str + The file name to be written to. + data : bytes + The data to be written to `filename`. + """ + if not os.path.exists(os.path.dirname(filename)): + os.makedirs(os.path.dirname(filename)) + + if os.path.exists(filename): + with open(filename, 'rb') as fd: + original_data = fd.read() + else: + original_data = None + + if original_data != data: + print("Converting schema {0}".format( + os.path.basename(filename))) + with open(filename, 'wb') as fd: + fd.write(data) + + +def write_header(o, content, level): + """ + Write a reStructuredText header to the file. + + Parameters + ---------- + o : output stream + + content : str + The content of the header + + level : int + The level of the header + """ + levels = '=-~^.' + if level >= len(levels): + o.write('**{0}**\n\n'.format(content)) + else: + o.write(content) + o.write('\n') + o.write(levels[level] * len(content)) + o.write('\n\n') + + +def format_range(var_middle, var_end, minimum, maximum, + exclusiveMinimum, exclusiveMaximum): + """ + Formats an mathematical description of a range, for example, ``0 ≤ + x ≤ 2``. + + Parameters + ---------- + var_middle : str or None + The string to put in the middle of an expression, such as + the ``x`` in ``0 ≤ x ≤ 2``. + + var_end : str or None + The string to put at one end of a single comparision, such as + the ``x`` in ``x ≤ 0``. + + minimum : number + The minimum value. + + maximum : number + The maximum value. + + exclusiveMinimum : bool + If `True`, the range excludes the minimum value. + + exclusiveMaximum : bool + If `True`, the range excludes the maximum value + + Returns + ------- + expr : str + The formatted range expression + """ + if minimum is not None and maximum is not None: + part = '{0} '.format(minimum) + if exclusiveMinimum: + part += '<' + else: + part += '≤' + part += ' {0} '.format(var_middle) + if exclusiveMaximum: + part += '<' + else: + part += '≤' + part += ' {0}'.format(maximum) + elif minimum is not None: + if var_end is not None: + part = '{0} '.format(var_end) + else: + part = '' + if exclusiveMinimum: + part += '> {0}'.format(minimum) + else: + part += '≥ {0}'.format(minimum) + elif maximum is not None: + if var_end is not None: + part = '{0} '.format(var_end) + else: + part = '' + if exclusiveMaximum: + part += '< {0}'.format(maximum) + else: + part += '≤ {0}'.format(maximum) + else: + return None + return part + + +def format_type(schema, root): + """ + Creates an English/mathematical description of a schema fragment. + + Parameters + ---------- + schema : JSON schema fragment + + root : str + The JSON path to the schema fragment. + """ + if 'anyOf' in schema: + return ' :soft:`or` '.join( + format_type(x, root) for x in schema['anyOf']) + + elif 'allOf' in schema: + return ' :soft:`and` '.join( + format_type(x, root) for x in schema['allOf']) + + elif '$ref' in schema: + ref = schema['$ref'] + if ref.startswith('#/'): + return ':ref:`{0} <{1}/{2}>`'.format(ref[2:], root, ref[2:]) + else: + basename = os.path.basename(ref) + if "tag:stsci.edu:asdf" in ref or "tag:astropy.org:astropy" in ref: + return '`{0} <{1}>`'.format(basename, ref) + else: + return ':doc:`{0} <{1}>`'.format(basename, ref) + + else: + type = schema.get('type') + if isinstance(type, list): + parts = [' or '.join(type)] + + elif type is None: + parts = ['any'] + + else: + parts = [type] + + if type == 'string': + range = format_range('*len*', '*len*', schema.get('minLength'), + schema.get('maxLength'), False, False) + if range is not None or 'pattern' in schema or 'format' in schema: + parts.append('(') + if range is not None: + parts.append(range) + if 'pattern' in schema: + pattern = schema['pattern'].encode('unicode_escape') + pattern = pattern.decode('ascii') + parts.append(':soft:`regex` :regexp:`{0}`'.format(pattern)) + if 'format' in schema: + parts.append(':soft:`format` {0}'.format(schema['format'])) + parts.append(')') + + elif type in ('integer', 'number'): + range = format_range('*x*', '', schema.get('minimum'), + schema.get('maximum'), + schema.get('exclusiveMinimum'), + schema.get('exclusiveMaximum')) + if range is not None: + parts.append(range) + # TODO: multipleOf + + elif type == 'object': + range = format_range('*len*', '*len*', schema.get('minProperties'), + schema.get('maxProperties'), False, False) + if range is not None: + parts.append(range) + # TODO: Dependencies + # TODO: Pattern properties + + elif type == 'array': + items = schema.get('items') + if schema.get('items') and isinstance(items, dict): + if schema.get('uniqueItems'): + parts.append(':soft:`of unique`') + else: + parts.append(':soft:`of`') + parts.append('(') + parts.append(format_type(items, root)) + parts.append(')') + range = format_range('*len*', '*len*', schema.get('minItems'), + schema.get('maxItems'), False, False) + if range is not None: + parts.append(range) + + if 'enum' in schema: + parts.append(':soft:`from`') + parts.append(json.dumps(schema['enum'])) + + return ' '.join(parts) + + +def reindent(content, indent): + """ + Reindent a string to the given number of spaces. + """ + content = textwrap.dedent(content) + lines = [] + for line in content.split('\n'): + lines.append(indent + line) + return '\n'.join(lines) + + +def recurse(o, name, schema, path, level, required=False): + """ + Convert a schema fragment to reStructuredText. + + Parameters + ---------- + o : output stream + + name : str + Name of the entry + + schema : schema fragment + + path : list of str + Path to schema fragment + + level : int + Indentation level + + required : bool + If `True` the entry is required by the schema and will be + documented as such. + """ + indent = ' ' * max(level, 0) + o.write('\n\n') + o.write(indent) + o.write('.. _{0}:\n\n'.format(os.path.join(*path))) + if level == 0: + write_header(o, name, level) + else: + if name != 'items': + o.write(indent) + o.write(':entry:`{0}`\n\n'.format(name)) + + o.write(indent) + if path[0].startswith("tag:stsci.edu:asdf"): + o.write(format_type(schema, path[0])) + else: + o.write(":soft:`Type:` ") + o.write(format_type(schema, path[0])) + o.write('.') + if required: + o.write(' Required.') + o.write('\n\n') + + o.write(reindent(schema.get('title', ''), indent)) + o.write('\n\n') + + o.write(reindent(schema.get('description', ''), indent)) + o.write('\n\n') + + if 'default' in schema: + o.write(indent) + o.write(':soft:`Default:` {0}'.format( + json.dumps(schema['default']))) + o.write('\n\n') + + if 'definitions' in schema: + o.write(indent) + o.write(":category:`Definitions:`\n\n") + for key, val in schema['definitions'].items(): + recurse(o, key, val, path + ['definitions', key], level + 1) + + if 'anyOf' in schema and len(schema['anyOf']) > 1: + o.write(indent) + o.write(':category:`Any of:`\n\n') + for i, subschema in enumerate(schema['anyOf']): + recurse(o, '—', subschema, path + ['anyOf', str(i)], level + 1) + + elif 'allOf' in schema and len(schema['allOf']) > 1: + o.write(indent) + o.write(':category:`All of:`\n\n') + for i, subschema in enumerate(schema['allOf']): + recurse(o, i, subschema, path + ['allOf', str(i)], level + 1) + + if schema.get('type') == 'object': + o.write(indent) + o.write(':category:`Properties:`\n\n') + for key, val in schema.get('properties', {}).items(): + recurse(o, key, val, path + ['properties', key], level + 1, + key in schema.get('required', [])) + + elif schema.get('type') == 'array': + o.write(indent) + o.write(':category:`Items:`\n\n') + items = schema.get('items') + if isinstance(items, dict): + recurse(o, 'items', items, path + ['items'], level + 1) + elif isinstance(items, list): + for i, val in enumerate(items): + name = 'index[{0}]'.format(i) + recurse(o, name, val, path + [str(i)], level + 1) + + if 'examples' in schema: + o.write(indent) + o.write(":category:`Examples:`\n\n") + for description, example in schema['examples']: + o.write(reindent(description + "::\n\n", indent)) + o.write(reindent(example, indent + ' ')) + o.write('\n\n') + + +def convert_schema_to_rst(src, dst): + """ + Convert a YAML schema to reStructuredText. + """ + with open(src, 'rb') as fd: + schema = yaml.safe_load(fd) + with open(src, 'rb') as fd: + yaml_content = fd.read() + + o = io.StringIO() + + id = schema.get('id', '#') + name = os.path.basename(src[:-5]) + if 'title' in schema: + name += ': ' + schema['title'].strip() + recurse(o, name, schema, [id], 0) + + #o.write(".. only:: html\n\n :download:`Original schema in YAML <{0}>`\n". + #os.path.basename(src))) + + write_if_different(dst, yaml_content) + write_if_different(dst[:-5] + ".rst", o.getvalue().encode('utf-8')) + + +def construct_mapping(self, node, deep=False): + """ + Make sure the properties are written out in the same order as the + original file. + """ + if not isinstance(node, yaml.MappingNode): + raise yaml.constructor.ConstructorError(None, None, + "expected a mapping node, but found %s" % node.id, + node.start_mark) + mapping = OrderedDict() + for key_node, value_node in node.value: + key = self.construct_object(key_node, deep=deep) + try: + hash(key) + except TypeError as exc: + raise yaml.constructor.ConstructorError( + "while constructing a mapping", node.start_mark, + "found unacceptable key (%s)" % exc, key_node.start_mark) + value = self.construct_object(value_node, deep=deep) + mapping[key] = value + return mapping + + +yaml.SafeLoader.add_constructor( + 'tag:yaml.org,2002:map', construct_mapping) + + +def main(src, dst): + for root, dirs, files in os.walk(src): + for fname in files: + if not fname.endswith(".yaml"): + continue + src_path = os.path.join(root, fname) + dst_path = os.path.join( + dst, os.path.relpath(src_path, src)) + + convert_schema_to_rst(src_path, dst_path) + + +def decode_filename(fname): + return fname + + +if __name__ == '__main__': + src = decode_filename(sys.argv[-2]) + dst = decode_filename(sys.argv[-1]) + + sys.exit(main(src, dst)) diff --git a/docs/conf.py b/docs/conf.py index bc10969b..df0498ef 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -140,8 +140,8 @@ # Grouping the document tree into LaTeX files. List of tuples # (source start file, target name, title, author, documentclass [howto/manual]). -latex_documents = [('index', project + '.tex', project + u' Documentation', - author, 'manual')] +#latex_documents = [('index', project + '.tex', project + u' Documentation', +# author, 'manual')] # -- Options for manual page output -------------------------------------------- @@ -167,5 +167,7 @@ edit_on_github_source_root = "" edit_on_github_doc_root = "docs" -#sys.path.insert(0, os.path.join(os.path.dirname('__file__'), 'sphinxext')) -#extensions += ['example'] + +sys.path.insert(0, os.path.join(os.path.dirname('__file__'), 'sphinxext')) +extensions += ['category'] + diff --git a/docs/gwcs/ifu.rst b/docs/gwcs/ifu.rst new file mode 100644 index 00000000..caf2bd06 --- /dev/null +++ b/docs/gwcs/ifu.rst @@ -0,0 +1,122 @@ +An IFU Example - managing a discontiguous WCS +============================================= + +An IFU image represents the projection of several slices on a detector. +Between the slices there are pixels which don't belong to any slice. +In general each slice has a unique WCS transform. +There are two ways to represent this kind of transforms in GWCS depending on +the way the instrument is calibrated and the available information. + + +Using a pixel to slice mapping +------------------------------ + +In this case a pixel map associating each pixel with a slice label (number or string) +is available. The image below represents the projection of the slits of an IFU on a detector +with a size (500, 1000). Slices are labeled from 1 to 6, while label 0 is reserved for pixels +between the slices. + +.. image:: ifu-regions.png + +There are several models in GWCS which are useful in creating a WCS. + +Given (x, y) pixel indices, `~gwcs.selector.LabelMapperArray` returns labels (int or str) +associated with these indices. `~gwcs.selector.RegionsSelector` +maps labels with transforms. It uses the `~gwcs.selector.LabelMapperArray` to map +these transforms to pixel indices. + +A step by step example of constructing the WCS for an IFU with 6 slits follows. + +First, import the usual packages. + + >>> import numpy as np + >>> from astropy.modeling import models + >>> from astropy import coordinates as coord + >>> from astropy import units as u + >>> from gwcs import wcs, selector + >>> from gwcs import coordinate_frames as cf + +The output frame is common for all slits and is a composite frame with two subframes, +`~gwcs.coordinate_frames.CelestialFrame` and `~gwcs.coordinate_frames.SpectralFrame`. + + >>> sky_frame = cf.CelestialFrame(name='icrs', reference_frame=coord.ICRS(), axes_order=(0, 2)) + >>> spec_frame = cf.SpectralFrame(name='wave', unit=(u.micron,), axes_order=(1,), axes_names=('lambda',)) + >>> cframe = cf.CompositeFrame([sky_frame, spec_frame], name='world') + >>> det = cf.Frame2D(name='detector') + +All slices have the same input and output frames, however each slices has a different model transforming +from pixels to world coordinates (RA, lambda, dec). For the sake of brevity this example uses a simple +shift transform for each slice. Detailed examples of how to create more realistic transforms +are available in :ref:`imaging_example`. + + >>> transforms = {} + >>> for i in range(1, 7): + ... transforms[i] = models.Mapping([0, 0, 1]) | models.Shift(i * 0.1) & models.Shift(i * 0.2) & models.Scale(i * 0.1) + +One way to initialize `~gwcs.selector.LabelMapperArray` is to pass it the shape of the array and the vertices +of each slit on the detector {label: vertices} see :meth: `~gwcs.selector.LabelMapperArray.from_vertices`. +In this example the mask is an array with the size of the detector where each item in the array +corresponds to a pixel on the detector and its value is the slice number (label) this pixel +belongs to. + +Assuming the array is stored in +`ASDF `__ format, create the mask: + +.. doctest-skip-all + + >>> import asdf + >>> f = asdf.open('mask.asdf') + >>> data = f.tree['mask'] + >>> mask = selector.LabelMapperArray(data) + +Create the pixel to world transform for the entire IFU: + + >>> regions_transform = selector.RegionsSelector(inputs=['x','y'], + ... outputs=['ra', 'dec', 'lam'], + ... selector=transforms, + ... label_mapper=mask, + ... undefined_transform_value=np.nan) + +The WCS object now can evaluate simultaneously the transforms of all slices. + + >>> wifu = wcs.WCS(forward_transform=regions_transform, output_frame=cframe, input_frame=det) + >>> y, x = mask.mapper.shape + >>> y, x = np.mgrid[:y, :x] + >>> r, d, l = wifu(x, y) + +or of single slices. + +The :meth:`~gwcs.selector.RegionsSelector.set_input` method returns the forward_transform for +a specific label. + + + >>> wifu.forward_transform.set_input(4)(1, 2) + (1.4, 1.8, 0.8) + +Custom model storing transforms in a dictionary +----------------------------------------------- + +In case a pixel to slice mapping is not available, one can write a custom mdoel +storing transforms in a dictionary. The model would look like this: + +.. code:: + + from astropy.modeling.core import Model + from astropy.modeling.parameters import Parameter + + + class CustomModel(Model): + + inputs = ('label', 'x', 'y') + outputs = ('xout', 'yout') + + + def __init__(self, labels, transforms): + super().__init__() + self.labels = labels + self.models = models + + def evaluate(self, label, x, y): + index = self.labels.index(label) + return self.models[index](x, y) + diff --git a/docs/gwcs/imaging_with_distortion.rst b/docs/gwcs/imaging_with_distortion.rst new file mode 100644 index 00000000..47b950f7 --- /dev/null +++ b/docs/gwcs/imaging_with_distortion.rst @@ -0,0 +1,78 @@ +.. _imaging_example: + +Adding distortion to the imaging example +======================================== + +Let's expand the WCS created in :ref:`getting-started` by adding a polynomial +distortion correction. + +Because the polynomial models in `~astropy.modeling` do not support units yet, +this example will use transforms without units. At the end the units +associated with the output frame are used to create a `~astropy.coordinates.SkyCoord` object. + +The imaging example without units: + + >>> import numpy as np + >>> from astropy.modeling import models + >>> from astropy import coordinates as coord + >>> from astropy import units as u + >>> from gwcs import wcs + >>> from gwcs import coordinate_frames as cf + + >>> 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) + >>> rotation.inverse = models.AffineTransformation2D(np.linalg.inv(matrix)) + >>> 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"), + ... 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) + ... ] + >>> wcsobj = wcs.WCS(pipeline) + >>> print(wcsobj) + From Transform + -------- ---------------- + detector linear_transform + icrs None + +First create distortion corrections represented by a polynomial +model of fourth degree. The example uses the astropy `~astropy.modeling.polynomial.Polynomial2D` +and `~astropy.modeling.mappings.Mapping` models. + + + >>> poly_x = models.Polynomial2D(4) + >>> poly_x.parameters = np.arange(15) * .1 + >>> poly_y = models.Polynomial2D(4) + >>> poly_y.parameters = np.arange(15) * .2 + >>> distortion = models.Mapping((0, 1, 0, 1)) | poly_x & poly_y + >>> distortion.name = "distortion" + +Create an intermediate frame for distortion free coordinates. + + + >>> undistorted_frame = cf.Frame2D(name="undistorted_frame", unit=(u.pix, u.pix), + ... axes_names=("undist_x", "undist_y")) + + +Using the example in :ref:`getting-started`, add the distortion correction to +the WCS pipeline and initialize the WCS. + + >>> pipeline = [(detector_frame, distortion), + ... (undistorted_frame, det2sky), + ... (sky_frame, None) + ... ] + >>> wcsobj = wcs.WCS(pipeline) + >>> print(wcsobj) + From Transform + ----------------- ---------------- + detector distortion + undistorted_frame linear_transform + icrs None + diff --git a/docs/gwcs/points_to_wcs.rst b/docs/gwcs/points_to_wcs.rst new file mode 100644 index 00000000..5c265b0e --- /dev/null +++ b/docs/gwcs/points_to_wcs.rst @@ -0,0 +1,83 @@ + +.. _wcs_from_points_example: + +Fitting a WCS to input pixels & sky positions +============================================= + +Suppose we have an image where we have centroid positions for a number of sources, and we have matched these +positions to an external catalog to obtain (RA, Dec). If this data is missing or has inaccurate WCS information, +it is useful to fit or re-fit a GWCS object with this matched list of coordinate pairs to be able to transform +between pixel and sky. + +This example shows how to use the `~gwcs.wcstools.wcs_from_points` tool to fit a WCS to a matched set of +pixel and sky positions. Along with arrays of the (x,y) pixel position in the image and the matched sky coordinates, +the fiducial point for the projection must be supplied as a `~astropy.coordinates.SkyCoord` object. Additionally, +the projection type must be specified from the available projections in `~astropy.modeling.projections.projcode`. + +Geometric distortion can also be fit to the input coordinates - the distortion type (2D polynomial, chebyshev, legendre) and +the degree can be supplied to fit this component of the model. + +The following example will show how to fit a WCS, including a 4th degree 2D polynomial, to a set of input pixel positions of +sources in an image and their corresponding positions on the sky obtained from a catalog. + +Import the wcs_from_points function, + + >>> from gwcs.wcstools import wcs_from_points + +along with some useful general imports. + + >>> from astropy.coordinates import SkyCoord + >>> from astropy.io import ascii + >>> import astropy.units as u + >>> import numpy as np + +A collection of 20 matched coordinate pairs in x, y, RA, and Dec stored in two arrays, will be used to fit the WCS information. The function requires tuples of arrays. + + >>> xy = (np.array([2810.156, 2810.156, 650.236, 1820.927, 3425.779, 2750.369, + ... 212.422, 1146.91 , 27.055, 2100.888, 648.149, 22.212, + ... 2003.314, 727.098, 248.91 , 409.998, 1986.931, 128.925, + ... 1106.654, 1502.67 ]), + ... np.array([1670.347, 1670.347, 360.325, 165.663, 900.922, 700.148, + ... 1416.235, 1372.364, 398.823, 580.316, 317.952, 733.984, + ... 339.024, 234.29 , 1241.608, 293.545, 1794.522, 1365.706, + ... 583.135, 25.306])) + >>> radec = (np.array([246.75001315, 246.75001315, 246.72033646, 246.72303144, + ... 246.74164072, 246.73540614, 246.73379121, 246.73761455, + ... 246.7179495 , 246.73051123, 246.71970072, 246.7228646 , + ... 246.72647213, 246.7188386 , 246.7314031 , 246.71821002, + ... 246.74785534, 246.73265223, 246.72579817, 246.71943263]), + ... np.array([43.48690547, 43.48690547, 43.46792989, 43.48075238, + ... 43.49560501, 43.48903538, 43.46045875, 43.47030776, + ... 43.46132376, 43.48252763, 43.46802566, 43.46035331, + ... 43.48218262, 43.46908299, 43.46131665, 43.46560591, + ... 43.47791234, 43.45973025, 43.47208325, 43.47779988])) + + +We can now choose the reference point on the sky for the projection. This is passed in +as a `~astropy.coordinates.SkyCoord` object so that information about the celestial frame and units is given as well. +The input world coordinates are passed in as unitless arrays, and so are assumed to be of the same unit and frame +as the fiducial point. + + >>> proj_point = SkyCoord(246.7368408, 43.480712949, frame = 'icrs', unit = (u.deg,u.deg)) + +We can now call the function that returns a GWCS object corresponding to the best fit parameters +that relate the input pixels and sky coordinates with a TAN projection centered at the reference point +we specified, with a distortion model (degree 4 polynomial). This function will return a GWCS object that +can be used to transform between coordinate frames. + + >>> gwcs_obj = wcs_from_points(xy, radec, proj_point) + +This GWCS object contains parameters for a TAN projection, rotation, scale, skew and a polynomial fit to x and y +that represent the best-fit to the input coordinates. With WCS information associated with the data now, we can +easily work in both pixel and sky space, and transform between frames. + +The GWCS object, which by default when called executes for forward transformation, +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 + (246.72158004206716, 43.46075091731673) + diff --git a/docs/gwcs/pure_asdf.rst b/docs/gwcs/pure_asdf.rst new file mode 100644 index 00000000..26c644a4 --- /dev/null +++ b/docs/gwcs/pure_asdf.rst @@ -0,0 +1,137 @@ +.. _pure_asdf: + +Listing of ``imaging_wcs.asdf`` +=============================== + +Listing of ``imaging_wcs.asdf``:: + + + #ASDF 1.0.0 + #ASDF_STANDARD 1.2.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.2.0.dev1526} + history: + extensions: + - !core/extension_metadata-1.0.0 + extension_class: asdf.extension.BuiltinExtension + software: {name: asdf, version: 2.2.0.dev1526} + - !core/extension_metadata-1.0.0 + extension_class: astropy.io.misc.asdf.extension.AstropyExtension + software: {name: astropy, version: 3.2.dev23222} + - !core/extension_metadata-1.0.0 + extension_class: astropy.io.misc.asdf.extension.AstropyAsdfExtension + software: {name: astropy, version: 3.2.dev23222} + - !core/extension_metadata-1.0.0 + extension_class: gwcs.extension.GWCSExtension + software: {name: gwcs, version: 0.10.dev417} + wcs: ! + name: '' + steps: + - ! + frame: ! + axes_names: [x, y] + name: detector + unit: [!unit/unit-1.0.0 pixel, !unit/unit-1.0.0 pixel] + transform: !transform/compose-1.1.0 + forward: + - !transform/remap_axes-1.1.0 + mapping: [0, 1, 0, 1] + - !transform/concatenate-1.1.0 + forward: + - !transform/polynomial-1.1.0 + coefficients: !core/ndarray-1.0.0 + data: + - [0.0, 0.5, 0.6000000000000001, 0.7000000000000001, 0.8] + - [0.1, 0.9, 1.0, 1.1, 0.0] + - [0.2, 1.2000000000000002, 1.3, 0.0, 0.0] + - [0.30000000000000004, 1.4000000000000001, 0.0, 0.0, 0.0] + - [0.4, 0.0, 0.0, 0.0, 0.0] + datatype: float64 + shape: [5, 5] + - !transform/polynomial-1.1.0 + coefficients: !core/ndarray-1.0.0 + data: + - [0.0, 1.0, 1.2000000000000002, 1.4000000000000001, 1.6] + - [0.2, 1.8, 2.0, 2.2, 0.0] + - [0.4, 2.4000000000000004, 2.6, 0.0, 0.0] + - [0.6000000000000001, 2.8000000000000003, 0.0, 0.0, 0.0] + - [0.8, 0.0, 0.0, 0.0, 0.0] + datatype: float64 + shape: [5, 5] + - ! + frame: ! + axes_names: [undist_x, undist_y] + name: undistorted_frame + unit: [!unit/unit-1.0.0 pixel, !unit/unit-1.0.0 pixel] + transform: !transform/compose-1.1.0 + forward: + - !transform/compose-1.1.0 + forward: + - !transform/compose-1.1.0 + forward: + - !transform/concatenate-1.1.0 + forward: + - !transform/shift-1.2.0 {offset: -2048.0} + - !transform/shift-1.2.0 {offset: -1024.0} + - !transform/affine-1.2.0 + inverse: !transform/affine-1.2.0 + matrix: !core/ndarray-1.0.0 + data: + - [65488.318039522, 30828.31712434267] + - [26012.509548778366, -66838.34993781192] + datatype: float64 + shape: [2, 2] + translation: !core/ndarray-1.0.0 + data: [0.0, 0.0] + datatype: float64 + shape: [2] + matrix: !core/ndarray-1.0.0 + data: + - [1.290551569736e-05, 5.9525007864732e-06] + - [5.0226382102765e-06, -1.2644844123757e-05] + datatype: float64 + shape: [2, 2] + translation: !core/ndarray-1.0.0 + data: [0.0, 0.0] + datatype: float64 + shape: [2] + - !transform/gnomonic-1.1.0 {direction: pix2sky} + - !transform/rotate3d-1.2.0 {phi: 5.63056810618, psi: 180.0, theta: -72.05457184279} + inverse: !transform/compose-1.1.0 + forward: + - !transform/rotate3d-1.2.0 {direction: celestial2native, phi: 5.63056810618, + psi: 180.0, theta: -72.05457184279} + - !transform/compose-1.1.0 + forward: + - !transform/gnomonic-1.1.0 {direction: sky2pix} + - !transform/compose-1.1.0 + forward: + - !transform/affine-1.2.0 + matrix: !core/ndarray-1.0.0 + data: + - [65488.318039522, 30828.31712434267] + - [26012.509548778366, -66838.34993781192] + datatype: float64 + shape: [2, 2] + translation: !core/ndarray-1.0.0 + data: [0.0, 0.0] + datatype: float64 + shape: [2] + - !transform/concatenate-1.1.0 + forward: + - !transform/shift-1.2.0 {offset: 2048.0} + - !transform/shift-1.2.0 {offset: 1024.0} + name: linear_transform + - ! + frame: ! + axes_names: [lon, lat] + name: icrs + reference_frame: ! + frame_attributes: {} + unit: [!unit/unit-1.0.0 deg, !unit/unit-1.0.0 deg] + ... + + diff --git a/docs/gwcs/schemas/celestial_frame-1.0.0.rst b/docs/gwcs/schemas/celestial_frame-1.0.0.rst new file mode 100644 index 00000000..a389a625 --- /dev/null +++ b/docs/gwcs/schemas/celestial_frame-1.0.0.rst @@ -0,0 +1,78 @@ + + +.. _http://stsci.edu/schemas/gwcs/celestial_frame-1.0.0: + +celestial_frame-1.0.0: Represents a celestial frame. +==================================================== + +:soft:`Type:` object :soft:`and` :doc:`frame-1.0.0 `. + +Represents a celestial frame. + + + + +:category:`All of:` + + + + .. _http://stsci.edu/schemas/gwcs/celestial_frame-1.0.0/allOf/0: + + :entry:`0` + + :soft:`Type:` object. + + + + + + :category:`Properties:` + + + + .. _http://stsci.edu/schemas/gwcs/celestial_frame-1.0.0/allOf/0/properties/axes_names: + + :entry:`axes_names` + + :soft:`Type:` any. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/celestial_frame-1.0.0/allOf/0/properties/axes_order: + + :entry:`axes_order` + + :soft:`Type:` any. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/celestial_frame-1.0.0/allOf/0/properties/unit: + + :entry:`unit` + + :soft:`Type:` any. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/celestial_frame-1.0.0/allOf/1: + + :entry:`1` + + :soft:`Type:` :doc:`frame-1.0.0 `. + + + + + diff --git a/docs/gwcs/schemas/composite_frame-1.0.0.rst b/docs/gwcs/schemas/composite_frame-1.0.0.rst new file mode 100644 index 00000000..12964afa --- /dev/null +++ b/docs/gwcs/schemas/composite_frame-1.0.0.rst @@ -0,0 +1,14 @@ + + +.. _http://stsci.edu/schemas/gwcs/composite_frame-1.0.0: + +composite_frame-1.0.0: Represents a set of frames. +================================================== + +:soft:`Type:` object. + +Represents a set of frames. + + + + diff --git a/docs/gwcs/schemas/frame-1.0.0.rst b/docs/gwcs/schemas/frame-1.0.0.rst new file mode 100644 index 00000000..b0a20a01 --- /dev/null +++ b/docs/gwcs/schemas/frame-1.0.0.rst @@ -0,0 +1,197 @@ + + +.. _http://stsci.edu/schemas/gwcs/frame-1.0.0: + +frame-1.0.0: The base class of all coordinate frames. +===================================================== + +:soft:`Type:` object. + +The base class of all coordinate frames. + + +These objects are designed to be nested in arbitrary ways to build up +transformation pipelines out of a number of low-level pieces. + + +:category:`Properties:` + + + + .. _http://stsci.edu/schemas/gwcs/frame-1.0.0/properties/name: + + :entry:`name` + + :soft:`Type:` string. Required. + + + + A user-friendly name for the frame. + + + + + .. _http://stsci.edu/schemas/gwcs/frame-1.0.0/properties/axes_order: + + :entry:`axes_order` + + :soft:`Type:` array :soft:`of` ( integer ). + + + + The order of the axes. + + + :category:`Items:` + + + + .. _http://stsci.edu/schemas/gwcs/frame-1.0.0/properties/axes_order/items: + + :soft:`Type:` integer. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/frame-1.0.0/properties/axes_names: + + :entry:`axes_names` + + :soft:`Type:` array :soft:`of` ( string :soft:`or` null ). + + + + The name of each axis in this frame. + + + :category:`Items:` + + + + .. _http://stsci.edu/schemas/gwcs/frame-1.0.0/properties/axes_names/items: + + :soft:`Type:` string :soft:`or` null. + + + + + + :category:`Any of:` + + + + .. _http://stsci.edu/schemas/gwcs/frame-1.0.0/properties/axes_names/items/anyOf/0: + + :entry:`—` + + :soft:`Type:` string. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/frame-1.0.0/properties/axes_names/items/anyOf/1: + + :entry:`—` + + :soft:`Type:` null. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/frame-1.0.0/properties/reference_frame: + + :entry:`reference_frame` + + :soft:`Type:` `baseframe-1.0.0 `. + + + + The reference frame. + + + + + .. _http://stsci.edu/schemas/gwcs/frame-1.0.0/properties/unit: + + :entry:`unit` + + :soft:`Type:` array :soft:`of` ( `unit-1.0.0 ` ). + + + + Units for each axis. + + + :category:`Items:` + + + + .. _http://stsci.edu/schemas/gwcs/frame-1.0.0/properties/unit/items: + + :soft:`Type:` `unit-1.0.0 `. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/frame-1.0.0/properties/axis_physical_types: + + :entry:`axis_physical_types` + + :soft:`Type:` array :soft:`of` ( string ). + + + + An iterable of strings describing the physical type for each world axis. + These should be names from the VO UCD1+ controlled Vocabulary + (http://www.ivoa.net/documents/latest/UCDlist.html). + + + :category:`Items:` + + + + .. _http://stsci.edu/schemas/gwcs/frame-1.0.0/properties/axis_physical_types/items: + + :soft:`Type:` string. + + + + + +:category:`Examples:` + +A celestial frame in the ICRS reference frame. +:: + + ! + axes_names: [lon, lat] + name: CelestialFrame + reference_frame: ! + frame_attributes: {} + unit: [!unit/unit-1.0.0 deg, !unit/unit-1.0.0 deg] + + +A pixel frame in three dimensions +:: + + ! + axes_names: [raster position, slit position, wavelength] + axes_order: [0, 1, 2] + axes_type: [SPATIAL, SPATIAL, SPECTRAL] + name: pixel + naxes: 3 + unit: [!unit/unit-1.0.0 pixel, !unit/unit-1.0.0 pixel, !unit/unit-1.0.0 pixel] + + diff --git a/docs/gwcs/schemas/frame2d-1.0.0.rst b/docs/gwcs/schemas/frame2d-1.0.0.rst new file mode 100644 index 00000000..3209f042 --- /dev/null +++ b/docs/gwcs/schemas/frame2d-1.0.0.rst @@ -0,0 +1,89 @@ + + +.. _http://stsci.edu/schemas/gwcs/frame2d-1.0.0: + +frame2d-1.0.0: Represents a 2D frame. +===================================== + +:soft:`Type:` object :soft:`and` :doc:`frame-1.0.0 `. + +Represents a 2D frame. + + + + +:category:`All of:` + + + + .. _http://stsci.edu/schemas/gwcs/frame2d-1.0.0/allOf/0: + + :entry:`0` + + :soft:`Type:` object. + + + + + + :category:`Properties:` + + + + .. _http://stsci.edu/schemas/gwcs/frame2d-1.0.0/allOf/0/properties/axes_names: + + :entry:`axes_names` + + :soft:`Type:` any. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/frame2d-1.0.0/allOf/0/properties/axes_order: + + :entry:`axes_order` + + :soft:`Type:` any. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/frame2d-1.0.0/allOf/0/properties/unit: + + :entry:`unit` + + :soft:`Type:` any. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/frame2d-1.0.0/allOf/1: + + :entry:`1` + + :soft:`Type:` :doc:`frame-1.0.0 `. + + + + + +:category:`Examples:` + +A two dimensional spatial frame +:: + + ! + axes_names: [lon, lat] + name: Frame2D + unit: [!unit/unit-1.0.0 pixel, !unit/unit-1.0.0 pixel] + + diff --git a/docs/gwcs/schemas/index.rst b/docs/gwcs/schemas/index.rst new file mode 100644 index 00000000..af42777f --- /dev/null +++ b/docs/gwcs/schemas/index.rst @@ -0,0 +1,35 @@ +.. _gwcs-schemas: + +GWCS Schema Definitions +======================= + +WCS object +---------- + +.. toctree:: + :maxdepth: 1 + + wcs-1.0.0.rst + step-1.0.0.rst + +Coordinate Frames +----------------- + +.. toctree:: + :maxdepth: 1 + + celestial_frame-1.0.0.rst + frame-1.0.0.rst + spectral_frame-1.0.0.rst + frame2d-1.0.0.rst + temporal_frame-1.0.0.rst + composite_frame-1.0.0.rst + +Transforms +---------- + +.. toctree:: + :maxdepth: 1 + + label_mapper-1.0.0.rst + regions_selector-1.0.0.rst diff --git a/docs/gwcs/schemas/label_mapper-1.0.0.rst b/docs/gwcs/schemas/label_mapper-1.0.0.rst new file mode 100644 index 00000000..8b5626cd --- /dev/null +++ b/docs/gwcs/schemas/label_mapper-1.0.0.rst @@ -0,0 +1,340 @@ + + +.. _http://stsci.edu/schemas/gwcs/label_mapper-1.0.0: + +label_mapper-1.0.0: Represents a mapping from a coordinate value to a label. +============================================================================ + +:soft:`Type:` `transform-1.1.0 ` :soft:`and` object. + +Represents a mapping from a coordinate value to a label. + + +A label mapper instance maps inputs to a label. It is used together +with +[regions_selector](ref:http://stsci.edu/schemas/gwcs/regions_selector-1.0.0). The +[label_mapper](ref:http://stsci.edu/schemas/gwcs/label_mapper-1.0.0) +returns the label corresponding to given inputs. The +[regions_selector](ref:http://stsci.edu/schemas/gwcs/regions_selector-1.0.0) +returns the transform corresponding to this label. This maps inputs +(e.g. pixels on a detector) to transforms uniquely. + + +:category:`All of:` + + + + .. _http://stsci.edu/schemas/gwcs/label_mapper-1.0.0/allOf/0: + + :entry:`0` + + :soft:`Type:` `transform-1.1.0 `. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/label_mapper-1.0.0/allOf/1: + + :entry:`1` + + :soft:`Type:` object. + + + + + + :category:`Properties:` + + + + .. _http://stsci.edu/schemas/gwcs/label_mapper-1.0.0/allOf/1/properties/mapper: + + :entry:`mapper` + + :soft:`Type:` `ndarray-1.0.0 ` :soft:`or` `transform-1.1.0 ` :soft:`or` object. Required. + + + + A mapping of inputs to labels. + In the general case this is a `astropy.modeling.core.Model`. + + It could be a numpy array with the shape of the detector/observation. + Pixel values are of type integer or string and represent + region labels. Pixels which are not within any region have value ``no_label``. + + It could be a dictionary which maps tuples to labels or floating point numbers to labels. + + + :category:`Any of:` + + + + .. _http://stsci.edu/schemas/gwcs/label_mapper-1.0.0/allOf/1/properties/mapper/anyOf/0: + + :entry:`—` + + :soft:`Type:` `ndarray-1.0.0 `. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/label_mapper-1.0.0/allOf/1/properties/mapper/anyOf/1: + + :entry:`—` + + :soft:`Type:` `transform-1.1.0 `. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/label_mapper-1.0.0/allOf/1/properties/mapper/anyOf/2: + + :entry:`—` + + :soft:`Type:` object. + + + + + + :category:`Properties:` + + + + .. _http://stsci.edu/schemas/gwcs/label_mapper-1.0.0/allOf/1/properties/mapper/anyOf/2/properties/labels: + + :entry:`labels` + + :soft:`Type:` array :soft:`of` ( number :soft:`or` array :soft:`of` ( number ) ). + + + + + + :category:`Items:` + + + + .. _http://stsci.edu/schemas/gwcs/label_mapper-1.0.0/allOf/1/properties/mapper/anyOf/2/properties/labels/items: + + :soft:`Type:` number :soft:`or` array :soft:`of` ( number ). + + + + + + :category:`Any of:` + + + + .. _http://stsci.edu/schemas/gwcs/label_mapper-1.0.0/allOf/1/properties/mapper/anyOf/2/properties/labels/items/anyOf/0: + + :entry:`—` + + :soft:`Type:` number. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/label_mapper-1.0.0/allOf/1/properties/mapper/anyOf/2/properties/labels/items/anyOf/1: + + :entry:`—` + + :soft:`Type:` array :soft:`of` ( number ). + + + + + + :category:`Items:` + + + + .. _http://stsci.edu/schemas/gwcs/label_mapper-1.0.0/allOf/1/properties/mapper/anyOf/2/properties/labels/items/anyOf/1/items: + + :soft:`Type:` number. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/label_mapper-1.0.0/allOf/1/properties/mapper/anyOf/2/properties/models: + + :entry:`models` + + :soft:`Type:` array :soft:`of` ( `transform-1.1.0 ` ). + + + + + + :category:`Items:` + + + + .. _http://stsci.edu/schemas/gwcs/label_mapper-1.0.0/allOf/1/properties/mapper/anyOf/2/properties/models/items: + + :soft:`Type:` `transform-1.1.0 `. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/label_mapper-1.0.0/allOf/1/properties/inputs: + + :entry:`inputs` + + :soft:`Type:` array :soft:`of` ( string ). + + + + Names of inputs. + + + :category:`Items:` + + + + .. _http://stsci.edu/schemas/gwcs/label_mapper-1.0.0/allOf/1/properties/inputs/items: + + :soft:`Type:` string. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/label_mapper-1.0.0/allOf/1/properties/inputs_mapping: + + :entry:`inputs_mapping` + + :soft:`Type:` `transform-1.1.0 `. + + + + [mapping](ref:http://stsci.edu/schemas/asdf/transform/remap-axes-1.1.0) + + + + + .. _http://stsci.edu/schemas/gwcs/label_mapper-1.0.0/allOf/1/properties/atol: + + :entry:`atol` + + :soft:`Type:` number. + + + + absolute tolerance to compare keys in mapper. + + + + + .. _http://stsci.edu/schemas/gwcs/label_mapper-1.0.0/allOf/1/properties/no_label: + + :entry:`no_label` + + :soft:`Type:` number :soft:`or` string. + + + + Fill in value for missing output. + + + :category:`Any of:` + + + + .. _http://stsci.edu/schemas/gwcs/label_mapper-1.0.0/allOf/1/properties/no_label/anyOf/0: + + :entry:`—` + + :soft:`Type:` number. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/label_mapper-1.0.0/allOf/1/properties/no_label/anyOf/1: + + :entry:`—` + + :soft:`Type:` string. + + + + + +:category:`Examples:` + +Map array indices are to labels.:: + + ! + mapper: !core/ndarray-1.0.0 + data: + - [1, 0, 2] + - [1, 0, 2] + - [1, 0, 2] + datatype: int64 + shape: [3, 3] + no_label: 0 + + +Map numbers dictionary to transforms which return labels.:: + + ! + atol: 1.0e-08 + inputs: [x, y] + inputs_mapping: !transform/remap_axes-1.1.0 + mapping: [0] + n_inputs: 2 + mapper: !!omap + - !!omap + labels: [-1.67833272, -1.9580548, -1.118888] + - !!omap + models: + - !transform/shift-1.1.0 {offset: 6.0} + - !transform/shift-1.1.0 {offset: 2.0} + - !transform/shift-1.1.0 {offset: 4.0} + no_label: 0 + + +Map a number within a range of numbers to transforms which return labels.:: + + ! + mapper: !!omap + - !!omap + labels: + - [3.2, 4.1] + - [2.67, 2.98] + - [1.95, 2.3] + - !!omap + models: + - !transform/shift-1.1.0 {offset: 6.0} + - !transform/shift-1.1.0 {offset: 2.0} + - !transform/shift-1.1.0 {offset: 4.0} + inputs: [x, y] + inputs_mapping: !transform/remap_axes-1.1.0 + mapping: [0] + n_inputs: 2 + + diff --git a/docs/gwcs/schemas/regions_selector-1.0.0.rst b/docs/gwcs/schemas/regions_selector-1.0.0.rst new file mode 100644 index 00000000..f724c13e --- /dev/null +++ b/docs/gwcs/schemas/regions_selector-1.0.0.rst @@ -0,0 +1,234 @@ + + +.. _http://stsci.edu/schemas/gwcs/regions_selector-1.0.0: + +regions_selector-1.0.0: Represents a discontinuous transform. +============================================================= + +:soft:`Type:` `transform-1.1.0 ` :soft:`and` object. + +Represents a discontinuous transform. + + +Maps regions to transgorms and evaluates the transforms with the corresponding inputs. + + +:category:`All of:` + + + + .. _http://stsci.edu/schemas/gwcs/regions_selector-1.0.0/allOf/0: + + :entry:`0` + + :soft:`Type:` `transform-1.1.0 `. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/regions_selector-1.0.0/allOf/1: + + :entry:`1` + + :soft:`Type:` object. + + + + + + :category:`Properties:` + + + + .. _http://stsci.edu/schemas/gwcs/regions_selector-1.0.0/allOf/1/properties/label_mapper: + + :entry:`label_mapper` + + :soft:`Type:` :doc:`label_mapper-1.0.0 <./label_mapper-1.0.0>`. Required. + + + + An instance of + [label_mapper-1.1.0](ref:http://stsci.edu/schemas/gwcs/label_mapper-1.0.0) + + + + + .. _http://stsci.edu/schemas/gwcs/regions_selector-1.0.0/allOf/1/properties/inputs: + + :entry:`inputs` + + :soft:`Type:` array :soft:`of` ( string ). Required. + + + + Names of inputs. + + + :category:`Items:` + + + + .. _http://stsci.edu/schemas/gwcs/regions_selector-1.0.0/allOf/1/properties/inputs/items: + + :soft:`Type:` string. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/regions_selector-1.0.0/allOf/1/properties/outputs: + + :entry:`outputs` + + :soft:`Type:` array :soft:`of` ( string ). Required. + + + + Names of outputs. + + + :category:`Items:` + + + + .. _http://stsci.edu/schemas/gwcs/regions_selector-1.0.0/allOf/1/properties/outputs/items: + + :soft:`Type:` string. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/regions_selector-1.0.0/allOf/1/properties/selector: + + :entry:`selector` + + :soft:`Type:` object. Required. + + + + A mapping of regions to trransforms. + + + :category:`Properties:` + + + + .. _http://stsci.edu/schemas/gwcs/regions_selector-1.0.0/allOf/1/properties/selector/properties/labels: + + :entry:`labels` + + :soft:`Type:` array :soft:`of` ( integer or string ). + + + + An array of unique region labels. + + + :category:`Items:` + + + + .. _http://stsci.edu/schemas/gwcs/regions_selector-1.0.0/allOf/1/properties/selector/properties/labels/items: + + :soft:`Type:` integer or string. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/regions_selector-1.0.0/allOf/1/properties/selector/properties/transforms: + + :entry:`transforms` + + :soft:`Type:` array :soft:`of` ( `transform-1.1.0 ` ). + + + + A transform for each region. The order should match the order of labels. + + + :category:`Items:` + + + + .. _http://stsci.edu/schemas/gwcs/regions_selector-1.0.0/allOf/1/properties/selector/properties/transforms/items: + + :soft:`Type:` `transform-1.1.0 `. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/regions_selector-1.0.0/allOf/1/properties/undefined_transform_value: + + :entry:`undefined_transform_value` + + :soft:`Type:` number. + + + + Value to be returned if there's no transform defined for the inputs. + + +:category:`Examples:` + +Create a regions_selector schema for 2 regions, labeled "1" and "2".:: + + ! + inputs: [x, y] + label_mapper: ! + mapper: !core/ndarray-1.0.0 + datatype: int8 + data: + - [0, 1, 1, 0, 2, 0] + - [0, 1, 1, 0, 2, 0] + - [0, 1, 1, 0, 2, 0] + - [0, 1, 1, 0, 2, 0] + - [0, 1, 1, 0, 2, 0] + datatype: int64 + shape: [5, 6] + no_label: 0 + outputs: [ra, dec, lam] + selector: !!omap + - !!omap + labels: [1, 2] + - !!omap + transforms: + - !transform/compose-1.1.0 + forward: + - !transform/remap_axes-1.1.0 + mapping: [0, 1, 1] + - !transform/concatenate-1.1.0 + forward: + - !transform/concatenate-1.1.0 + forward: + - !transform/shift-1.1.0 {offset: 1.0} + - !transform/shift-1.1.0 {offset: 2.0} + - !transform/shift-1.1.0 {offset: 3.0} + - !transform/compose-1.1.0 + forward: + - !transform/remap_axes-1.1.0 + mapping: [0, 1, 1] + - !transform/concatenate-1.1.0 + forward: + - !transform/concatenate-1.1.0 + forward: + - !transform/scale-1.1.0 {factor: 2.0} + - !transform/scale-1.1.0 {factor: 3.0} + - !transform/scale-1.1.0 {factor: 3.0} + undefined_transform_value: .nan + + diff --git a/docs/gwcs/schemas/spectral_frame-1.0.0.rst b/docs/gwcs/schemas/spectral_frame-1.0.0.rst new file mode 100644 index 00000000..d1adfb37 --- /dev/null +++ b/docs/gwcs/schemas/spectral_frame-1.0.0.rst @@ -0,0 +1,93 @@ + + +.. _http://stsci.edu/schemas/gwcs/spectral_frame-1.0.0: + +spectral_frame-1.0.0: Represents a spectral frame. +================================================== + +:soft:`Type:` object :soft:`and` :doc:`frame-1.0.0 `. + +Represents a spectral frame. + + + + +:category:`All of:` + + + + .. _http://stsci.edu/schemas/gwcs/spectral_frame-1.0.0/allOf/0: + + :entry:`0` + + :soft:`Type:` object. + + + + + + :category:`Properties:` + + + + .. _http://stsci.edu/schemas/gwcs/spectral_frame-1.0.0/allOf/0/properties/reference_position: + + :entry:`reference_position` + + :soft:`Type:` any :soft:`from` ["geocenter", "barycenter", "heliocenter"]. + + + + The position of the reference frame. + + + :soft:`Default:` "geocenter" + + + + .. _http://stsci.edu/schemas/gwcs/spectral_frame-1.0.0/allOf/0/properties/axes_names: + + :entry:`axes_names` + + :soft:`Type:` any. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/spectral_frame-1.0.0/allOf/0/properties/axes_order: + + :entry:`axes_order` + + :soft:`Type:` any. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/spectral_frame-1.0.0/allOf/0/properties/unit: + + :entry:`unit` + + :soft:`Type:` any. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/spectral_frame-1.0.0/allOf/1: + + :entry:`1` + + :soft:`Type:` :doc:`frame-1.0.0 `. + + + + + diff --git a/docs/gwcs/schemas/step-1.0.0.rst b/docs/gwcs/schemas/step-1.0.0.rst new file mode 100644 index 00000000..352991b9 --- /dev/null +++ b/docs/gwcs/schemas/step-1.0.0.rst @@ -0,0 +1,101 @@ + + +.. _http://stsci.edu/schemas/gwcs/step-1.0.0: + +step-1.0.0: Describes a single step of a WCS transform pipeline. +================================================================ + +:soft:`Type:` object. + +Describes a single step of a WCS transform pipeline. + + + + +:category:`Properties:` + + + + .. _http://stsci.edu/schemas/gwcs/step-1.0.0/properties/frame: + + :entry:`frame` + + :soft:`Type:` string :soft:`or` :doc:`frame-1.0.0 `. Required. + + + + The frame of the inputs to the transform. + + + :category:`Any of:` + + + + .. _http://stsci.edu/schemas/gwcs/step-1.0.0/properties/frame/anyOf/0: + + :entry:`—` + + :soft:`Type:` string. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/step-1.0.0/properties/frame/anyOf/1: + + :entry:`—` + + :soft:`Type:` :doc:`frame-1.0.0 `. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/step-1.0.0/properties/transform: + + :entry:`transform` + + :soft:`Type:` `transform-1.1.0 ` :soft:`or` null. + + + + The transform from this step to the next one. The + last step in a WCS should not have a transform, but + exists only to describe the frames and units of the + final output axes. + + + :soft:`Default:` null + + :category:`Any of:` + + + + .. _http://stsci.edu/schemas/gwcs/step-1.0.0/properties/transform/anyOf/0: + + :entry:`—` + + :soft:`Type:` `transform-1.1.0 `. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/step-1.0.0/properties/transform/anyOf/1: + + :entry:`—` + + :soft:`Type:` null. + + + + + +:category:`Examples:` + diff --git a/docs/gwcs/schemas/temporal_frame-1.0.0.rst b/docs/gwcs/schemas/temporal_frame-1.0.0.rst new file mode 100644 index 00000000..e5f551d6 --- /dev/null +++ b/docs/gwcs/schemas/temporal_frame-1.0.0.rst @@ -0,0 +1,91 @@ + + +.. _http://stsci.edu/schemas/gwcs/temporal_frame-1.0.0: + +temporal_frame-1.0.0: Represents a temporal frame. +================================================== + +:soft:`Type:` object :soft:`and` :doc:`frame-1.0.0 `. + +Represents a temporal frame. + + + + +:category:`All of:` + + + + .. _http://stsci.edu/schemas/gwcs/temporal_frame-1.0.0/allOf/0: + + :entry:`0` + + :soft:`Type:` object. + + + + + + :category:`Properties:` + + + + .. _http://stsci.edu/schemas/gwcs/temporal_frame-1.0.0/allOf/0/properties/reference_time: + + :entry:`reference_time` + + :soft:`Type:` `time-1.1.0 `. + + + + The reference time. + + + + + .. _http://stsci.edu/schemas/gwcs/temporal_frame-1.0.0/allOf/0/properties/axes_names: + + :entry:`axes_names` + + :soft:`Type:` any. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/temporal_frame-1.0.0/allOf/0/properties/axes_order: + + :entry:`axes_order` + + :soft:`Type:` any. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/temporal_frame-1.0.0/allOf/0/properties/unit: + + :entry:`unit` + + :soft:`Type:` any. + + + + + + + + .. _http://stsci.edu/schemas/gwcs/temporal_frame-1.0.0/allOf/1: + + :entry:`1` + + :soft:`Type:` :doc:`frame-1.0.0 `. + + + + + diff --git a/docs/gwcs/schemas/wcs-1.0.0.rst b/docs/gwcs/schemas/wcs-1.0.0.rst new file mode 100644 index 00000000..6b05de50 --- /dev/null +++ b/docs/gwcs/schemas/wcs-1.0.0.rst @@ -0,0 +1,60 @@ + + +.. _http://stsci.edu/schemas/gwcs/wcs-1.0.0: + +wcs-1.0.0: A system for describing generalized world coordinate transformations. +================================================================================ + +:soft:`Type:` object. + +A system for describing generalized world coordinate transformations. + + +ASDF WCS is a way of specifying transformations (usually from detector space to world coordinate space and back) by using the transformations in the `transform-schema` module. + + +:category:`Properties:` + + + + .. _http://stsci.edu/schemas/gwcs/wcs-1.0.0/properties/name: + + :entry:`name` + + :soft:`Type:` string. Required. + + + + A descriptive name for this WCS. + + + + + .. _http://stsci.edu/schemas/gwcs/wcs-1.0.0/properties/steps: + + :entry:`steps` + + :soft:`Type:` array :soft:`of` ( :doc:`step-1.0.0 ` ). Required. + + + + A list of steps in the forward transformation from detector to + world coordinates. + The inverse transformation is determined automatically by + reversing this list, and inverting each of the individual + transforms according to the rules described in + [inverse](ref:http://stsci.edu/schemas/asdf/transform/transform-1.1.0/properties/inverse). + + + :category:`Items:` + + + + .. _http://stsci.edu/schemas/gwcs/wcs-1.0.0/properties/steps/items: + + :soft:`Type:` :doc:`step-1.0.0 `. + + + + + diff --git a/docs/gwcs/selector_model.rst b/docs/gwcs/selector_model.rst deleted file mode 100644 index 1c829244..00000000 --- a/docs/gwcs/selector_model.rst +++ /dev/null @@ -1,107 +0,0 @@ -Discontinuous WCS - An IFU Example -================================== - -There are a couple of models in GWCS which help with managing discontinuous WCSs. -Given (x, y) pixel indices, `~gwcs.selector.LabelMapperArray` returns a label (int or str) -which uniquely identifies a location in a frame. `~gwcs.selector.RegionsSelector` -maps different transforms to different locations in the frame of `~gwcs.selector._LabelMapper`. - -An example use case is an IFU observation. The locations on the detector where slices are -projected are labeled with integer numbers. Given an (x, y) pixel, the `~gwcs.selector.LabelMapperArray` -returns the label of the IFU slice which this pixel belongs to. Assuming each slice has its own WCS -transformation, `~gwcs.selector.RegionsSelector` takes as input an instance of `~gwcs.selector.LabelMapperArray` -and maps labels to transforms. A step by step example of constructing the WCS for an IFU with 6 slits follows. - - >>> import numpy as np - >>> from astropy.modeling.models import Shift, Scale, Mapping - >>> from astropy import coordinates as coord - >>> from astropy import units as u - >>> from gwcs import wcs, selector - >>> from gwcs import coordinate_frames as cf - -The output frame is common for all slits and is a composite frame with two subframes, -`~gwcs.coordinate_frames.CelestialFrame` and `~gwcs.coordinate_frames.SpectralFrame`. - - >>> sky_frame = cf.CelestialFrame(name='icrs', reference_frame=coord.ICRS(), axes_order=(0, 2)) - >>> spec_frame = cf.SpectralFrame(name='wave', unit=(u.micron,), axes_order=(1,), axes_names=('lambda',)) - >>> cframe = cf.CompositeFrame([sky_frame, spec_frame], name='world') - -The input frame can be a string (default is 'detector') or a subclass of `~gwcs.coordinate_frames.CoordinateFrame` -in which allows additional functionality (see example below). - - >>> det = cf.Frame2D(name='detector') - -All slits have the same input and output frames, however each slit has a different model transforming -from pixels to world coordinates (RA, lambda, dec). For the sake of brevity this example uses a simple -shift transform for each slit. Detailed examples of how to create more realistic transforms -are here (imaging) and here the ref: spectral_example. - - >>> transforms = {} - >>> for i in range(1, 7): - ... transforms[i] = Mapping([0, 0, 1]) | Shift(i * 0.1) & Shift(i * 0.2) & Scale(i * 0.1) - -One way to initialize `~gwcs.selector.LabelMapperArray` is to pass it the shape of the array and the vertices -of each slit on the detector {label: vertices} see :meth: `~gwcs.selector.LabelMapperArray.from_vertices`. -In this example the mask is an array with the size of the detector where each item in the array -corresponds to a pixel on the detector and its value is the slice number (label) this pixel -belongs to. - -.. image:: ifu-regions.png - -The image above shows the projection of the 6 slits on the detector. Pixels, with a label of 0 do -not belong to any slit. Assuming the array is stored in -`ASDF `__ format, create the mask: - -.. doctest-skip-all - - >>> from asdf import AsdfFile - >>> f = AsdfFile.open('mask.asdf') - >>> data = f.tree['mask'] - >>> mask = selector.LabelMapperArray(data) - -For more information on using the `ASDF standard `__ format -see `asdf `__ - -Create the pixel to world transform for the entire IFU: - - >>> regions_transform = selector.RegionsSelector(inputs=['x','y'], - ... outputs=['ra', 'dec', 'lam'], - ... selector=transforms, - ... label_mapper=mask, - ... undefined_transform_value=np.nan) - -The WCS object now can evaluate simultaneously the transforms of all slices - - >>> wifu = wcs.WCS(forward_transform=regions_transform, output_frame=cframe, input_frame=det) - >>> x, y = mask.mapper.shape - >>> x, y = np.mgrid[:x, :y] - >>> r, d, l = wifu(x, y) - -or of single slices. - - >>> wifu.forward_transform.set_input(4)(1, 2) - (1.4, 1.8, 0.8) - -The :meth:`~gwcs.selector.RegionsSelector.set_input` method returns the forward_transform for -a specific label. - -The above commands return numerical values. The :meth: `~gwcs.coordinate_frames.CoordinateFrame.coordinates` -is functionally equivalent to the above commands bu returns coordinate objects: - - >>> wifu(10, 200) - (10.3, 10.6, 60.00000000000001) - >>> wifu.available_frames - ['detector', 'world'] - - >>> wifu.output_frame.coordinates(10, 200) - (, ) - -Frames provide additional information: - - >>> print(wifu.output_frame.axes_type) - [u'SPATIAL', u'SPECTRAL', u'SPATIAL'] - >>> print(wifu.output_frame.axes_names) - [u'ra', 'lambda', u'dec'] - >>> print(wifu.output_frame.unit) - [Unit("deg"), Unit("micron"), Unit("deg")] diff --git a/docs/gwcs/using_wcs.rst b/docs/gwcs/using_wcs.rst index 2a8efe30..d6ac1b2d 100644 --- a/docs/gwcs/using_wcs.rst +++ b/docs/gwcs/using_wcs.rst @@ -1,95 +1,95 @@ +.. _user_api: + Using the WCS object ==================== - -Let's expand the WCS created in :ref:`getting-started` by adding a distortion. - -First create polynomial transform to represent distortion: - - >>> import numpy as np - >>> from astropy.modeling.models import (Polynomial2D, Shift, Scale, Rotation2D, - ... Pix2Sky_TAN, RotateNative2Celestial, Mapping) - >>> polyx = Polynomial2D(4) - >>> polyx.parameters = np.arange(15) * .1 - >>> polyy = Polynomial2D(4) - >>> polyy.parameters = np.arange(15) * .2 - >>> distortion = (Mapping((0, 1, 0, 1)) | polyx & polyy).rename("distortion") - >>> det2sky = (Shift(-10.5) & Shift(-13.2) | Rotation2D(0.0023) | \ - ... Scale(.01) & Scale(.04) | Pix2Sky_TAN() | \ - ... RotateNative2Celestial(5.6, -72.05, 180)).rename("det2sky") - -Create an intermediate frame. The distortion transforms positions on the -detector into this frame. - - >>> from astropy import units as u - >>> from astropy import coordinates as coord - >>> from gwcs import coordinate_frames as cf - >>> focal_frame = cf.Frame2D(name="focal_frame", unit=(u.arcsec, u.arcsec)) - >>> detector_frame = cf.Frame2D(name="detector", unit=(u.pix, u.pix)) - >>> sky_frame = cf.CelestialFrame(reference_frame=coord.ICRS(), name='icrs') - -Create the WCS pipeline and initialize the WCS: - - >>> from gwcs import wcs - >>> pipeline = [(detector_frame, distortion), - ... (focal_frame, det2sky), - ... (sky_frame, None) - ... ] - >>> wcsobj = wcs.WCS(pipeline) - >>> print(wcsobj) - From Transform - ----------- ---------- - detector distortion - focal_frame det2sky - icrs None + +This section uses the ``imaging_wcs.asdf`` created in :ref:`imaging_example` +to read in a WCS object and demo its methods. + +.. doctest-skip:: + + >>> import asdf + >>> asdf_file = asdf.open("imaging_wcs.asdf") + >>> wcsobj = asdf_file.tree["wcs"] + >>> print(wcsobj) # doctest: +SKIP + From Transform + ----------------- ---------------- + detector distortion + undistorted_frame linear_transform + icrs None To see what frames are defined: +.. doctest-skip:: + >>> print(wcsobj.available_frames) - ['detector', 'focal_frame', 'icrs'] + ['detector', 'undistorted_frame', 'icrs'] >>> wcsobj.input_frame - - >>> wcsobj.output_frame - )> + + >>> wcsobj.output_frame + )> Because the ``output_frame`` is a `~gwcs.coordinate_frames.CoordinateFrame` object we can get -the result of the WCS transform as an `astropy.coordinates.SkyCoord` object and transform +the result of the WCS transform as an `~astropy.coordinates.SkyCoord` object and transform them to other standard coordinate frames supported by `astropy.coordinates`. +.. doctest-skip:: + >>> skycoord = wcsobj(1, 2, with_units=True) - >>> print(skycoord) # doctest: +SKIP + >>> print(skycoord) - >>> print(skycoord.transform_to("galactic")) # doctest: +SKIP + (5.52886119, -72.05285219)> + >>> print(skycoord.transform_to("galactic")) + (306.11346489, -44.89382103)> +The WCS object has an attribute :attr:`~gwcs.WCS.bounding_box` +(default value of ``None``) which describes the range of +acceptable values for each input axis. + +.. doctest-skip:: + + >>> wcsobj.bounding_box = ((0, 2048), (0, 1000)) + >>> wcsobj((2,3), (1020, 980)) + array([nan, 133.48248429]), array([nan, -11.24021056]) + +The WCS object accepts a boolean flag called ``with_bounding_box`` with default value of +``True``. Output values which are outside the ``bounding_box`` are set to ``NaN``. +There are cases when this is not desirable and ``with_bounding_box=False`` should be passes. + +Calling the :meth:`~gwcs.WCS.footprint` returns the footprint on the sky. + +.. doctest-skip:: + + >>> wcsobj.footprint() + Some methods allow managing the transforms in a more detailed manner. Transforms between frames can be retrieved and evaluated separately. - >>> distortion = wcsobj.get_transform('detector', 'focal_frame') - >>> distortion(1, 2) # doctest: +FLOAT_CMP +.. doctest-skip:: + + + >>> dist = wcsobj.get_transform('detector', 'undistorted_frame') + >>> dist(1, 2) # doctest: +FLOAT_CMP (47.8, 95.60) Transforms in the pipeline can be replaced by new transforms. - >>> from astropy.modeling.models import Shift - >>> new_transform = Shift(1) & Shift(1.5) | distortion +.. doctest-skip:: + + >>> new_transform = models.Shift(1) & models.Shift(1.5) | distortion >>> wcsobj.set_transform('detector', 'focal_frame', new_transform) - >>> wcsobj(1, 2) - (10.338562883899195, -42.33182878519405) + >>> wcsobj(1, 2) # doctest: +FLOAT_CMP + (5.5583005430002785, -72.06028278184611) + A transform can be inserted before or after a frame in the pipeline. - >>> from astropy.modeling.models import Scale - >>> scale = Scale(2) & Scale(1) +.. doctest-skip:: + + >>> scale = models.Scale(2) & models.Scale(1) >>> wcsobj.insert_transform('icrs', scale, after=False) - >>> wcsobj(1, 2) - (20.67712576779839, -42.33182878519405) + >>> wcsobj(1, 2) # doctest: +FLOAT_CMP + (11.116601086000557, -72.06028278184611) -The WCS object has an attribute ``domain`` which describes the range of -acceptable values for each input axis. - - >>> wcsobj.bounding_box = ((0, 2048), (0, 1000)) diff --git a/docs/gwcs/wcs_ape.rst b/docs/gwcs/wcs_ape.rst new file mode 100644 index 00000000..f878fb92 --- /dev/null +++ b/docs/gwcs/wcs_ape.rst @@ -0,0 +1,14 @@ +.. _ape14: + +Common Interface for World Coordinate System - APE 14 +===================================================== + +To improve interoperability between packages, the Astropy Project and other interested parties have collaboratively defined a standardized application +programming interface (API) for world coordinate system objects to be used +in Python. This API is described in the Astropy Proposal for Enhancements (APE) 14: +`A shared Python interface for World Coordinate Systems +`_. + +The base classes that define the low- (`~astropy.wcs.wcsapi.BaseLowLevelWCS`) and high- (:class:`~astropy.wcs.wcsapi.BaseHighLevelWCS`) level APIs are in astropy. +GWCS implements both APIs. Once a gWCS object is created the API methods will be available. It is recommended that applications use the ``Common API`` to +ensure transparent use of ``GWCS`` and ``FITSWCS`` objects. diff --git a/docs/gwcs/wcs_validation.rst b/docs/gwcs/wcs_validation.rst new file mode 100644 index 00000000..b02cb558 --- /dev/null +++ b/docs/gwcs/wcs_validation.rst @@ -0,0 +1,20 @@ +.. _wcs_validation: + +WCS validation +============== + + +The WCS is validated when an object is read in or written to a file. +However, this happens transparently to the end user and knowing +the details of the validation machinery is not necessary to use or +construct a WCS object. + +GWCS uses the `Advanced Scientific Data Format `__ (ASDF) +to validate the transforms, coordinate frames and the overall WCS object structure. +ASDF makes use of abstract data type +definitions called ``schemas``. The serialization and deserialization happens in classes, +referred to as ``tags``. Most of the transform schemas live in the ``asdf-standard`` package while most of the transform tags live in ``astropy``. :ref:`gwcs-schemas` are available for the WCS object, coordinate frames and some WCS specific transforms. + +Packages using GWCS may create their own transforms and schemas and register them as an ``Asdf Extension``. If those are of general use, it is recommended they be included in astropy. + + diff --git a/docs/gwcs/wcstools.rst b/docs/gwcs/wcstools.rst index bc69efd6..d90f2d65 100644 --- a/docs/gwcs/wcstools.rst +++ b/docs/gwcs/wcstools.rst @@ -1,8 +1,14 @@ WCS User Tools ============== +`~gwcs.wcstools.grid_from_bounding_box` is a function which returns a grid of input points based on the bounding_box of the WCS. -The `~gwcs.wcstools` module contains functions of general usability related to WCS. + >>> from gwcs.wcstools import grid_from_bounding_box + >>> bounding_box = ((0, 4096), (0, 2048)) + >>> x, y = grid_from_bounding_box(bounding_box) + >>> ra, dec = w(x, y) # doctest: +SKIP + +The `~gwcs.wcstools` module contains functions of general usability. `~gwcs.wcstools.wcs_from_fiducial` is a function which given a fiducial in some coordinate system, returns a WCS object. @@ -23,9 +29,3 @@ Any additional transforms are prepended to the projection and sky rotation. >>> w(2048, 1024) # doctest: +FLOAT_CMP (5.46, -72.2) -`~gwcs.wcstools.grid_from_bounding_box` is a function which returns a grid of input points based on the bounding_box of the WCS. - - >>> from gwcs.wcstools import grid_from_bounding_box - >>> bounding_box = ((0, 4096), (0, 2048)) - >>> x, y = grid_from_bounding_box(bounding_box) - >>> ra, dec = w(x, y) \ No newline at end of file diff --git a/docs/index.rst b/docs/index.rst index eba3efd4..d31eb62b 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -4,36 +4,87 @@ GWCS Documentation `GWCS `__ is a package for managing the World Coordinate System (WCS) of astronomical data. -Introduction ------------- - - -`GWCS `__ takes a general approach to -the problem of expressing transformations between pixel and world coordinates. -The term "World Coordinate System" encapsulates this transformation, although -that term is sometimes misunderstood to mean "the FITS implementation of WCS" -(See `the discussion in APE14 `_ -for more on this topic). - -GWCS supports a data model which includes the entire transformation pipeline from -input coordinates (detector by default) to world coordinates. -Transformations can be chained, joined or combined with arithmetic operators -using the flexible framework of compound models in `~astropy.modeling`. -In the case of a celestial output frame `~astropy.coordinates` provides -further transformations between standard coordinate frames. -Spectral output coordinates are instances of `~astropy.units.Quantity` and are -transformed to other units with the tools in that package. The goal is to provide -a flexible toolkit which is easily extendable by adding new transforms and frames. - +Introduction & Motivation for GWCS +---------------------------------- + +The mapping from ‘pixel’ coordinates to corresponding ‘real-world’ coordinates (e.g. celestial coordinates, +spectral wavelength) is crucial to relating astronomical data to the phenomena they describe. Images and +other types of data often come encoded with information that describes this mapping – this is referred +to as the ‘World Coordinate System’ or WCS. The term WCS is often used to refer specifically +to the most widely used 'FITS implementation of WCS', but here unless specified WCS refers to +the broader concept of relating pixel ⟷ world. (See the discussion in `APE14 `__ +for more on this topic). + +The FITS WCS standard, currently the most widely used method of encoding WCS in data, describes a +set of required FITS header keywords and allowed values that describe how pixel ⟷ world transformations +should be done. This current paradigm of encoding data with only instructions on how to relate pixel to world, separate +from the transformation machinery itself, has several limitations: + +* Limited flexibility. WCS keywords and their values are rigidly defined so that the instructions are unambiguous. + This places limitations on, for example, describing geometric distortion in images since only a handful of distortion models are defined + in the FITS standard (and therefore can be encoded in FITS headers as WCS information). +* Separation of data from transformation pipelines. The machinery that transforms pixel ⟷ world + does not exist along side the data – there is merely a roadmap for how one *would* do the transformation. + External packages and libraries (e.g wcslib, or its Python interface astropy.wcs) must be + written to interpret the instructions and execute the transformation. These libraries + don’t allow easy access to coordinate frames along the course of the full pixel to world + transformation pipeline. Additionally, since these libraries can only interpret FITS WCS + information, any custom ‘WCS’ definitions outside of FITS require the user to write their own transformation pipelines. +* Incompatibility with varying file formats. New file formats that are becoming more widely + used in place of FITS to store astronomical data, like the ASDF format, also require a + method of encoding WCS information. FITS WCS and the accompanying libraries are adapted for + FITS only. A more flexible interface would be agnostic to file type, as long as the necessary + information is present. + +The `GWCS `__ package and GWCS object is a generalized WCS +implementation that mitigates these limitations. The goal of the GWCS package is to provide a +flexible toolkit for expressing and evaluating transformations between pixel and world coordinates, +as well as intermediate frames along the course of this transformation.The GWCS object supports a +data model which includes the entire transformation pipeline from input pixel coordinates to +world coordinates (and vice versa). The basis of the GWCS object is astropy `modeling `__. +Models that describe the pixel ⟷ world transformations can be chained, joined or combined with arithmetic operators +using the flexible framework of compound models in modeling. This approach allows for easy +access to intermediate frames. In the case of a celestial output frame `coordinates `__. provides further transformations between +standard celestial coordinate frames. Spectral output coordinates are instances of Quantity +and can be transformed to other units with the tools in that package. ``Time`` coordinates are instances of ``Time``. +GWCS supports transforms initialized with `Quantity `__ objects ensuring automatic +unit conversion. + +Pixel Conventions and Definitions +--------------------------------- + +This API assumes that integer pixel values fall at the center of pixels (as +assumed in the FITS-WCS standard, see Section 2.1.4 of `Greisen et al., 2002, +A&A 446, 747 `_), while at the same +time matching the Python 0-index philosophy. That is, the first pixel is +considered pixel ``0``, but pixel coordinates ``(0, 0)`` are the *center* of +that pixel. Hence the first pixel spans pixel values ``-0.5`` to ``0.5``. + +There are two main conventions for ordering pixel coordinates. In the context of +2-dimensional imaging data/arrays, one can either think of the pixel coordinates +as traditional Cartesian coordinates (which we call ``x`` and ``y`` here), which +are usually given with the horizontal coordinate (``x``) first, and the vertical +coordinate (``y``) second, meaning that pixel coordinates would be given as +``(x, y)``. Alternatively, one can give the coordinates by first giving the row +in the data, then the column, i.e. ``(row, column)``. While the former is a more +common convention when e.g. plotting (think for example of the Matplotlib +``scatter(x, y)`` method), the latter is the convention used when accessing +values from e.g. Numpy arrays that represent images (``image[row, column]``). + +The GWCS object assumes Cartesian order ``(x, y)``, however the :ref:`ape14` accepts both conventions. +The order of the pixel coordinates (``(x, y)`` vs ``(row, column)``) in the ``Common API`` depends on the method or property used, and this can normally be +determined from the property or method name. Properties and methods containing +``pixel`` assume ``(x, y)`` ordering, while properties and methods containing +``array`` assume ``(row, column)`` ordering. Installation ------------ `gwcs `__ requires: -- `numpy `__ 1.7 or later +- `numpy `__ -- `astropy `__ 1.2 or later +- `astropy `__ - `asdf `__ @@ -47,93 +98,228 @@ To install the latest release:: pip install gwcs -GWCS is also available as part of `astroconda `__. +The latest release of GWCS is also available as part of `astroconda `__. .. _getting-started: Getting Started -=============== +--------------- -The WCS data model represents a pipeline of transformations from some -initial coordinate frame to a standard coordinate frame. -It is implemented as a list of steps executed in order. Each step defines a +The WCS data model represents a pipeline of transformations between two +coordinate frames, the final one usually a physical coordinate system. +It is represented as a list of steps executed in order. Each step defines a starting coordinate frame and the transform to the next frame in the pipeline. The last step has no transform, only a frame which is the output frame of the total transform. As a minimum a WCS object has an ``input_frame`` (defaults to "detector"), an ``output_frame`` and the transform between them. -As an example, consider a typical WCS of an image without distortion. +The WCS is validated using the `ASDF Standard `__ +and serialized to file using the `asdf `__ package. +There are two ways to save the WCS to a file: + +- `Save a WCS object as a pure ASDF file`_ + +- `Save a WCS object as an ASDF extension in a FITS file`_ + + +A step by step example of constructing an imaging GWCS object. +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The following example shows how to construct a GWCS object equivalent to +a FITS imaging WCS without distortion, defined in this FITS imaging header:: + + WCSAXES = 2 / Number of coordinate axes + WCSNAME = '47 Tuc ' / Coordinate system title + CRPIX1 = 2048.0 / Pixel coordinate of reference point + CRPIX2 = 1024.0 / Pixel coordinate of reference point + PC1_1 = 1.290551569736E-05 / Coordinate transformation matrix element + PC1_2 = 5.9525007864732E-06 / Coordinate transformation matrix element + PC2_1 = 5.0226382102765E-06 / Coordinate transformation matrix element + PC2_2 = -1.2644844123757E-05 / Coordinate transformation matrix element + CDELT1 = 1.0 / [deg] Coordinate increment at reference point + CDELT2 = 1.0 / [deg] Coordinate increment at reference point + CUNIT1 = 'deg' / Units of coordinate increment and value + CUNIT2 = 'deg' / Units of coordinate increment and value + CTYPE1 = 'RA---TAN' / TAN (gnomonic) projection + SIP distortions + CTYPE2 = 'DEC--TAN' / TAN (gnomonic) projection + SIP distortions + CRVAL1 = 5.63056810618 / [deg] Coordinate value at reference point + CRVAL2 = -72.05457184279 / [deg] Coordinate value at reference point + LONPOLE = 180.0 / [deg] Native longitude of celestial pole + LATPOLE = -72.05457184279 / [deg] Native latitude of celestial pole + RADESYS = 'ICRS' / Equatorial coordinate system + The following imports are generally useful: >>> import numpy as np - >>> from astropy.modeling.models import (Shift, Scale, Rotation2D, - ... Pix2Sky_TAN, RotateNative2Celestial, Mapping, Polynomial2D) + >>> from astropy.modeling import models >>> from astropy import coordinates as coord >>> from astropy import units as u >>> from gwcs import wcs >>> from gwcs import coordinate_frames as cf The ``forward_transform`` is constructed as a combined model using `astropy.modeling`. -The frames are subclasses of `~gwcs.coordinate_frames.CoordinateFrame` (although strings are -acceptable too). +The ``frames`` are subclasses of `~gwcs.coordinate_frames.CoordinateFrame`. Although strings are +acceptable as ``coordinate_frames`` it is recommended this is used only in testing/debugging. + +Using the `~astropy.modeling` package create a combined model to transform +detector coordinates to ICRS following the FITS WCS standard convention. + +First, create a transform which shifts the input ``x`` and ``y`` coordinates by ``CRPIX``: + + >>> shift_by_crpix = models.Shift(-2048*u.pix) & models.Shift(-1024*u.pix) -Create a transform to convert detector coordinates to ICRS. +Create a transform which rotates the inputs using the ``PC matrix``. - >>> det2sky = (Shift(-10.5) & Shift(-13.2) | Rotation2D(0.0023) | \ - ... Scale(.01) & Scale(.04) | Pix2Sky_TAN() | \ - ... RotateNative2Celestial(5.6, -72.05, 180)).rename("detector2sky") + >>> 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)} -Create a coordinate frame associated with the detector and a celestial frame. +Create a tangent projection and a rotation on the sky using ``CRVAL``. - >>> 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') + >>> 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" -Initialize the WCS: +Create a ``detector`` coordinate frame and a ``celestial`` ICRS frame. - >>> wcsobj = wcs.WCS(forward_transform=det2sky, input_frame=detector_frame, output_frame=sky_frame) + >>> 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)) + +This WCS pipeline has only one step - from ``detector`` to ``sky``: + + >>> pipeline = [(detector_frame, det2sky), + ... (sky_frame, None) + ... ] + >>> wcsobj = wcs.WCS(pipeline) >>> print(wcsobj) - From Transform - -------- ------------ - detector detector2sky - icrs None + From Transform + -------- ---------------- + detector linear_transform + icrs None To convert a pixel (x, y) = (1, 2) to sky coordinates, call the WCS object as a function: - >>> sky = wcsobj(1, 2) + >>> sky = wcsobj(1*u.pix, 2*u.pix, with_units=True) >>> print(sky) - (5.284139265842838, -72.49775640633504) + 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. - >>> wcsobj.invert(*sky) # doctest: +FLOAT_CMP - (1.0, 2.0) + >>> wcsobj.invert(sky) + (, ) + +.. _save_as_asdf: + +Save a WCS object as a pure ASDF file ++++++++++++++++++++++++++++++++++++++ + +.. doctest-skip:: + + >>> from asdf import AsdfFile + >>> tree = {"wcs": wcsobj} + >>> wcs_file = AsdfFile(tree) + >>> wcs_file.write_to("imaging_wcs.asdf") + + +:ref:`pure_asdf` + + +Save a WCS object as an ASDF extension in a FITS file ++++++++++++++++++++++++++++++++++++++++++++++++++++++ +.. doctest-skip:: + + >>> from astropy.io import fits + >>> from asdf import fits_embed + >>> hdul = fits.open("example_imaging.fits") + >>> hdul.info() + Filename: example_imaging.fits + No. Name Ver Type Cards Dimensions Format + 0 PRIMARY 1 PrimaryHDU 775 () + 1 SCI 1 ImageHDU 71 (600, 550) float32 + >>> tree = {"sci", hdul.data, + ... "wcs": wcsobj} + >>> fa = fits.embed.AsdfInFits(hdul, tree) + >>> fa.write_to("imaging_with_wcs_in_asdf.fits") + >>> fits.info("imaging_with_wcs_in_asdf.fits") + Filename: example_with_wcs.asdf + No. Name Ver Type Cards Dimensions Format + 0 PRIMARY 1 PrimaryHDU 775 () + 1 SCI 1 ImageHDU 71 (600, 550) float32 + 2 ASDF 1 BinTableHDU 11 1R x 1C [5086B] + +Reading a WCS object from a file +++++++++++++++++++++++++++++++++ + + +`ASDF `__ is used to read a WCS object from a + pure ASDF file or from an ASDF extension in a FITS file. + + +.. doctest-skip:: + + >>> import asdf + >>> asdf_file = asdf.open("imaging_wcs.asdf") + >>> wcsobj = asdf_file.tree['wcs'] + + +.. doctest-skip:: + + >>> import asdf + >>> fa = asdf.open("imaging_with_wcs_in_asdf.fits") + >>> wcsobj = fa.tree["wcs"] + +Other Examples +-------------- + +.. toctree:: + :maxdepth: 2 + + gwcs/imaging_with_distortion.rst + gwcs/ifu.rst + + + Using `gwcs` ------------ - .. toctree:: :maxdepth: 2 + gwcs/wcs_ape.rst gwcs/using_wcs.rst - gwcs/selector_model.rst gwcs/wcstools.rst - + gwcs/pure_asdf.rst + gwcs/wcs_validation.rst + gwcs/schemas/index.rst + gwcs/points_to_wcs.rst See also -------- -- The `modeling package in astropy +- `The modeling package in astropy `__ -- The `Coordinates package in astropy +- `The coordinates package in astropy `__ -- The `Advanced Scientific Data Format (ASDF) standard +- `The Advanced Scientific Data Format (ASDF) standard `__ and its `Python implementation `__ @@ -144,5 +330,6 @@ Reference/API .. automodapi:: gwcs.wcs .. automodapi:: gwcs.coordinate_frames -.. automodapi:: gwcs.selector .. automodapi:: gwcs.wcstools +.. automodapi:: gwcs.selector + diff --git a/docs/rtd-pip-requirements b/docs/rtd-pip-requirements deleted file mode 100644 index 0a5f9620..00000000 --- a/docs/rtd-pip-requirements +++ /dev/null @@ -1,6 +0,0 @@ -numpy -Cython -scipy > 0.18 -astropy_helpers -astropy -asdf diff --git a/docs/sphinxext/__init__.py b/docs/sphinxext/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/docs/sphinxext/category.py b/docs/sphinxext/category.py new file mode 100644 index 00000000..1f532f98 --- /dev/null +++ b/docs/sphinxext/category.py @@ -0,0 +1,119 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +# -*- coding: utf-8 -*- + +#from __future__ import absolute_import, division, unicode_literals, print_function + + +from docutils import nodes + + +class category(nodes.Element): + pass + + +def category_role(name, rawtext, text, lineno, inliner, + options={}, content=[]): + node = category(text) + return [node], [] + + +def visit_category_node_html(self, node): + self.body.append('') + self.body.append(node.rawsource) + + +def depart_category_node_html(self, node): + self.body.append('') + + +def visit_category_node_latex(self, node): + self.body.append('\n\n\\textsf{%s}' % node.rawsource) + + +def depart_category_node_latex(self, node): + pass + + +class entry(nodes.Element): + pass + + +def entry_role(name, rawtext, text, lineno, inliner, + options={}, content=[]): + node = entry(text) + return [node], [] + + +def visit_entry_node_html(self, node): + self.body.append('') + self.body.append(node.rawsource) + + +def depart_entry_node_html(self, node): + self.body.append('') + + +def visit_entry_node_latex(self, node): + self.body.append('\n\n\\texttt{%s}' % + node.rawsource.replace('_', r'\_')) + + +def depart_entry_node_latex(self, node): + pass + + +class soft(nodes.Element): + pass + + +def soft_role(name, rawtext, text, lineno, inliner, + options={}, content=[]): + node = soft(text) + return [node], [] + + +def visit_soft_node_html(self, node): + self.body.append('') + self.body.append(node.rawsource) + + +def depart_soft_node_html(self, node): + self.body.append('') + + +def visit_soft_node_latex(self, node): + self.body.append('\\textcolor[gray]{0.5}{%s}' % node.rawsource) + + +def depart_soft_node_latex(self, node): + pass + + +def setup(app): + app.add_node( + category, + html=(visit_category_node_html, depart_category_node_html), + latex=(visit_category_node_latex, depart_category_node_latex) + ) + + app.add_role( + 'category', category_role) + + app.add_node( + entry, + html=(visit_entry_node_html, depart_entry_node_html), + #latex=(visit_entry_node_latex, depart_entry_node_latex), + override=True + ) + + app.add_role( + 'entry', entry_role) + + app.add_node( + soft, + html=(visit_soft_node_html, depart_soft_node_html), + #latex=(visit_soft_node_latex, depart_soft_node_latex) + ) + + app.add_role( + 'soft', soft_role) diff --git a/ez_setup.py b/ez_setup.py deleted file mode 100644 index 800c31ef..00000000 --- a/ez_setup.py +++ /dev/null @@ -1,414 +0,0 @@ -#!/usr/bin/env python - -""" -Setuptools bootstrapping installer. - -Maintained at https://github.com/pypa/setuptools/tree/bootstrap. - -Run this script to install or upgrade setuptools. - -This method is DEPRECATED. Check https://github.com/pypa/setuptools/issues/581 for more details. -""" - -import os -import shutil -import sys -import tempfile -import zipfile -import optparse -import subprocess -import platform -import textwrap -import contextlib - -from distutils import log - -try: - from urllib.request import urlopen -except ImportError: - from urllib2 import urlopen - -try: - from site import USER_SITE -except ImportError: - USER_SITE = None - -# 33.1.1 is the last version that supports setuptools self upgrade/installation. -DEFAULT_VERSION = "33.1.1" -DEFAULT_URL = "https://pypi.io/packages/source/s/setuptools/" -DEFAULT_SAVE_DIR = os.curdir -DEFAULT_DEPRECATION_MESSAGE = "ez_setup.py is deprecated and when using it setuptools will be pinned to {0} since it's the last version that supports setuptools self upgrade/installation, check https://github.com/pypa/setuptools/issues/581 for more info; use pip to install setuptools" - -MEANINGFUL_INVALID_ZIP_ERR_MSG = 'Maybe {0} is corrupted, delete it and try again.' - -log.warn(DEFAULT_DEPRECATION_MESSAGE.format(DEFAULT_VERSION)) - - -def _python_cmd(*args): - """ - Execute a command. - - Return True if the command succeeded. - """ - args = (sys.executable,) + args - return subprocess.call(args) == 0 - - -def _install(archive_filename, install_args=()): - """Install Setuptools.""" - with archive_context(archive_filename): - # installing - log.warn('Installing Setuptools') - if not _python_cmd('setup.py', 'install', *install_args): - log.warn('Something went wrong during the installation.') - log.warn('See the error message above.') - # exitcode will be 2 - return 2 - - -def _build_egg(egg, archive_filename, to_dir): - """Build Setuptools egg.""" - with archive_context(archive_filename): - # building an egg - log.warn('Building a Setuptools egg in %s', to_dir) - _python_cmd('setup.py', '-q', 'bdist_egg', '--dist-dir', to_dir) - # returning the result - log.warn(egg) - if not os.path.exists(egg): - raise IOError('Could not build the egg.') - - -class ContextualZipFile(zipfile.ZipFile): - - """Supplement ZipFile class to support context manager for Python 2.6.""" - - def __enter__(self): - return self - - def __exit__(self, type, value, traceback): - self.close() - - def __new__(cls, *args, **kwargs): - """Construct a ZipFile or ContextualZipFile as appropriate.""" - if hasattr(zipfile.ZipFile, '__exit__'): - return zipfile.ZipFile(*args, **kwargs) - return super(ContextualZipFile, cls).__new__(cls) - - -@contextlib.contextmanager -def archive_context(filename): - """ - Unzip filename to a temporary directory, set to the cwd. - - The unzipped target is cleaned up after. - """ - tmpdir = tempfile.mkdtemp() - log.warn('Extracting in %s', tmpdir) - old_wd = os.getcwd() - try: - os.chdir(tmpdir) - try: - with ContextualZipFile(filename) as archive: - archive.extractall() - except zipfile.BadZipfile as err: - if not err.args: - err.args = ('', ) - err.args = err.args + ( - MEANINGFUL_INVALID_ZIP_ERR_MSG.format(filename), - ) - raise - - # going in the directory - subdir = os.path.join(tmpdir, os.listdir(tmpdir)[0]) - os.chdir(subdir) - log.warn('Now working in %s', subdir) - yield - - finally: - os.chdir(old_wd) - shutil.rmtree(tmpdir) - - -def _do_download(version, download_base, to_dir, download_delay): - """Download Setuptools.""" - py_desig = 'py{sys.version_info[0]}.{sys.version_info[1]}'.format(sys=sys) - tp = 'setuptools-{version}-{py_desig}.egg' - egg = os.path.join(to_dir, tp.format(**locals())) - if not os.path.exists(egg): - archive = download_setuptools(version, download_base, - to_dir, download_delay) - _build_egg(egg, archive, to_dir) - sys.path.insert(0, egg) - - # Remove previously-imported pkg_resources if present (see - # https://bitbucket.org/pypa/setuptools/pull-request/7/ for details). - if 'pkg_resources' in sys.modules: - _unload_pkg_resources() - - import setuptools - setuptools.bootstrap_install_from = egg - - -def use_setuptools( - version=DEFAULT_VERSION, download_base=DEFAULT_URL, - to_dir=DEFAULT_SAVE_DIR, download_delay=15): - """ - Ensure that a setuptools version is installed. - - Return None. Raise SystemExit if the requested version - or later cannot be installed. - """ - to_dir = os.path.abspath(to_dir) - - # prior to importing, capture the module state for - # representative modules. - rep_modules = 'pkg_resources', 'setuptools' - imported = set(sys.modules).intersection(rep_modules) - - try: - import pkg_resources - pkg_resources.require("setuptools>=" + version) - # a suitable version is already installed - return - except ImportError: - # pkg_resources not available; setuptools is not installed; download - pass - except pkg_resources.DistributionNotFound: - # no version of setuptools was found; allow download - pass - except pkg_resources.VersionConflict as VC_err: - if imported: - _conflict_bail(VC_err, version) - - # otherwise, unload pkg_resources to allow the downloaded version to - # take precedence. - del pkg_resources - _unload_pkg_resources() - - return _do_download(version, download_base, to_dir, download_delay) - - -def _conflict_bail(VC_err, version): - """ - Setuptools was imported prior to invocation, so it is - unsafe to unload it. Bail out. - """ - conflict_tmpl = textwrap.dedent(""" - The required version of setuptools (>={version}) is not available, - and can't be installed while this script is running. Please - install a more recent version first, using - 'easy_install -U setuptools'. - - (Currently using {VC_err.args[0]!r}) - """) - msg = conflict_tmpl.format(**locals()) - sys.stderr.write(msg) - sys.exit(2) - - -def _unload_pkg_resources(): - sys.meta_path = [ - importer - for importer in sys.meta_path - if importer.__class__.__module__ != 'pkg_resources.extern' - ] - del_modules = [ - name for name in sys.modules - if name.startswith('pkg_resources') - ] - for mod_name in del_modules: - del sys.modules[mod_name] - - -def _clean_check(cmd, target): - """ - Run the command to download target. - - If the command fails, clean up before re-raising the error. - """ - try: - subprocess.check_call(cmd) - except subprocess.CalledProcessError: - if os.access(target, os.F_OK): - os.unlink(target) - raise - - -def download_file_powershell(url, target): - """ - Download the file at url to target using Powershell. - - Powershell will validate trust. - Raise an exception if the command cannot complete. - """ - target = os.path.abspath(target) - ps_cmd = ( - "[System.Net.WebRequest]::DefaultWebProxy.Credentials = " - "[System.Net.CredentialCache]::DefaultCredentials; " - '(new-object System.Net.WebClient).DownloadFile("%(url)s", "%(target)s")' - % locals() - ) - cmd = [ - 'powershell', - '-Command', - ps_cmd, - ] - _clean_check(cmd, target) - - -def has_powershell(): - """Determine if Powershell is available.""" - if platform.system() != 'Windows': - return False - cmd = ['powershell', '-Command', 'echo test'] - with open(os.path.devnull, 'wb') as devnull: - try: - subprocess.check_call(cmd, stdout=devnull, stderr=devnull) - except Exception: - return False - return True -download_file_powershell.viable = has_powershell - - -def download_file_curl(url, target): - cmd = ['curl', url, '--location', '--silent', '--output', target] - _clean_check(cmd, target) - - -def has_curl(): - cmd = ['curl', '--version'] - with open(os.path.devnull, 'wb') as devnull: - try: - subprocess.check_call(cmd, stdout=devnull, stderr=devnull) - except Exception: - return False - return True -download_file_curl.viable = has_curl - - -def download_file_wget(url, target): - cmd = ['wget', url, '--quiet', '--output-document', target] - _clean_check(cmd, target) - - -def has_wget(): - cmd = ['wget', '--version'] - with open(os.path.devnull, 'wb') as devnull: - try: - subprocess.check_call(cmd, stdout=devnull, stderr=devnull) - except Exception: - return False - return True -download_file_wget.viable = has_wget - - -def download_file_insecure(url, target): - """Use Python to download the file, without connection authentication.""" - src = urlopen(url) - try: - # Read all the data in one block. - data = src.read() - finally: - src.close() - - # Write all the data in one block to avoid creating a partial file. - with open(target, "wb") as dst: - dst.write(data) -download_file_insecure.viable = lambda: True - - -def get_best_downloader(): - downloaders = ( - download_file_powershell, - download_file_curl, - download_file_wget, - download_file_insecure, - ) - viable_downloaders = (dl for dl in downloaders if dl.viable()) - return next(viable_downloaders, None) - - -def download_setuptools( - version=DEFAULT_VERSION, download_base=DEFAULT_URL, - to_dir=DEFAULT_SAVE_DIR, delay=15, - downloader_factory=get_best_downloader): - """ - Download setuptools from a specified location and return its filename. - - `version` should be a valid setuptools version number that is available - as an sdist for download under the `download_base` URL (which should end - with a '/'). `to_dir` is the directory where the egg will be downloaded. - `delay` is the number of seconds to pause before an actual download - attempt. - - ``downloader_factory`` should be a function taking no arguments and - returning a function for downloading a URL to a target. - """ - # making sure we use the absolute path - to_dir = os.path.abspath(to_dir) - zip_name = "setuptools-%s.zip" % version - url = download_base + zip_name - saveto = os.path.join(to_dir, zip_name) - if not os.path.exists(saveto): # Avoid repeated downloads - log.warn("Downloading %s", url) - downloader = downloader_factory() - downloader(url, saveto) - return os.path.realpath(saveto) - - -def _build_install_args(options): - """ - Build the arguments to 'python setup.py install' on the setuptools package. - - Returns list of command line arguments. - """ - return ['--user'] if options.user_install else [] - - -def _parse_args(): - """Parse the command line for options.""" - parser = optparse.OptionParser() - parser.add_option( - '--user', dest='user_install', action='store_true', default=False, - help='install in user site package') - parser.add_option( - '--download-base', dest='download_base', metavar="URL", - default=DEFAULT_URL, - help='alternative URL from where to download the setuptools package') - parser.add_option( - '--insecure', dest='downloader_factory', action='store_const', - const=lambda: download_file_insecure, default=get_best_downloader, - help='Use internal, non-validating downloader' - ) - parser.add_option( - '--version', help="Specify which version to download", - default=DEFAULT_VERSION, - ) - parser.add_option( - '--to-dir', - help="Directory to save (and re-use) package", - default=DEFAULT_SAVE_DIR, - ) - options, args = parser.parse_args() - # positional arguments are ignored - return options - - -def _download_args(options): - """Return args for download_setuptools function from cmdline args.""" - return dict( - version=options.version, - download_base=options.download_base, - downloader_factory=options.downloader_factory, - to_dir=options.to_dir, - ) - - -def main(): - """Install or upgrade setuptools and EasyInstall.""" - options = _parse_args() - archive = download_setuptools(**_download_args(options)) - return _install(archive, _build_install_args(options)) - -if __name__ == '__main__': - sys.exit(main()) diff --git a/gwcs/__init__.py b/gwcs/__init__.py index 82866834..696c5516 100644 --- a/gwcs/__init__.py +++ b/gwcs/__init__.py @@ -13,13 +13,13 @@ # Affiliated packages may add whatever they like to this file, but # should keep this content at the top. # ---------------------------------------------------------------------------- -from ._astropy_init import * +from ._astropy_init import * # noqa # ---------------------------------------------------------------------------- # For egg_info test builds to pass, put package imports here. if not _ASTROPY_SETUP_: - from .wcs import * - from .wcstools import * - from .coordinate_frames import * - from .selector import * + from .wcs import * # noqa + from .wcstools import * # noqa + from .coordinate_frames import * # noqa + from .selector import * # noqa diff --git a/gwcs/api.py b/gwcs/api.py new file mode 100644 index 00000000..9eca32eb --- /dev/null +++ b/gwcs/api.py @@ -0,0 +1,212 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +""" +This module contains a mixin class which exposes the WCS API defined +in astropy APE 14 (https://doi.org/10.5281/zenodo.1188875). + +""" +from astropy.wcs.wcsapi import BaseHighLevelWCS, BaseLowLevelWCS +from astropy.modeling import separable + +from . import utils + +__all__ = ["GWCSAPIMixin"] + + +class GWCSAPIMixin(BaseHighLevelWCS, BaseLowLevelWCS): + """ + A mix-in class that is intended to be inherited by the + :class:`~gwcs.wcs.WCS` class and provides the low- and high-level + WCS API described in the astropy APE 14 + (https://doi.org/10.5281/zenodo.1188875). + """ + + # Low Level APE 14 API + @property + def pixel_n_dim(self): + """ + The number of axes in the pixel coordinate system. + """ + return self.input_frame.naxes + + @property + def world_n_dim(self): + """ + The number of axes in the world coordinate system. + """ + return self.output_frame.naxes + + @property + def world_axis_physical_types(self): + """ + An iterable of strings describing the physical type for each world axis. + These should be names from the VO UCD1+ controlled Vocabulary + (http://www.ivoa.net/documents/latest/UCDlist.html). If no matching UCD + type exists, this can instead be ``"custom:xxx"``, where ``xxx`` is an + arbitrary string. Alternatively, if the physical type is + unknown/undefined, an element can be `None`. + """ + return self.output_frame.axis_physical_types + + @property + def world_axis_units(self): + """ + An iterable of strings given the units of the world coordinates for each + axis. + The strings should follow the `IVOA VOUnit standard + `_ (though as noted in the VOUnit + specification document, units that do not follow this standard are still + allowed, but just not recommended). + """ + return tuple(unit.to_string(format='vounit') for unit in self.output_frame.unit) + + def pixel_to_world_values(self, *pixel_arrays): + """ + Convert pixel coordinates to world coordinates. + + This method takes ``pixel_n_dim`` scalars or arrays as input, and pixel + coordinates should be zero-based. Returns ``world_n_dim`` scalars or + arrays in units given by ``world_axis_units``. Note that pixel + coordinates are assumed to be 0 at the center of the first pixel in each + dimension. If a pixel is in a region where the WCS is not defined, NaN + can be returned. The coordinates should be specified in the ``(x, y)`` + order, where for an image, ``x`` is the horizontal coordinate and ``y`` + is the vertical coordinate. + """ + result = self(*pixel_arrays, with_units=False) + return result + + def array_index_to_world_values(self, *index_arrays): + """ + Convert array indices to world coordinates. + This is the same as `~BaseLowLevelWCS.pixel_to_world_values` except that + the indices should be given in ``(i, j)`` order, where for an image + ``i`` is the row and ``j`` is the column (i.e. the opposite order to + `~BaseLowLevelWCS.pixel_to_world_values`). + """ + result = self(*index_arrays[::-1], with_units=False) + return result + + def world_to_pixel_values(self, *world_arrays): + """ + Convert world coordinates to pixel coordinates. + + This method takes ``world_n_dim`` scalars or arrays as input in units + given by ``world_axis_units``. Returns ``pixel_n_dim`` scalars or + arrays. Note that pixel coordinates are assumed to be 0 at the center of + the first pixel in each dimension. If a world coordinate does not have a + matching pixel coordinate, NaN can be returned. The coordinates should + be returned in the ``(x, y)`` order, where for an image, ``x`` is the + horizontal coordinate and ``y`` is the vertical coordinate. + """ + result = self.invert(*world_arrays, with_units=False) + return result + + + def world_to_array_index_values(self, *world_arrays): + """ + Convert world coordinates to array indices. + This is the same as `~BaseLowLevelWCS.world_to_pixel_values` except that + the indices should be returned in ``(i, j)`` order, where for an image + ``i`` is the row and ``j`` is the column (i.e. the opposite order to + `~BaseLowLevelWCS.pixel_to_world_values`). The indices should be + returned as rounded integers. + """ + result = self.invert(*world_arrays, with_units=False)[::-1] + return result # astype(int) + + @property + def array_shape(self): + """ + The shape of the data that the WCS applies to as a tuple of + length `~BaseLowLevelWCS.pixel_n_dim`. + If the WCS is valid in the context of a dataset with a particular + shape, then this property can be used to store the shape of the + data. This can be used for example if implementing slicing of WCS + objects. This is an optional property, and it should return `None` + if a shape is not known or relevant. + The shape should be given in ``(row, column)`` order (the convention + for arrays in Python). + """ + return self._array_shape + + @array_shape.setter + def array_shape(self, value): + self._array_shape = value + + @property + def pixel_bounds(self): + """ + The bounds (in pixel coordinates) inside which the WCS is defined, + as a list with `~BaseLowLevelWCS.pixel_n_dim` ``(min, max)`` tuples. + The bounds should be given in ``[(xmin, xmax), (ymin, ymax)]`` + order. WCS solutions are sometimes only guaranteed to be accurate + within a certain range of pixel values, for example when defining a + WCS that includes fitted distortions. This is an optional property, + and it should return `None` if a shape is not known or relevant. + """ + return self.bounding_box + + @property + def axis_correlation_matrix(self): + """ + Returns an (`~BaseLowLevelWCS.world_n_dim`, + `~BaseLowLevelWCS.pixel_n_dim`) matrix that indicates using booleans + whether a given world coordinate depends on a given pixel coordinate. + This defaults to a matrix where all elements are `True` in the absence of + any further information. For completely independent axes, the diagonal + would be `True` and all other entries `False`. + """ + + return separable.separability_matrix(self.forward_transform) + + def serialized_classes(self): + """ + Indicates whether Python objects are given in serialized form or as + actual Python objects. + """ + return False + + def world_axis_object_classes(self): + raise NotImplementedError() + + def world_axis_object_components(self): + raise NotImplementedError() + + + # High level APE 14 API + + def low_level_wcs(self): + """ + Returns a reference to the underlying low-level WCS object. + """ + return self + + def pixel_to_world(self, *pixel_arrays): + """ + Convert pixel values to world coordinates. + """ + if not self.forward_transform.uses_quantity: + kwargs = {'with_units': True} + return self(*pixel_arrays, **kwargs) + + def array_index_to_world(self, *index_arrays): + """ + Convert array indices to world coordinates (represented by Astropy + objects). + """ + return self(*(index_arrays[::-1]), with_units=True) + + def world_to_pixel(self, *world_objects): + """ + Convert world coordinates to pixel values. + """ + return self.invert(*world_objects) + + def world_to_array_index(self, *world_objects): + """ + Convert world coordinates (represented by Astropy objects) to array + indices. + """ + result = self.invert(*world_objects, with_units=True)[::-1] + return tuple([utils._toindex(r) for r in result]) + diff --git a/gwcs/coordinate_frames.py b/gwcs/coordinate_frames.py index f961d92e..25133ed2 100644 --- a/gwcs/coordinate_frames.py +++ b/gwcs/coordinate_frames.py @@ -4,10 +4,13 @@ """ import numpy as np +from astropy.utils.misc import isiterable import astropy.time from astropy import units as u from astropy import utils as astutil from astropy import coordinates as coord +from astropy.wcs.wcsapi.low_level_api import (validate_physical_types, + VALID_UCDS) __all__ = ['Frame2D', 'CelestialFrame', 'SpectralFrame', 'CompositeFrame', @@ -47,7 +50,7 @@ class CoordinateFrame: def __init__(self, naxes, axes_type, axes_order, reference_frame=None, reference_position=None, unit=None, axes_names=None, - name=None): + name=None, axis_physical_types=None): self._naxes = naxes self._axes_order = tuple(axes_order) if isinstance(axes_type, str): @@ -65,6 +68,8 @@ def __init__(self, naxes, axes_type, axes_order, reference_frame=None, raise ValueError("Number of units does not match number of axes.") else: self._unit = tuple([u.Unit(au) for au in unit]) + else: + self._unit = tuple("" for na in range(naxes)) if axes_names is not None: if isinstance(axes_names, str): axes_names = (axes_names,) @@ -75,6 +80,7 @@ def __init__(self, naxes, axes_type, axes_order, reference_frame=None, else: axes_names = tuple([""] * naxes) self._axes_names = axes_names + if name is None: self._name = self.__class__.__name__ else: @@ -88,6 +94,60 @@ def __init__(self, naxes, axes_type, axes_order, reference_frame=None, raise ValueError("Length of axes_order does not match number of axes.") super(CoordinateFrame, self).__init__() + self._axis_physical_types = self._set_axis_physical_types(axis_physical_types) + + def _set_axis_physical_types(self, pht=None): + """ + Set the physical type of the coordinate axes using VO UCD1+ v1.23 definitions. + + """ + if pht is not None: + if isinstance(pht, str): + pht = (pht,) + elif not isiterable(pht): + raise TypeError("axis_physical_types must be of type string or iterable of strings") + if len(pht) != self.naxes: + raise ValueError('"axis_physical_types" must be of length {}'.format(self.naxes)) + ph_type = [] + for axt in pht: + if axt not in VALID_UCDS and not axt.startswith("custom:"): + ph_type.append("custom:{}".format(axt)) + else: + ph_type.append(axt) + validate_physical_types(ph_type) + return tuple(ph_type) + + if isinstance(self, CelestialFrame): + if isinstance(self.reference_frame, coord.Galactic): + ph_type = "pos.galactic.lon", "pos.galactic.lat" + elif isinstance(self.reference_frame, (coord.GeocentricTrueEcliptic, + coord.GCRS, + coord.PrecessedGeocentric)): + ph_type = "pos.bodyrc.lon", "pos.bodyrc.lat" + elif isinstance(self.reference_frame, coord.builtin_frames.BaseRADecFrame): + ph_type = "pos.eq.ra", "pos.eq.dec" + elif isinstance(self.reference_frame, coord.builtin_frames.BaseEclipticFrame): + ph_type = "pos.ecliptic.lon", "pos.ecliptic.lat" + else: + ph_type = tuple("custom:{}".format(t) for t in self.axes_names) + elif isinstance(self, SpectralFrame): + if self.unit[0].physical_type == "frequency": + ph_type = ("em.freq",) + elif self.unit[0].physical_type == "length": + ph_type = ("em.wl",) + elif isinstance(self, TemporalFrame): + ph_type = ("time",) + elif isinstance(self, Frame2D): + if all(self.axes_names): + ph_type = self.axes_names + else: + ph_type = self.axes_type + ph_type = tuple("custom:{}".format(t) for t in ph_type) + else: + ph_type = tuple("custom:{}".format(t) for t in self.axes_type) + validate_physical_types(ph_type) + return ph_type + def __repr__(self): fmt = '<{0}(name="{1}", unit={2}, axes_names={3}, axes_order={4}'.format( @@ -153,7 +213,9 @@ def axes_type(self): def coordinates(self, *args): """ Create world coordinates object""" - raise NotImplementedError("Subclasses may implement this") + args = [args[i] for i in self.axes_order] + coo = tuple([arg * un if not hasattr(arg, "to") else arg.to(un) for arg, un in zip(args, self.unit)]) + return coo def coordinate_to_quantity(self, *coords): """ @@ -162,6 +224,10 @@ def coordinate_to_quantity(self, *coords): # NoOp leaves it to the model to handle return coords + @property + def axis_physical_types(self): + return self._axis_physical_types + class CelestialFrame(CoordinateFrame): """ @@ -173,8 +239,6 @@ class CelestialFrame(CoordinateFrame): A dimension in the input data that corresponds to this axis. reference_frame : astropy.coordinates.builtin_frames A reference frame. - reference_position : str - Reference position. unit : str or units.Unit instance or iterable of those Units on axes. axes_names : list @@ -185,7 +249,7 @@ class CelestialFrame(CoordinateFrame): def __init__(self, axes_order=None, reference_frame=None, unit=None, axes_names=None, - name=None): + name=None, axis_physical_types=None): naxes = 2 if reference_frame is not None: if not isinstance(reference_frame, str): @@ -205,12 +269,13 @@ def __init__(self, axes_order=None, reference_frame=None, if unit is None: unit = tuple([u.degree] * naxes) axes_type = ['SPATIAL'] * naxes + super(CelestialFrame, self).__init__(naxes=naxes, axes_type=axes_type, axes_order=axes_order, reference_frame=reference_frame, unit=unit, axes_names=axes_names, - name=name) + name=name, axis_physical_types = axis_physical_types) def coordinates(self, *args): """ @@ -221,8 +286,10 @@ def coordinates(self, *args): args : float inputs to wcs.input_frame """ - # Reorder axes if necessary. - return coord.SkyCoord(*args, unit=self.unit, frame=self._reference_frame) + if isinstance(args[0], coord.SkyCoord): + return args[0].transform_to(self.reference_frame) + else: + return coord.SkyCoord(*args, unit=self.unit, frame=self.reference_frame) def coordinate_to_quantity(self, *coords): if len(coords) == 2: @@ -232,7 +299,7 @@ def coordinate_to_quantity(self, *coords): else: raise ValueError("Unexpected number of coordinates in " "input to frame {} : " - "expected 2, got {}".format(name, len(coords))) + "expected 2, got {}".format(self.name, len(coords))) if isinstance(arg, coord.SkyCoord): arg = arg.transform_to(self._reference_frame) @@ -268,18 +335,23 @@ class SpectralFrame(CoordinateFrame): Spectral axis name. name : str Name for this frame. + reference_position : str + Reference position - one of `STANDARD_REFERENCE_POSITION` + """ def __init__(self, axes_order=(0,), reference_frame=None, unit=None, - axes_names=None, name=None, reference_position=None): + axes_names=None, name=None, axis_physical_types=None, + reference_position=None): super(SpectralFrame, self).__init__(naxes=1, axes_type="SPECTRAL", axes_order=axes_order, axes_names=axes_names, reference_frame=reference_frame, unit=unit, name=name, - reference_position=reference_position) - - def coordinates(self, *args): + reference_position=reference_position, + axis_physical_types=axis_physical_types) + + def coordinates(self, *args, equivalencies=[]): if hasattr(args[0], 'unit'): - return args[0] + return args[0].to(self.unit[0], equivalencies=equivalencies) if np.isscalar(args): return args * self.unit[0] else: @@ -317,12 +389,12 @@ class TemporalFrame(CoordinateFrame): def __init__(self, axes_order=(0,), reference_time=None, reference_frame=astropy.time.Time, unit=None, - axes_names=None, name=None): + axes_names=None, name=None, axis_physical_types=None): super().__init__(naxes=1, axes_type="TIME", axes_order=axes_order, axes_names=axes_names, reference_frame=reference_frame, unit=unit, name=name, - reference_position=reference_time) + reference_position=reference_time, axis_physical_types=axis_physical_types) def coordinates(self, *args): if np.isscalar(args): @@ -374,23 +446,27 @@ def __init__(self, frames, name=None): unit = list(range(naxes)) axes_names = list(range(naxes)) axes_order = [] + ph_type = list(range(naxes)) for frame in frames: axes_order.extend(frame.axes_order) for frame in frames: - for ind, axtype, un, n in zip(frame.axes_order, frame.axes_type, - frame.unit, frame.axes_names): + for ind, axtype, un, n, pht in zip(frame.axes_order, frame.axes_type, + frame.unit, frame.axes_names, frame.axis_physical_types): axes_type[ind] = axtype axes_names[ind] = n unit[ind] = un + ph_type[ind] = pht if len(np.unique(axes_order)) != len(axes_order): raise ValueError("Incorrect numbering of axes, " "axes_order should contain unique numbers, " "got {}.".format(axes_order)) + super(CompositeFrame, self).__init__(naxes, axes_type=axes_type, axes_order=axes_order, unit=unit, axes_names=axes_names, name=name) - + self._axis_physical_types = tuple(ph_type) + @property def frames(self): return self._frames @@ -400,9 +476,13 @@ def __repr__(self): def coordinates(self, *args): coo = [] - for frame in self.frames: - fargs = [args[i] for i in frame.axes_order] - coo.append(frame.coordinates(*fargs)) + if len(args) == len(self.frames): + for frame, arg in zip(self.frames, args): + coo.append(frame.coordinates(arg)) + else: + for frame in self.frames: + fargs = [args[i] for i in frame.axes_order] + coo.append(frame.coordinates(*fargs)) return coo def coordinate_to_quantity(self, *coords): @@ -476,11 +556,12 @@ class Frame2D(CoordinateFrame): """ def __init__(self, axes_order=(0, 1), unit=(u.pix, u.pix), axes_names=('x', 'y'), - name=None): + name=None, axis_physical_types=None): super(Frame2D, self).__init__(naxes=2, axes_type=["SPATIAL", "SPATIAL"], axes_order=axes_order, name=name, - axes_names=axes_names, unit=unit) + axes_names=axes_names, unit=unit, + axis_physical_types=axis_physical_types) def coordinates(self, *args): args = [args[i] for i in self.axes_order] @@ -496,7 +577,7 @@ def coordinate_to_quantity(self, *coords): else: raise ValueError("Unexpected number of coordinates in " "input to frame {} : " - "expected 2, got {}".format(name, len(coords))) + "expected 2, got {}".format(self.name, len(coords))) for i in range(2): if not hasattr(coords[i], 'unit'): diff --git a/gwcs/gwcs_types.py b/gwcs/gwcs_types.py index 88bee452..d4acf749 100644 --- a/gwcs/gwcs_types.py +++ b/gwcs/gwcs_types.py @@ -9,7 +9,7 @@ import six from astropy.io.misc.asdf.tags.transform.basic import TransformType -from asdf.asdftypes import ExtensionTypeMeta, CustomType +from asdf.types import ExtensionTypeMeta, CustomType from astropy.io.misc.asdf.types import AstropyTypeMeta diff --git a/gwcs/schemas/stsci.edu/gwcs/frame-1.0.0.yaml b/gwcs/schemas/stsci.edu/gwcs/frame-1.0.0.yaml index f02c066e..2544dcd9 100644 --- a/gwcs/schemas/stsci.edu/gwcs/frame-1.0.0.yaml +++ b/gwcs/schemas/stsci.edu/gwcs/frame-1.0.0.yaml @@ -68,5 +68,15 @@ properties: items: $ref: "tag:stsci.edu:asdf/unit/unit-1.0.0" + axis_physical_types: + description: | + An iterable of strings describing the physical type for each world axis. + These should be names from the VO UCD1+ controlled Vocabulary + (http://www.ivoa.net/documents/latest/UCDlist.html). + type: array + items: + type: + string + required: [name] additionalProperties: true diff --git a/gwcs/schemas/stsci.edu/gwcs/label_mapper-1.0.0.yaml b/gwcs/schemas/stsci.edu/gwcs/label_mapper-1.0.0.yaml index 1dddea37..5a95344f 100644 --- a/gwcs/schemas/stsci.edu/gwcs/label_mapper-1.0.0.yaml +++ b/gwcs/schemas/stsci.edu/gwcs/label_mapper-1.0.0.yaml @@ -11,7 +11,7 @@ description: | [regions_selector](ref:http://stsci.edu/schemas/gwcs/regions_selector-1.0.0). The [label_mapper](ref:http://stsci.edu/schemas/gwcs/label_mapper-1.0.0) returns the label corresponding to given inputs. The - [regions_selector](ref:http://stsci.edu/schemas/gwcsregions_selector-1.0.0) + [regions_selector](ref:http://stsci.edu/schemas/gwcs/regions_selector-1.0.0) returns the transform corresponding to this label. This maps inputs (e.g. pixels on a detector) to transforms uniquely. diff --git a/gwcs/selector.py b/gwcs/selector.py index 3c0dde55..26245cce 100644 --- a/gwcs/selector.py +++ b/gwcs/selector.py @@ -189,7 +189,7 @@ def __init__(self, mapper, inputs_mapping=None, name=None, **kwargs): super(LabelMapperArray, self).__init__(mapper, _no_label, name=name, **kwargs) def evaluate(self, *args): - args = [_toindex(a) for a in args] + args = tuple([_toindex(a) for a in args]) try: result = self._mapper[args[::-1]] except IndexError as e: diff --git a/gwcs/tags/wcs.py b/gwcs/tags/wcs.py index 9720b62c..4ab306fc 100644 --- a/gwcs/tags/wcs.py +++ b/gwcs/tags/wcs.py @@ -98,6 +98,9 @@ def _from_tree(cls, node, ctx): kwargs['unit'] = tuple( yamlutil.tagged_tree_to_custom_tree(node['unit'], ctx)) + if 'axis_physical_type' in node: + kwargs['axis_physical_type'] = tuple(node['axis_physical_type']) + return kwargs @classmethod @@ -112,7 +115,7 @@ def _to_tree(cls, frame, ctx): node['axes_type'] = frame.axes_type node['naxes'] = frame.naxes - if frame.axes_order != (0, 1): + if frame.axes_order is not None: node['axes_order'] = list(frame.axes_order) if frame.axes_names is not None: @@ -126,6 +129,9 @@ def _to_tree(cls, frame, ctx): node['unit'] = yamlutil.custom_tree_to_tagged_tree( list(frame.unit), ctx) + if frame.axis_physical_types is not None: + node['axis_physical_types'] = list(frame.axis_physical_types) + return node @classmethod diff --git a/gwcs/tests/setup_package.py b/gwcs/tests/setup_package.py index f9b6ff38..edc713e6 100644 --- a/gwcs/tests/setup_package.py +++ b/gwcs/tests/setup_package.py @@ -2,5 +2,5 @@ def get_package_data(): return { - _ASTROPY_PACKAGE_NAME_ + '.tests': ['coveragerc', 'data/*.hdr'] + _ASTROPY_PACKAGE_NAME_ + '.tests': ['coveragerc', 'data/*.hdr'] # noqa } diff --git a/gwcs/tests/test_api.py b/gwcs/tests/test_api.py new file mode 100644 index 00000000..6c45a310 --- /dev/null +++ b/gwcs/tests/test_api.py @@ -0,0 +1,183 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +""" +Tests the API defined in astropy APE 14 (https://doi.org/10.5281/zenodo.1188875). +""" +import numpy as np +from numpy.testing import assert_allclose, assert_array_equal + +from astropy import coordinates as coord +from astropy.modeling import models +from astropy import units as u + +from .. import wcs +from .. import coordinate_frames as cf + +import pytest + +# frames +sky_frame = cf.CelestialFrame(reference_frame=coord.ICRS(), axes_order=(0, 1)) +detector = cf.Frame2D(name='detector', axes_order=(0, 1)) + +spec1 = cf.SpectralFrame(name='freq', unit=[u.Hz, ], axes_order=(2, )) +spec2 = cf.SpectralFrame(name='wave', unit=[u.m, ], axes_order=(2, ), axes_names=('lambda', )) + +comp1 = cf.CompositeFrame([sky_frame, spec1]) + +# transforms + +m1 = models.Shift(1) & models.Shift(2) +m2 = models.Scale(2) +m = m1 & m2 + +pipe = [(detector, m1), + (sky_frame, None) + ] + +example_wcs = wcs.WCS(pipe) + +def create_example_wcs(): + example_wcs = [wcs.WCS([(detector, m1), + (sky_frame, None)]), + wcs.WCS([(detector, m2), + (spec1, None)]), + wcs.WCS([(detector, m), + (comp1, None)]) + ] + + pixel_world_ndim = [(2, 2), (2, 1), (2, 3)] + physical_types = [("pos.eq.ra", "pos.eq.dec"), ("em.freq",), ("pos.eq.ra", "pos.eq.dec", "em.freq")] + world_units = [("deg", "deg"), ("Hz",), ("deg", "deg", "Hz")] + + return example_wcs, pixel_world_ndim, physical_types, world_units + +# x, y inputs - scalar and array +x, y = 1, 2 +xarr, yarr = np.ones((3, 4)), np.ones((3, 4)) + 1 + +# ra, dec inputs - scalar, arrays and SkyCoord objects +ra, dec = 2, 4 +sky = coord.SkyCoord(ra * u.deg, dec*u.deg, frame = sky_frame.reference_frame) +raarr = np.ones((3, 4)) * ra +decarr = np.ones((3, 4)) * dec +skyarr = coord.SkyCoord(raarr * u.deg, decarr*u.deg, frame = sky_frame.reference_frame) + +ex_wcs, dims, physical_types, world_units = create_example_wcs() + +@pytest.mark.parametrize(("wcsobj", "ndims"), zip(ex_wcs, dims)) +def test_pixel_n_dim(wcsobj, ndims): + assert wcsobj.pixel_n_dim == ndims[0] + + +@pytest.mark.parametrize(("wcsobj", "ndims"), zip(ex_wcs, dims)) +def test_world_n_dim(wcsobj, ndims): + assert wcsobj.world_n_dim == ndims[1] + + +@pytest.mark.parametrize(("wcsobj", "physical_types"), zip(ex_wcs, physical_types)) +def test_world_axis_physical_types(wcsobj, physical_types): + assert wcsobj.world_axis_physical_types == physical_types + + +@pytest.mark.parametrize(("wcsobj", "world_units"), zip(ex_wcs, world_units)) +def test_world_axis_units(wcsobj, world_units): + assert wcsobj.world_axis_units == world_units + + +@pytest.mark.parametrize(("x", "y"), zip((x, xarr), (y, yarr))) +def test_pixel_to_world_values(x, y): + wcsobj = example_wcs + assert_allclose(wcsobj.pixel_to_world_values(x, y), wcsobj(x, y, with_units=False)) + + +@pytest.mark.parametrize(("x", "y"), zip((x, xarr), (y, yarr))) +def test_array_index_to_world_values(x, y): + wcsobj = example_wcs + assert_allclose(wcsobj.array_index_to_world_values(x, y), wcsobj(y, x, with_units=False)) + + +@pytest.mark.parametrize(("sky", "ra", "dec"), zip((sky, skyarr), (ra, raarr), (dec, decarr))) +def test_world_to_pixel_values(sky, ra, dec): + wcsobj = example_wcs + assert_allclose(wcsobj.world_to_pixel_values(sky), wcsobj.invert(ra, dec, with_units=False)) + + +@pytest.mark.parametrize(("sky", "ra", "dec"), zip((sky, skyarr), (ra, raarr), (dec, decarr))) +def test_world_to_array_index_values(sky, ra, dec): + wcsobj = example_wcs + assert_allclose(wcsobj.world_to_array_index_values(sky), + wcsobj.invert(ra, dec, with_units=False)[::-1]) + + +def test_world_axis_object_components(): + wcsobj = example_wcs + with pytest.raises(NotImplementedError): + wcsobj.world_axis_object_components() + + +def test_world_axis_object_classes(): + wcsobj = example_wcs + with pytest.raises(NotImplementedError): + wcsobj.world_axis_object_classes() + + +def test_array_shape(): + wcsobj = example_wcs + assert wcsobj.array_shape is None + + wcsobj.array_shape = (2040, 1020) + assert_array_equal(wcsobj.array_shape, (2040, 1020)) + + +def test_pixel_bounds(): + wcsobj = example_wcs + assert wcsobj.pixel_bounds is None + + wcsobj.bounding_box = ((-0.5, 2039.5), (-0.5, 1019.5)) + assert_array_equal(wcsobj.pixel_bounds, wcsobj.bounding_box) + + +def test_axis_correlation_matrix(): + wcsobj = example_wcs + assert_array_equal(wcsobj.axis_correlation_matrix, np.identity(2)) + + +def test_serialized_classes(): + wcsobj = example_wcs + assert wcsobj.serialized_classes() == False + + +def test_low_level_wcs(): + wcsobj = example_wcs + assert id(wcsobj.low_level_wcs()) == id(wcsobj) + + +def test_pixel_to_world(): + wcsobj = example_wcs + comp = wcsobj(x, y, with_units=True) + comp = wcsobj.output_frame.coordinates(comp) + result = wcsobj.pixel_to_world(x, y) + assert isinstance(comp, coord.SkyCoord) + assert isinstance(result, coord.SkyCoord) + assert_allclose(comp.data.lon, result.data.lon) + assert_allclose(comp.data.lat, result.data.lat) + + +def test_array_index_to_world(): + wcsobj = example_wcs + comp = wcsobj(x, y, with_units=True) + comp = wcsobj.output_frame.coordinates(comp) + result = wcsobj.array_index_to_world(y, x) + assert isinstance(comp, coord.SkyCoord) + assert isinstance(result, coord.SkyCoord) + assert_allclose(comp.data.lon, result.data.lon) + assert_allclose(comp.data.lat, result.data.lat) + + +def test_world_to_pixel(): + wcsobj = example_wcs + assert_allclose(wcsobj.world_to_pixel(sky), wcsobj.invert(ra, dec, with_units=False)) + + +def test_world_to_array_index(): + wcsobj = example_wcs + assert_allclose(wcsobj.world_to_array_index(sky), wcsobj.invert(ra, dec, with_units=False)[::-1]) diff --git a/gwcs/tests/test_coordinate_systems.py b/gwcs/tests/test_coordinate_systems.py index 0082174c..557848a1 100644 --- a/gwcs/tests/test_coordinate_systems.py +++ b/gwcs/tests/test_coordinate_systems.py @@ -19,7 +19,7 @@ # For now remove OffsetFrame, issue #55 try: coord_frames.remove("SkyOffsetFrame") -except ValueError as ve: +except ValueError: pass @@ -84,9 +84,9 @@ def test_celestial_attributes_length(frame): CelestialFrame attributes from reference_frame. """ fr = getattr(coord, frame) - if issubclass(fr, coord.BaseCoordinateFrame): + if issubclass(fr.__class__, coord.BaseCoordinateFrame): cel = cf.CelestialFrame(reference_frame=fr()) - assert(len(cel.axes_names) == len(cel.axes_type) == len(cel.unit) == + assert(len(cel.axes_names) == len(cel.axes_type) == len(cel.unit) == \ len(cel.axes_order) == cel.naxes) @@ -174,7 +174,13 @@ def test_coordinate_to_quantity_celestial(inp): assert_quantity_allclose(lon, 10 * u.deg) assert_quantity_allclose(lat, 20 * u.deg) + with pytest.raises(ValueError): + cel.coordinate_to_quantity(10 * u.deg, 2 * u.deg, 3 * u.deg) + + with pytest.raises(ValueError): + cel.coordinate_to_quantity((1, 2)) + @pytest.mark.parametrize('inp', [ (100,), (100 * u.nm,), @@ -272,3 +278,70 @@ def test_coordinate_to_quantity_frame_2d(): result = frame.coordinate_to_quantity(*inp) for output, exp in zip(result, expected): assert_quantity_allclose(output, exp) + + +def test_coordinate_to_quantity_error(): + frame = cf.Frame2D(unit=(u.one, u.arcsec)) + with pytest.raises(ValueError): + frame.coordinate_to_quantity(1) + + with pytest.raises(ValueError): + comp1.coordinate_to_quantity((1, 1), 2) + + frame = cf.TemporalFrame(unit=u.s) + with pytest.raises(ValueError): + frame.coordinate_to_quantity(1) + + +def test_axis_physical_type(): + assert icrs.axis_physical_types == ("pos.eq.ra", "pos.eq.dec") + assert spec1.axis_physical_types == ("em.freq",) + assert spec2.axis_physical_types == ("em.wl",) + assert comp1.axis_physical_types == ("pos.eq.ra", "pos.eq.dec", "em.freq") + assert comp2.axis_physical_types == ("custom:x", "custom:y", "em.wl") + assert comp.axis_physical_types == ('pos.eq.ra', 'pos.eq.dec', 'em.freq', 'em.wl') + + spec3 = cf.SpectralFrame(name='waven', axes_order=(1,), + axis_physical_types='em.wavenumber') + assert spec3.axis_physical_types == ('em.wavenumber',) + + t = cf.TemporalFrame(reference_time=Time("2018-01-01T00:00:00"), unit=u.s) + assert t.axis_physical_types == ('time',) + + fr2d = cf.Frame2D(name='d', axes_names=("x", "y")) + assert fr2d.axis_physical_types == ('custom:x', 'custom:y') + + fr2d = cf.Frame2D(name='d', axes_names=None) + assert fr2d.axis_physical_types == ('custom:SPATIAL', 'custom:SPATIAL') + + fr2d = cf.Frame2D(name='d', axis_physical_types=("pos.x", "pos.y")) + assert fr2d.axis_physical_types == ('custom:pos.x', 'custom:pos.y') + + with pytest.raises(ValueError): + cf.CelestialFrame(reference_frame=coord.ICRS(), axis_physical_types=("pos.eq.ra",)) + + fr = cf.CelestialFrame(reference_frame=coord.ICRS(), axis_physical_types=("ra", "dec")) + assert fr.axis_physical_types == ("custom:ra", "custom:dec") + + fr = cf.CelestialFrame(reference_frame=coord.BarycentricTrueEcliptic()) + assert fr.axis_physical_types == ('pos.ecliptic.lon', 'pos.ecliptic.lat') + + frame = cf.CoordinateFrame(name='custom_frame', axes_type=("SPATIAL",), axes_order=(0,), axis_physical_types="length", axes_names="x", naxes=1) + assert frame.axis_physical_types == ("custom:length",) + frame = cf.CoordinateFrame(name='custom_frame', axes_type=("SPATIAL",), axes_order=(0,), axis_physical_types=("length",), axes_names="x", naxes=1) + assert frame.axis_physical_types == ("custom:length",) + with pytest.raises(ValueError): + cf.CoordinateFrame(name='custom_frame', axes_type=("SPATIAL",), axes_order=(0,), axis_physical_types=("length", "length"), naxes=1) + + +def test_base_frame(): + with pytest.raises(ValueError): + cf.CoordinateFrame(name='custom_frame', + axes_type=("SPATIAL",), + naxes=1, axes_order=(0,), + axes_names=("x", "y")) + frame = cf.CoordinateFrame(name='custom_frame', axes_type=("SPATIAL",), axes_order=(0,), axes_names="x", naxes=1) + assert frame.naxes == 1 + assert frame.axes_names == ("x",) + + frame.coordinate_to_quantity(1, 2) diff --git a/gwcs/tests/test_region.py b/gwcs/tests/test_region.py index e39fd332..7e4d2ffb 100644 --- a/gwcs/tests/test_region.py +++ b/gwcs/tests/test_region.py @@ -1,3 +1,4 @@ + # Licensed under a 3-clause BSD style license - see LICENSE.rst """ Test regions @@ -37,7 +38,7 @@ def test_LabelMapperArray_from_vertices_string(): assert(np.sort(labels) == np.sort(mask_labels)).all() -#### These tests below check the scanning algorithm for two shapes ########## +# These tests below check the scanning algorithm for two shapes def polygon1(shape=(9, 9)): ar = np.zeros(shape) ar[1, 2] = 1 diff --git a/gwcs/tests/test_utils.py b/gwcs/tests/test_utils.py index 18e32415..0d1ac5c6 100644 --- a/gwcs/tests/test_utils.py +++ b/gwcs/tests/test_utils.py @@ -29,7 +29,7 @@ def test_fits_transform(): hdr = fits.Header.fromfile(get_pkg_data_filename('data/simple_wcs2.hdr')) gw1 = gwutils.make_fitswcs_transform(hdr) w1 = fitswcs.WCS(hdr) - assert_allclose(gw1(1, 2), w1.wcs_pix2world(1, 2, 1), atol=10**-8) + assert_allclose(gw1(1, 2), w1.wcs_pix2world(1, 2, 1), atol=10 ** -8) def test_lon_pole(): @@ -86,10 +86,10 @@ def test_get_axes(): def test_isnumerical(): - sky = coord.SkyCoord(1*u.deg, 2*u.deg) + sky = coord.SkyCoord(1 * u.deg, 2 * u.deg) assert not gwutils.isnumerical(sky) - assert not gwutils.isnumerical(2*u.m) + assert not gwutils.isnumerical(2 * u.m) assert gwutils.isnumerical(np.float(0)) assert gwutils.isnumerical(np.array(0)) diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py index 83ff4c75..8ebe7dda 100644 --- a/gwcs/tests/test_wcs.py +++ b/gwcs/tests/test_wcs.py @@ -1,6 +1,7 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst +import warnings import numpy as np -from numpy.testing import assert_allclose +from numpy.testing import assert_allclose, assert_equal from astropy.modeling import models from astropy import coordinates as coord from astropy.io import fits @@ -10,7 +11,7 @@ from astropy import wcs as astwcs from .. import wcs -from ..wcstools import * +from ..wcstools import (wcs_from_fiducial, grid_from_bounding_box, wcs_from_points) from .. import coordinate_frames as cf from .. import utils from ..utils import CoordinateFrameError @@ -30,8 +31,15 @@ (icrs, None) ] +# Create some data. +nx, ny = (5, 2) +x = np.linspace(0, 1, nx) +y = np.linspace(0, 1, ny) +xv, yv = np.meshgrid(x, y) + # Test initializing a WCS + def test_create_wcs(): """ Test initializing a WCS object. @@ -45,7 +53,7 @@ def test_create_wcs(): # use a pipeline to initialize pipe = [(detector, m1), (icrs, None)] gw4 = wcs.WCS(forward_transform=pipe) - assert(gw1.available_frames == gw2.available_frames == + assert(gw1.available_frames == gw2.available_frames == \ gw3.available_frames == gw4.available_frames == ['detector', 'icrs']) res = m(1, 2) assert_allclose(gw1(1, 2), res) @@ -86,7 +94,7 @@ def test_insert_transform(): def test_set_transform(): """ Test setting a transform between two frames in the pipeline.""" - w = wcs.WCS(input_frame=detector, output_frame=icrs, forward_transform=pipe) + w = wcs.WCS(forward_transform=pipe[:]) w.set_transform('detector', 'focal', models.Identity(2)) assert_allclose(w(1, 1), (2, -2)) with pytest.raises(CoordinateFrameError): @@ -97,7 +105,7 @@ def test_set_transform(): def test_get_transform(): """ Test getting a transform between two frames in the pipeline.""" - w = wcs.WCS(pipe) + w = wcs.WCS(pipe[:]) tr_forward = w.get_transform('detector', 'focal') tr_back = w.get_transform('icrs', 'detector') x, y = 1, 2 @@ -126,7 +134,7 @@ def test_backward_transform(): def test_return_coordinates(): """Test converting to coordinate objects or quantities.""" - w = wcs.WCS(pipe) + w = wcs.WCS(pipe[:]) x = 1 y = 2.3 numerical_result = (26.8, -0.6) @@ -244,6 +252,32 @@ def test_grid_from_bounding_box_step(): grid_from_bounding_box(bb, step=(1, 2, 1)) +def test_wcs_from_points(): + np.random.seed(0) + hdr = fits.Header.fromtextfile(get_pkg_data_filename("data/acs.hdr"), endcard=False) + with warnings.catch_warnings() as w: + warnings.simplefilter("ignore") + w = astwcs.WCS(hdr) + y, x = np.mgrid[:2046:20j, :4023:10j] + ra, dec = w.wcs_pix2world(x, y, 1) + fiducial = coord.SkyCoord(ra.mean()*u.deg, dec.mean()*u.deg, frame="icrs") + w = wcs_from_points(xy=(x, y), world_coordinates=(ra, dec), fiducial=fiducial) + newra, newdec = w(x, y) + assert_allclose(newra, ra) + assert_allclose(newdec, dec) + + n = np.random.randn(ra.size) + n.shape = ra.shape + nra = n * 10 ** -2 + ndec = n * 10 ** -2 + w = wcs_from_points(xy=(x + nra, y + ndec), + world_coordinates=(ra, dec), + fiducial=fiducial) + newra, newdec = w(x, y) + assert_allclose(newra, ra, atol=10**-6) + assert_allclose(newdec, dec, atol=10**-6) + + def test_grid_from_bounding_box_2(): bb = ((-0.5, 5.5), (-0.5, 4.5)) x, y = grid_from_bounding_box(bb) @@ -295,6 +329,58 @@ def test_available_frames(): assert w.available_frames == ['detector', 'focal', 'icrs'] +def test_footprint(): + icrs = cf.CelestialFrame(name='icrs', reference_frame=coord.ICRS(), + axes_order=(0, 1)) + spec = cf.SpectralFrame(name='freq', unit=[u.Hz, ], axes_order=(2, )) + world = cf.CompositeFrame([icrs, spec]) + transform = (models.Shift(10) & models.Shift(-1)) & models.Scale(2) + pipe = [('det', transform), (world, None)] + w = wcs.WCS(pipe) + + with pytest.raises(TypeError): + w.footprint() + + w.bounding_box = ((1,5), (1,3), (1, 6)) + + assert_equal(w.footprint(), np.array([[11, 0, 2], + [11, 0, 12], + [11, 2, 2], + [11, 2, 12], + [15, 0, 2], + [15, 0, 12], + [15, 2, 2], + [15, 2, 12]])) + assert_equal(w.footprint(axis_type='spatial'), np.array([[11., 0.], + [11., 2.], + [15., 2.], + [15., 0.]])) + + assert_equal(w.footprint(axis_type='spectral'), np.array([2, 12])) + + +def test_high_level_api(): + """ + Test WCS high level API. + """ + output_frame = cf.CompositeFrame(frames=[icrs, spec]) + transform = m1 & models.Scale(1.5) + w = wcs.WCS(forward_transform=transform, output_frame=output_frame) + + r, d, lam = w(xv, yv, xv) + world_coord = w.pixel_to_world(xv, yv, xv) + assert isinstance(world_coord[0], coord.SkyCoord) + assert isinstance(world_coord[1], u.Quantity) + assert_allclose(world_coord[0].data.lon.value, r) + assert_allclose(world_coord[0].data.lat.value, d) + assert_allclose(world_coord[1].value, lam) + + x1, y1, z1 = w.world_to_pixel(*world_coord) + assert_allclose(x1, xv) + assert_allclose(y1, yv) + assert_allclose(z1, xv) + + class TestImaging(object): def setup_class(self): hdr = fits.Header.fromtextfile(get_pkg_data_filename("data/acs.hdr"), endcard=False) @@ -333,10 +419,7 @@ def setup_class(self): self.wcs = wcs.WCS(input_frame=det, output_frame=sky_cs, forward_transform=pipeline) - nx, ny = (5, 2) - x = np.linspace(0, 1, nx) - y = np.linspace(0, 1, ny) - self.xv, self.yv = np.meshgrid(x, y) + self.xv, self.yv = xv, yv @pytest.mark.filterwarnings('ignore') def test_distortion(self): @@ -368,7 +451,7 @@ def test_backward(self): def test_footprint(self): bb = ((1, 4096), (1, 2048)) - footprint = (self.wcs.footprint(bb)).T + footprint = (self.wcs.footprint(bb)) fits_footprint = self.fitsw.calc_footprint(axes=(4096, 2048)) assert_allclose(footprint, fits_footprint) @@ -387,5 +470,12 @@ def test_units(self): def test_get_transform(self): with pytest.raises(wcs.CoordinateFrameError): - assert(self.wcs.get_transform('x_translation', 'sky_rotation').submodel_names == + assert(self.wcs.get_transform('x_translation', 'sky_rotation').submodel_names == \ self.wcs.forward_transform[1:].submodel_names) + + def test_pixel_to_world(self): + sky_coord = self.wcs.pixel_to_world(self.xv, self.yv) + ra, dec = self.fitsw.all_pix2world(self.xv, self.yv, 1) + assert isinstance(sky_coord, coord.SkyCoord) + assert_allclose(sky_coord.data.lon.value, ra) + assert_allclose(sky_coord.data.lat.value, dec) diff --git a/gwcs/utils.py b/gwcs/utils.py index 515c4d1e..8f477689 100644 --- a/gwcs/utils.py +++ b/gwcs/utils.py @@ -168,7 +168,7 @@ def read_wcs_from_header(header): try: wcs_info['WCSAXES'] = header['WCSAXES'] except KeyError: - p = re.compile('ctype[\d]*', re.IGNORECASE) + p = re.compile(r'ctype[\d]*', re.IGNORECASE) ctypes = header['CTYPE*'] keys = list(ctypes.keys()) for key in keys[::-1]: @@ -466,50 +466,3 @@ def isnumerical(val): elif isinstance(val, np.ndarray) and val.dtype not in dtypes: isnum = False return isnum - - -def separable_axes(wcsobj, start_frame=None, end_frame=None): - """ - Computes the separability of axes in ``end_frame``. - - Returns a 1D boolean array of size frame.naxes where True means - the axis is completely separable and False means the axis is nonseparable - from at least one other axis. - - Parameters - ---------- - wcsobj : `~gwcs.wcs.WCS` - WCS object - start_frame : `~gwcs.coordinate_frames.CoordinateFrame` - A frame in the WCS pipeline. - The transform between start_frame and the end frame is used to compute the - mapping inputs: outputs. - If None the input_frame is used as start_frame. - end_frame : `~gwcs.coordinate_frames.CoordinateFrame` - A frame in the WCS pipeline. - The transform between start_frame and the end frame is used to compute the - mapping inputs: outputs. - If None wcsobj.output_frame is used. - - See Also - -------- - input_axes : For each output axis return the input axes contributing to it. - - """ - if wcsobj is not None: - if start_frame is None: - start_frame = wcsobj.input_frame - else: - if start_frame not in wcsobj.available_frames: - raise ValueError("Unrecognized frame {0}".format(start_frame)) - if end_frame is None: - end_frame = wcsobj.output_frame - else: - if end_frame not in wcsobj.available_frames: - raise ValueError("Unrecognized frame {0}".format(end_frame)) - transform = wcsobj.get_transform(start_frame, end_frame) - else: - raise ValueError("A starting frame is needed to determine separability of axes.") - - sep = is_separable(transform) - return [sep[ax] for ax in end_frame.axes_order] diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 2e696ee4..66de205c 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -1,10 +1,11 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst import functools - +import itertools import numpy as np from astropy.modeling.core import Model from astropy.modeling import utils as mutils +from .api import GWCSAPIMixin from . import coordinate_frames from .utils import CoordinateFrameError from .utils import _toindex @@ -14,7 +15,7 @@ __all__ = ['WCS'] -class WCS: +class WCS(GWCSAPIMixin): """ Basic WCS class. @@ -36,9 +37,11 @@ class WCS: def __init__(self, forward_transform=None, input_frame='detector', output_frame=None, name=""): + #self.low_level_wcs = self self._available_frames = [] self._pipeline = [] self._name = name + self._array_shape = None self._initialize_wcs(forward_transform, input_frame, output_frame) def _initialize_wcs(self, forward_transform, input_frame, output_frame): @@ -53,14 +56,14 @@ def _initialize_wcs(self, forward_transform, input_frame, output_frame): super(WCS, self).__setattr__(_input_frame, inp_frame_obj) super(WCS, self).__setattr__(_output_frame, outp_frame_obj) - self._pipeline = [(_input_frame, forward_transform.copy()), - (_output_frame, None)] + self._pipeline = [(input_frame, forward_transform.copy()), + (output_frame, None)] elif isinstance(forward_transform, list): for item in forward_transform: name, frame_obj = self._get_frame_name(item[0]) super(WCS, self).__setattr__(name, frame_obj) - self._pipeline.append((name, item[1])) - + #self._pipeline.append((name, item[1])) + self._pipeline = forward_transform else: raise TypeError("Expected forward_transform to be a model or a " "(frame, transform) list, got {0}".format( @@ -183,8 +186,9 @@ def _get_frame_index(self, frame): """ if isinstance(frame, coordinate_frames.CoordinateFrame): frame = frame.name - return np.asarray(self._pipeline)[:, 0].tolist().index(frame) - + frame_names = [getattr(item[0], "name", item[0]) for item in self._pipeline] + return frame_names.index(frame) + def _get_frame_name(self, frame): """ Return the name of the frame and a ``CoordinateFrame`` object. @@ -293,7 +297,7 @@ def invert(self, *args, **kwargs): except (NotImplementedError, KeyError): result = self._invert(*args, **kwargs) - if with_units: + if with_units and self.input_frame: if self.input_frame.naxes == 1: return self.input_frame.coordinates(result) else: @@ -367,7 +371,7 @@ def available_frames(self): {frame_name: frame_object or None} """ if self._pipeline: - return [frame[0] for frame in self._pipeline] + return [getattr(frame[0], "name", frame[0]) for frame in self._pipeline] else: return None @@ -401,7 +405,7 @@ def unit(self): """The unit of the coordinates in the output coordinate system.""" if self._pipeline: try: - return getattr(self, self._pipeline[-1][0]).unit + return getattr(self, self._pipeline[-1][0].name).unit except AttributeError: return None else: @@ -411,7 +415,10 @@ def unit(self): def output_frame(self): """Return the output coordinate frame.""" if self._pipeline: - return getattr(self, self._pipeline[-1][0]) + frame = self._pipeline[-1][0] + if not isinstance(frame, str): + frame = frame.name + return getattr(self, frame) else: return None @@ -419,7 +426,10 @@ def output_frame(self): def input_frame(self): """Return the input coordinate frame.""" if self._pipeline: - return getattr(self, self._pipeline[0][0]) + frame = self._pipeline[0][0] + if not isinstance(frame, str): + frame = frame.name + return getattr(self, frame) else: return None @@ -519,9 +529,9 @@ def __repr__(self): self.output_frame, self.input_frame, self.forward_transform) return fmt - def footprint(self, bounding_box=None, center=False): + def footprint(self, bounding_box=None, center=False, axis_type="all"): """ - Return the footprint of the observation in world coordinates. + Return the footprint in world coordinates. Parameters ---------- @@ -529,21 +539,55 @@ def footprint(self, bounding_box=None, center=False): `prop: bounding_box` center : bool If `True` use the center of the pixel, otherwise use the corner. + axis_type : str + A supported ``output_frame.axes_type`` or "all" (default). + One of ['spatial', 'spectral', 'temporal'] or a custom type. Returns ------- - coord : array of coordinates in the output_frame. - The order is counter-clockwise starting with the bottom left corner. + coord : ndarray + Array of coordinates in the output_frame mapping + corners to the output frame. For spatial coordinates the order + is clockwise, starting from the bottom left corner. + """ + def _order_clockwise(v): + return np.asarray([[v[0][0], v[1][0]], [v[0][0], v[1][1]], + [v[0][1], v[1][1]], [v[0][1], v[1][0]]]).T + 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 else: bb = bounding_box - vertices = np.asarray([[bb[0][0], bb[1][0]], [bb[0][0], bb[1][1]], - [bb[0][1], bb[1][1]], [bb[0][1], bb[1][0]]]).T + + all_spatial = all([t.lower() == "spatial" for t in self.output_frame.axes_type]) + + if all_spatial: + vertices = _order_clockwise(bb) + else: + vertices = np.array(list(itertools.product(*bb))).T + if center: vertices = _toindex(vertices) - result = self.__call__(*vertices, **{'with_bounding_box': False}) - return np.asarray(result) + + result = np.asarray(self.__call__(*vertices, **{'with_bounding_box': False})) + + axis_type = axis_type.lower() + if axis_type == 'spatial' and all_spatial: + return result.T + + if axis_type != "all": + 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 axis_type == "spatial": + result = _order_clockwise(result) + else: + result.sort() + result = np.squeeze(result) + + return result.T diff --git a/gwcs/wcstools.py b/gwcs/wcstools.py index 6d45951a..cfc53c9d 100644 --- a/gwcs/wcstools.py +++ b/gwcs/wcstools.py @@ -1,17 +1,21 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst import functools +import warnings import numpy as np from astropy.modeling.core import Model from astropy.modeling import projections -from astropy.modeling import models +from astropy.modeling import models, fitting from astropy import coordinates as coord +from astropy import units as u from .wcs import WCS -from .coordinate_frames import * +from .coordinate_frames import * # noqa from .utils import UnsupportedTransformError, UnsupportedProjectionError from .utils import _compute_lon_pole -__all__ = ['wcs_from_fiducial', 'grid_from_bounding_box'] + +__all__ = ['wcs_from_fiducial', 'grid_from_bounding_box', 'wcs_from_points'] + def wcs_from_fiducial(fiducial, coordinate_frame=None, projection=None, @@ -85,7 +89,7 @@ def wcs_from_fiducial(fiducial, coordinate_frame=None, projection=None, if len(bounding_box) != forward_transform.n_outputs: raise ValueError("Expected the number of items in 'bounding_box' to be equal to the " "number of outputs of the forawrd transform.") - forward_transform.bounding_box = bonding_box[::-1] + forward_transform.bounding_box = bounding_box[::-1] return WCS(output_frame=coordinate_frame, forward_transform=forward_transform, name=name) @@ -205,3 +209,72 @@ def _bbox_to_pixel(bbox): return grid[0] else: return grid + + +def wcs_from_points(xy, world_coordinates, fiducial, + projection=projections.Sky2Pix_TAN(), + degree=4, polynomial_type="polynomial"): + """ + Given two matching sets of coordinates on detector and sky, compute the WCS. + + Notes + ----- + This function implements the following algorithm: + ``world_coordinates`` are transformed to a projection plane using the specified projection. + A polynomial fits ``xy`` and the projected coordinates. + The fitted polynomials and the projection transforms are combined into a tranform + from detector to sky. + The input coordinate frame is set to ``detector``. + The output coordinate frame is initialized based on the frame in the fiducial. + + + Parameters + ---------- + xy : tuple of 2 ndarrays + Points in the input cooridnate frame - x, y inputs. + world_coordinates : tuple of 2 ndarrays + Points in the output coordinate frame. + The order matches the order of ``xy``. + fiducial_point : `~astropy.coordinates.SkyCoord` + A fiducial point in the output coordinate frame. + projection : `~astropy.modeling.projections.Projection` + A projection type. One of the projections in `~astropy.modeling.projections.projcode`. + degree : int + Degree of Polynpomial model to be fit to data. + polynomial_type : str + one of "polynomial", "chebyshev", "legendre" + + Returns + ------- + wcsobj : `~gwcs.wcs.WCS` + a WCS object for this observation. + """ + supported_poly_types = {"polynomial": models.Polynomial2D, + "chebyshev": models.Chebyshev2D, + "legendre": models.Legendre2D + } + x, y = xy + lon, lat = world_coordinates + + if not isinstance(projection, projections.Projection): + raise UnsupportedProjectionError("Unsupported projection code {0}".format(projection)) + if polynomial_type not in supported_poly_types.keys(): + raise ValueError("Unsupported polynomial_type: {}. " + "Only one of {} is supported.".format(polynomial_type, supported_poly_types.keys())) + skyrot = models.RotateCelestial2Native(fiducial.data.lon, fiducial.data.lat, 180*u.deg) + trans = (skyrot | projection) + projection_x, projection_y = trans(lon, lat) + poly = supported_poly_types[polynomial_type](degree) + fitter = fitting.LevMarLSQFitter() + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + poly_x = fitter(poly, x, y, projection_x) + poly_y = fitter(poly, x, y, projection_y) + transform = models.Mapping((0, 1, 0, 1)) | poly_x & poly_y | projection.inverse | skyrot.inverse + + skyframe = CelestialFrame(reference_frame=fiducial.frame) + detector = Frame2D(name="detector") + pipeline = [(detector, transform), + (skyframe, None) + ] + return WCS(pipeline) diff --git a/readthedocs.yml b/readthedocs.yml new file mode 100644 index 00000000..134a9daf --- /dev/null +++ b/readthedocs.yml @@ -0,0 +1,5 @@ +conda: + file: .rtd-environment.yml +python: + version: 3.5 + setup_py_install: true diff --git a/setup.cfg b/setup.cfg index 13e18196..79e3686e 100644 --- a/setup.cfg +++ b/setup.cfg @@ -9,9 +9,11 @@ show-response = 1 [tool:pytest] minversion = 3.1 +testpaths = gwcs docs norecursedirs = build docs/_build doctest_plus = enabled asdf_schema_root = gwcs/schemas +addopts = --doctest-rst [ah_bootstrap] auto_use = True diff --git a/setup.py b/setup.py index 089be49d..ccfa7297 100755 --- a/setup.py +++ b/setup.py @@ -54,8 +54,7 @@ builtins._ASTROPY_PACKAGE_NAME_ = PACKAGENAME # VERSION should be PEP386 compatible (http://www.python.org/dev/peps/pep-0386) - -VERSION = '0.9.0' +VERSION = '0.10.0' # Indicates if this version is a release version RELEASE = 'dev' not in VERSION @@ -115,6 +114,9 @@ def get_package_data(): entry_points = dict(asdf_extensions='gwcs = gwcs.extension:GWCSExtension') +#import convert_schemas +#convert_schemas.main('gwcs/schemas/stsci.edu/gwcs', 'docs/gwcs/schemas') + setup(name=PACKAGENAME, version=VERSION, description=DESCRIPTION,