diff --git a/CHANGELOG.md b/CHANGELOG.md index b2a255c..521049b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,16 @@ # rc: minor bump +- Fix up `toasty tile-study` to handle FITS files properly (@pkgw, #61). The + input must be in a tangential projection, and only some basic data layouts + within the FITS container are supported. The `--placeholder-thumbnail` option + also must be used. +- Fix an off-by-one error in the computations used by `toasty tile-multi-tan` + (@pkgw, #61) +- Improve some internal APIs for processing studies. + + +# toasty 0.10.0 (2021-09-10) + - Add `toasty check-avm`, which opens up an image file and reports whether it contains AVM (Astronomy Visualization Metadata) tags. This requires that the `pyavm` module is installed (#59, @pkgw). diff --git a/README.md b/README.md index 30a1bd0..9ea1e85 100644 --- a/README.md +++ b/README.md @@ -15,8 +15,9 @@ maps can be viewed in software such as the [AAS] [WorldWide Telescope]. [AAS]: https://aas.org/ [WorldWide Telescope]: http://www.worldwidetelescope.org/ -[toasty] was originally written by [Chris Beaumont] and is currently maintained -as part of the AAS [WorldWide Telescope] project. +[toasty] was originally written by [Chris Beaumont], benefited from +contributions by Clara Brasseur (Space Telescope Science Institute), and is +currently maintained as part of the AAS [WorldWide Telescope] project. [Chris Beaumont]: https://chrisbeaumont.org/ @@ -75,11 +76,13 @@ and [PyPI](https://pypi.org/project/toasty/#history). - [healpy] if using [HEALPix] maps - [numpy] - [pillow] +- [pyavm] if using [Astronomy Visualization Metadata] tags - [pytest] to run the test suite - [PyYAML] - [tqdm] - [wwt_data_formats] >= 0.7 +[Astronomy Visualization Metadata]: https://virtualastronomy.org/avm_metadata.php [astropy]: https://www.astropy.org/ [azure-storage-blob]: https://github.com/Azure/azure-sdk-for-python/tree/master/sdk/storage/azure-storage-blob [cython]: https://cython.org/ @@ -88,6 +91,7 @@ and [PyPI](https://pypi.org/project/toasty/#history). [HEALPix]: https://healpix.jpl.nasa.gov/ [numpy]: https://numpy.org/ [pillow]: https://pillow.readthedocs.io/ +[pyavm]: https://astrofrog.github.io/pyavm/ [pytest]: https://docs.pytest.org/ [PyYAML]: https://github.com/yaml/pyyaml [tqdm]: https://tqdm.github.io/ diff --git a/ci/azure-build-and-test.yml b/ci/azure-build-and-test.yml index 9cfe11b..04e57a8 100644 --- a/ci/azure-build-and-test.yml +++ b/ci/azure-build-and-test.yml @@ -104,6 +104,8 @@ jobs: pyavm \ pytest-cov \ pyyaml \ + reproject \ + shapely \ tqdm pip install wwt_data_formats pip install -e . diff --git a/docs/api/toasty.builder.Builder.rst b/docs/api/toasty.builder.Builder.rst index 04468a5..8c03d89 100644 --- a/docs/api/toasty.builder.Builder.rst +++ b/docs/api/toasty.builder.Builder.rst @@ -23,9 +23,11 @@ Builder ~Builder.cascade ~Builder.create_wtml_folder ~Builder.default_tiled_study_astrometry + ~Builder.execute_study_tiling ~Builder.load_from_wwtl ~Builder.make_placeholder_thumbnail ~Builder.make_thumbnail_from_other + ~Builder.prepare_study_tiling ~Builder.set_name ~Builder.tile_base_as_study ~Builder.toast_base @@ -44,9 +46,11 @@ Builder .. automethod:: cascade .. automethod:: create_wtml_folder .. automethod:: default_tiled_study_astrometry + .. automethod:: execute_study_tiling .. automethod:: load_from_wwtl .. automethod:: make_placeholder_thumbnail .. automethod:: make_thumbnail_from_other + .. automethod:: prepare_study_tiling .. automethod:: set_name .. automethod:: tile_base_as_study .. automethod:: toast_base diff --git a/setup.py b/setup.py index b8374be..3d9e894 100644 --- a/setup.py +++ b/setup.py @@ -87,6 +87,7 @@ def get_long_desc(): 'pytest-cov', ], 'docs': [ + 'astropy', 'astropy-sphinx-theme', 'numpydoc', 'sphinx', diff --git a/toasty/builder.py b/toasty/builder.py index e6958d2..7e477a8 100644 --- a/toasty/builder.py +++ b/toasty/builder.py @@ -71,7 +71,74 @@ def _check_no_wcs_yet(self): raise Exception('order-of-operations error: you must apply WCS after applying tiling settings') + def prepare_study_tiling(self, image): + """ + Set up to tile the specified image as a WWT "study". + + Parameters + ---------- + image : `toasty.image.Image` + The image that will be tiled + + Returns + ------- + tiling : `toasty.study.StudyTiling` + The prepared tiling information + + Remarks + ------- + After calling this method, you should set up the WCS for the tiled + imagery, using :meth:`default_tiled_study_astrometry` as a backstop if + no real information is available. Then use :meth:`execute_study_tiling` + to actually perform the tiling process. + """ + + from .study import StudyTiling + + tiling = StudyTiling(image.width, image.height) + tiling.apply_to_imageset(self.imgset) + return tiling + + + def execute_study_tiling(self, image, tiling, **kwargs): + """ + Tile the specified image as a WWT "study". + + Parameters + ---------- + image : `toasty.image.Image` + The image that will be tiled + tiling : `toasty.study.StudyTiling` + The prepared tiling information + **kwargs + Arguments relayed to :meth:`toasty.study.StudyTiling.tile_image`, + such as ``cli_progress``. + + Returns + ------- + *self* + """ + + tiling.tile_image(image, self.pio, **kwargs) + return self + + def tile_base_as_study(self, image, **kwargs): + """ + Tile an image assuming that it is in the appropriate format for WWT's + "study" framework, namely that it uses a tangential (gnomonic) + projection on the sky. + + Use of this method is somewhat discouraged since it both analyzes and + performs the tiling all at once, which means that you can only correctly + set (and validate) the WCS information *after* doing all the work of + tiling. (Which in turn is because the proper way to apply WCS + information to an imageset depends on the tiling parameters.) It is + generally better to use :meth:`prepare_study_tiling` and + :meth:`execute_study_tiling`, applying the WCS metadata in between, so + that WCS errors can be caught and reported before doing the I/O. + """ + from .study import tile_study_image self._check_no_wcs_yet() diff --git a/toasty/cli.py b/toasty/cli.py index 24be39a..754f6c1 100644 --- a/toasty/cli.py +++ b/toasty/cli.py @@ -418,6 +418,13 @@ def tile_study_impl(settings): pio = PyramidIO(settings.outdir, default_format=img.default_format) builder = Builder(pio) + # First, prepare the tiling, because we can't correctly apply WCS until we + # know the details of the tiling parameters + + tiling = builder.prepare_study_tiling(img) + + # Now deal with the WCS + if settings.avm: # We don't *always* check for AVM because pyavm prints out a lot of junk # and may not be installed. @@ -433,22 +440,26 @@ def tile_study_impl(settings): except Exception: print(f'error: failed to read AVM tags of input `{settings.imgpath}`', file=sys.stderr) raise + + builder.apply_avm_info(avm, img.width, img.height) + elif img.wcs is not None: + # Study images must have negative parity. + img.ensure_negative_parity() + builder.apply_wcs_info(img.wcs, img.width, img.height) else: builder.default_tiled_study_astrometry() - # Do the thumbnail first since for large inputs it can be the memory high-water mark! + # Do the thumbnail before the main image since for large inputs it can be + # the memory high-water mark! if settings.placeholder_thumbnail: builder.make_placeholder_thumbnail() else: builder.make_thumbnail_from_other(img) - builder.tile_base_as_study(img, cli_progress=True) - - if settings.avm: - builder.apply_avm_info(avm, img.width, img.height) - elif img.wcs is not None: - builder.apply_wcs_info(img.wcs, img.width, img.height) + # Finally, actually tile the image + builder.execute_study_tiling(img, tiling, cli_progress=True) + # Final metadata builder.set_name(settings.name) builder.write_index_rel_wtml() diff --git a/toasty/multi_tan.py b/toasty/multi_tan.py index d6c7c2e..ce1f241 100644 --- a/toasty/multi_tan.py +++ b/toasty/multi_tan.py @@ -168,8 +168,8 @@ def compute_global_pixelization(self, builder): desc.sub_tiling = self._tiling.compute_for_subimage( desc.imin, desc.jmin, - desc.imax - desc.imin, - desc.jmax - desc.jmin, + desc.imax + 1 - desc.imin, + desc.jmax + 1 - desc.jmin, ) self._n_todo += desc.sub_tiling.count_populated_positions() diff --git a/toasty/study.py b/toasty/study.py index 1f9c84f..640e485 100644 --- a/toasty/study.py +++ b/toasty/study.py @@ -320,6 +320,12 @@ def tile_image(self, image, pio, cli_progress=False): if image.width != self._width: raise ValueError('width of image to be sampled does not match tiling') + # For tiled FITS, the overall input image and tiling coordinate system + # need to have negative (JPEG-like) parity, but the individal tiles need + # to be saved with positive (FITS-like) parity. The upshot is that we + # have to vertically flip the data layout in our tile buffers. + invert_into_tiles = pio.get_default_vertical_parity_sign() == 1 + # TODO: ideally make_maskable_buffer should be a method # on the Image class which could then avoid having to # manually transfer _format. @@ -328,10 +334,21 @@ def tile_image(self, image, pio, cli_progress=False): with tqdm(total=self.count_populated_positions(), disable=not cli_progress) as progress: for pos, width, height, image_x, image_y, tile_x, tile_y in self.generate_populated_positions(): + if invert_into_tiles: + flip_tile_y1 = 255 - tile_y + flip_tile_y0 = flip_tile_y1 - height + + if flip_tile_y0 == -1: + flip_tile_y0 = None # with a slice, -1 does the wrong thing + + by_idx = slice(flip_tile_y1, flip_tile_y0, -1) + else: + by_idx = slice(tile_y, tile_y + height) + iy_idx = slice(image_y, image_y + height) ix_idx = slice(image_x, image_x + width) - by_idx = slice(tile_y, tile_y + height) bx_idx = slice(tile_x, tile_x + width) + image.fill_into_maskable_buffer(buffer, iy_idx, ix_idx, by_idx, bx_idx) pio.write_image(pos, buffer) progress.update(1) diff --git a/toasty/tests/test_multi_tan.py b/toasty/tests/test_multi_tan.py index a72bc61..9a02954 100644 --- a/toasty/tests/test_multi_tan.py +++ b/toasty/tests/test_multi_tan.py @@ -1,14 +1,15 @@ # -*- mode: python; coding: utf-8 -*- -# Copyright 2019-2020 the AAS WorldWide Telescope project +# Copyright 2019-2021 the AAS WorldWide Telescope project # Licensed under the MIT License. from __future__ import absolute_import, division, print_function import numpy as np -from numpy import testing as nt +import numpy.testing as nt import os.path import pytest import sys +from xml.etree import ElementTree as etree from . import assert_xml_elements_equal, test_path from ..builder import Builder @@ -17,42 +18,48 @@ from .. import multi_tan +try: + from astropy.io import fits + HAS_ASTRO = True +except ImportError: + HAS_ASTRO = False + + +try: + import reproject + HAS_REPROJECT = True +except ImportError: + HAS_REPROJECT = False + + class TestMultiTan(object): WTML = """ - - + + - - CT - CU - TU - DT + + thumb.jpg """ - # Gross workaround for Python 2.7, where the XML serialization is - # apparently slightly different from Python >=3 (for Numpy types only, I - # think?). Fun times. We could avoid this by comparing floats as floats, - # not text, but then we basically have to learn how to deserialize WTML - # with the full semantics of the format. - - if sys.version_info.major == 2: - WTML = (WTML - .replace('Dec="0.74380165289257"', 'Dec="0.743801652893"') - .replace('RA="14.419753086419734"', 'RA="14.4197530864"') - .replace('CenterX="216.296296296296"', 'CenterX="216.296296296"') - .replace('CenterY="0.74380165289257"', 'CenterY="0.743801652893"') - ) + # Gross workaround for platform differences in the XML output. + + if sys.platform == 'darwin': + WTML = WTML.replace('Dec="0.7438249862258411"', 'Dec="0.743824986225841"') # Back to the non-gross stuff. @@ -67,6 +74,7 @@ def teardown_method(self, method): def work_path(self, *pieces): return os.path.join(self.work_dir, *pieces) + def test_basic(self): coll = collection.SimpleFitsCollection([test_path('wcs512.fits.gz')]) @@ -80,12 +88,56 @@ def test_basic(self): proc.compute_global_pixelization(builder) proc.tile(pio) + BARY_SLICES = [ + (slice(0, 128), slice(0, 128)), + (slice(0, 128), slice(128, None)), + (slice(128, None), slice(0, 128)), + (slice(128, None), slice(128, None)), + ] + + def maybe_test_barycenter(self, path, bary_expected): + """ + Check the barycenters of four 128x128 quadrants of a tile file. The idea + here is that if we introduce a problem with vertical flips in tiled FITS + processing, we'll detect it here. + """ + + if not HAS_ASTRO: + return + + with fits.open(path) as hdul: + data = hdul[0].data + + data[~np.isfinite(data)] = 0.0 + bary_observed = [] + + for islice in self.BARY_SLICES: + idata = data[islice] + yidx, xidx = np.indices((128, 128)) + xbary = (idata * xidx).sum() / idata.sum() + ybary = (idata * yidx).sum() / idata.sum() + bary_observed.append((xbary, ybary)) + + nt.assert_array_almost_equal(bary_observed, bary_expected, decimal=5) + + WCS512_BARYDATA = [ + (63.44949378800272, 64.40535387506924), + (63.24744175084746, 63.67473452789256), + (65.22950207855361, 63.35629429568745), + (62.027396724898814, 62.815937534782144) + ] + def test_basic_cli(self): """ Test the CLI interface. We don't go out of our way to validate the computations in detail -- that's for the unit tests that probe the module directly. """ + expected = etree.fromstring( + self.WTML + .replace('Thumbnail="thumb.jpg"', '') + .replace('thumb.jpg', '') + ) args = [ 'tile-multi-tan', @@ -94,3 +146,77 @@ def test_basic_cli(self): test_path('wcs512.fits.gz') ] cli.entrypoint(args) + + with open(self.work_path('basic_cli', 'index_rel.wtml'), 'rt', encoding='utf8') as f: + observed = etree.fromstring(f.read()) + + assert_xml_elements_equal(observed, expected) + + args = [ + 'cascade', + '--start', '1', + self.work_path('basic_cli'), + ] + cli.entrypoint(args) + + self.maybe_test_barycenter(self.work_path('basic_cli', '0', '0', '0_0.fits'), self.WCS512_BARYDATA) + + + def test_study_cli(self): + """ + Test tile-study on FITS. This should properly go in test_study.py, but + this file is the one that has the reference WTML information. + """ + expected = etree.fromstring(self.WTML) + + args = [ + 'tile-study', + '--placeholder-thumbnail', + '--outdir', self.work_path('study_cli'), + test_path('wcs512.fits.gz') + ] + cli.entrypoint(args) + + with open(self.work_path('study_cli', 'index_rel.wtml'), 'rt', encoding='utf8') as f: + observed = etree.fromstring(f.read()) + + assert_xml_elements_equal(observed, expected) + + args = [ + 'cascade', + '--start', '1', + self.work_path('study_cli'), + ] + cli.entrypoint(args) + + self.maybe_test_barycenter(self.work_path('study_cli', '0', '0', '0_0.fits'), self.WCS512_BARYDATA) + + + @pytest.mark.skipif('not HAS_REPROJECT') + def test_as_multi_wcs(self): + """ + Once again, this doesn't super belong here, but this is where we have + the reference data. We don't compare the WTML contents here since the + reprojection isn't going to preserve the WCS in detail. + """ + from .. import builder, collection, multi_wcs, pyramid + + reproject_function = reproject.reproject_interp + outdir = self.work_path('as_multi_wcs') + + pio = pyramid.PyramidIO(outdir, default_format='fits') + bld = builder.Builder(pio) + coll = collection.SimpleFitsCollection([test_path('wcs512.fits.gz')], hdu_index=0) + proc = multi_wcs.MultiWcsProcessor(coll) + proc.compute_global_pixelization(bld) + proc.tile(pio, reproject_function, cli_progress=False, parallel=1) + bld.write_index_rel_wtml() + + args = [ + 'cascade', + '--start', '1', + self.work_path('as_multi_wcs'), + ] + cli.entrypoint(args) + + self.maybe_test_barycenter(self.work_path('as_multi_wcs', '0', '0', '0_0.fits'), self.WCS512_BARYDATA)