From 7d8b5c165b416ab657144d9c6774b914cd9a5ea5 Mon Sep 17 00:00:00 2001 From: Peter Williams Date: Sat, 12 Feb 2022 11:19:56 -0500 Subject: [PATCH 01/11] toasty check-avm: add some spatial info diagnostics --- docs/cli/check-avm.rst | 14 +++++----- toasty/cli.py | 58 ++++++++++++++++++++++++++++++++++++---- toasty/tests/test_avm.py | 29 ++++++++++++-------- 3 files changed, 79 insertions(+), 22 deletions(-) diff --git a/docs/cli/check-avm.rst b/docs/cli/check-avm.rst index d98ecac..28ec76e 100644 --- a/docs/cli/check-avm.rst +++ b/docs/cli/check-avm.rst @@ -28,14 +28,16 @@ successful. .. _WCS: https://fits.gsfc.nasa.gov/fits_wcs.html -By default, only the presence of AVM data is reported. If the ``--print`` or -``-p`` option is given, the AVM data will be dumped out. +By default, the presence of AVM data is reported along with some information +about the agreement between the AVM and the image shape. If the ``--print`` or +``-p`` option is given, the actual AVM data will be dumped out. If the ``--exitcode`` option is given, the exit code of the command line tool -will reflect whether AVM data were found. Zero (success) means that data were -found, while nonzero (failure) indicates that no data were present, the data -were corrupt, or that some other failure occurred. Without this option, the exit -code will be zero if the image can be opened but contains no AVM data. +will reflect whether spatial AVM data were found. Zero (success) means that +spatial data were found, while nonzero (failure) indicates that no data were +present, the data were corrupt, or that some other failure occurred. Without +this option, the exit code will be zero if the image can be opened but contains +no AVM data. Notes diff --git a/toasty/cli.py b/toasty/cli.py index d78d857..cb3f1c2 100644 --- a/toasty/cli.py +++ b/toasty/cli.py @@ -105,6 +105,8 @@ def check_avm_getparser(parser): def check_avm_impl(settings): + from .image import ImageLoader + try: from pyavm import AVM, exceptions @@ -129,15 +131,61 @@ def check_avm_impl(settings): f"{settings.image}: AVM data are present but malformed: {e} ({e.__class__.__name__})" ) sys.exit(1) - except IOError as e: - die(f"error probing AVM tags of `{settings.image}`: {e}") + except Exception as e: + die( + f"error probing AVM tags of `{settings.image}`: {e} ({e.__class__.__name__})" + ) + + status_exitcode = 0 + + if settings.print_avm: + print_remark = "" + else: + print_remark = " (use `--print` to display them)" + + print(f"{settings.image}: AVM tags are present{print_remark}") + + if repr(avm.Spatial) == "": + print(f"{settings.image}: however, no spatial information") + status_exitcode = 1 + elif avm.Spatial.ReferenceDimension is None: + print( + f"{settings.image}: including spatial information, but no reference dimension" + ) + else: + img = ImageLoader().load_path(settings.image) + this_shape = (int(img.width), int(img.height)) + ref_shape = ( + int(avm.Spatial.ReferenceDimension[0]), + int(avm.Spatial.ReferenceDimension[1]), + ) + + if this_shape == ref_shape: + print( + f"{settings.image}: including spatial information for this image's dimensions" + ) + else: + print( + f"{settings.image}: including spatial information at a different scale" + ) + print( + f"{settings.image}: AVM is for {ref_shape[0]} x {ref_shape[1]}, image is {this_shape[0]} x {this_shape[1]}" + ) + this_aspect = this_shape[0] / this_shape[1] + ref_aspect = ref_shape[0] / ref_shape[1] + + # The threshold here approximates the one used by pyavm + if abs(this_aspect - ref_aspect) / (this_aspect + ref_aspect) >= 0.005: + print( + f"{settings.image}: WARNING: these aspect ratios are not compatible" + ) + status_exitcode = 1 if settings.print_avm: - print(f"{settings.image}: AVM tags are present") print() print(avm) - else: - print(f"{settings.image}: AVM tags are present (use `--print` to display them)") + + sys.exit(status_exitcode if settings.exitcode else 0) # "make_thumbnail" subcommand diff --git a/toasty/tests/test_avm.py b/toasty/tests/test_avm.py index a6b2027..72c5f2c 100644 --- a/toasty/tests/test_avm.py +++ b/toasty/tests/test_avm.py @@ -12,22 +12,24 @@ class TestAvm(object): def setup_method(self, method): from tempfile import mkdtemp + self.work_dir = mkdtemp() def teardown_method(self, method): from shutil import rmtree + rmtree(self.work_dir) def work_path(self, *pieces): return os.path.join(self.work_dir, *pieces) - @pytest.mark.skipif('not HAS_AVM') + @pytest.mark.skipif("not HAS_AVM") def test_check_cli_bad(self): args = [ - 'check-avm', - '--print', - '--exitcode', - test_path('badavm-type.png'), + "check-avm", + "--print", + "--exitcode", + test_path("badavm-type.png"), ] try: @@ -35,14 +37,19 @@ def test_check_cli_bad(self): except SystemExit as e: assert e.code == 1 else: - assert False, 'no error exit on bad AVM' + assert False, "no error exit on bad AVM" - @pytest.mark.skipif('not HAS_AVM') + @pytest.mark.skipif("not HAS_AVM") def test_check_cli_good(self): args = [ - 'check-avm', - '--print', - test_path('geminiann11015a.jpg'), + "check-avm", + "--print", + test_path("geminiann11015a.jpg"), ] - cli.entrypoint(args) + try: + cli.entrypoint(args) + except SystemExit as e: + assert e.code == 0 + else: + assert False, "error exit on good AVM" From 6e88c21d8c7e154a1de67c29aae24d6642fd7a89 Mon Sep 17 00:00:00 2001 From: Peter Williams Date: Sat, 12 Feb 2022 11:20:37 -0500 Subject: [PATCH 02/11] toasty/toast.py: fix a new doc warning --- toasty/toast.py | 189 +++++++++++++++++++++++++++++++----------------- 1 file changed, 122 insertions(+), 67 deletions(-) diff --git a/toasty/toast.py b/toasty/toast.py index a9bedbb..7d8d7a7 100644 --- a/toasty/toast.py +++ b/toasty/toast.py @@ -23,7 +23,7 @@ """ from __future__ import absolute_import, division, print_function -__all__ = ''' +__all__ = """ count_tiles_matching_filter create_single_tile generate_tiles @@ -36,7 +36,7 @@ toast_tile_area toast_tile_for_point toast_tile_get_coords -'''.split() +""".split() from collections import namedtuple from enum import Enum @@ -73,7 +73,12 @@ def _spherical_triangle_area(lat1, lon1, lat2, lon2, lat3, lon3): a = _arclength(lat2, lon2, lat3, lon3) b = _arclength(lat3, lon3, lat1, lon1) s = 0.5 * (a + b + c) - tane4 = np.sqrt(np.tan(0.5 * s) * np.tan(0.5 * (s - a)) * np.tan(0.5 * (s - b)) * np.tan(0.5 * (s - c))) + tane4 = np.sqrt( + np.tan(0.5 * s) + * np.tan(0.5 * (s - a)) + * np.tan(0.5 * (s - b)) + * np.tan(0.5 * (s - c)) + ) e = 4 * np.arctan(tane4) return e @@ -83,24 +88,26 @@ class ToastCoordinateSystem(Enum): Different TOAST coordinate systems that are in use. """ - ASTRONOMICAL = 'astronomical' + ASTRONOMICAL = "astronomical" """The default TOAST coordinate system, where the ``lat = lon = 0`` point lies at the middle right edge of the TOAST projection square.""" - PLANETARY = 'planetary' + PLANETARY = "planetary" """The planetary TOAST coordinate system. This is rotated 180 degrees in longitude from the astronomical system, such that the ``lat = lon = 0`` point lies at the middle left edge of the TOAST projection square.""" -Tile = namedtuple('Tile', 'pos corners increasing') +Tile = namedtuple("Tile", "pos corners increasing") -_level1_astronomical_lonlats = np.radians([ - [(0, -90), (90, 0), (0, 90), (180, 0)], - [(90, 0), (0, -90), (0, 0), (0, 90)], - [(180, 0), (0, 90), (270, 0), (0, -90)], - [(0, 90), (0, 0), (0, -90), (270, 0)], -]) +_level1_astronomical_lonlats = np.radians( + [ + [(0, -90), (90, 0), (0, 90), (180, 0)], + [(90, 0), (0, -90), (0, 0), (0, 90)], + [(180, 0), (0, 90), (270, 0), (0, -90)], + [(0, 90), (0, 0), (0, -90), (270, 0)], + ] +) _level1_astronomical_lonlats.flags.writeable = False @@ -109,7 +116,7 @@ def _create_level1_tiles(coordsys): if coordsys == ToastCoordinateSystem.PLANETARY: lonlats = lonlats.copy() - lonlats[...,0] = (lonlats[...,0] + np.pi) % TWOPI + lonlats[..., 0] = (lonlats[..., 0] + np.pi) % TWOPI return [ Tile(Pos(n=1, x=0, y=0), lonlats[0], True), @@ -155,11 +162,13 @@ def _equ_to_xyz(lat, lon): """ clat = np.cos(lat) - return np.array([ - np.cos(lon) * clat, - np.sin(lat), - np.sin(lon) * clat, - ]) + return np.array( + [ + np.cos(lon) * clat, + np.sin(lat), + np.sin(lon) * clat, + ] + ) def _left_of_half_space_score(point_a, point_b, test_point): @@ -267,7 +276,7 @@ def toast_tile_for_point(depth, lat, lon, coordsys=ToastCoordinateSystem.ASTRONO return Tile(Pos(n=0, x=0, y=0), (None, None, None, None), False) for tile in _create_level1_tiles(coordsys): - if _toast_tile_containment_score(tile, lat, lon) == 0.: + if _toast_tile_containment_score(tile, lat, lon) == 0.0: break while tile.pos.n < depth: @@ -283,7 +292,7 @@ def toast_tile_for_point(depth, lat, lon, coordsys=ToastCoordinateSystem.ASTRONO for child in _div4(tile): score = _toast_tile_containment_score(child, lat, lon) - if score == 0.: + if score == 0.0: tile = child break @@ -355,7 +364,7 @@ def toast_pixel_for_point(depth, lat, lon, coordsys=ToastCoordinateSystem.ASTRON # that is closest to the input position. lons, lats = toast_tile_get_coords(tile) - dist2 = (lons - lon)**2 + (lats - lat)**2 + dist2 = (lons - lon) ** 2 + (lats - lat) ** 2 min_y, min_x = np.unravel_index(np.argmin(dist2), (256, 256)) # Now, identify a postage stamp around that best-fit pixel and fit a biquadratic @@ -367,21 +376,23 @@ def toast_pixel_for_point(depth, lat, lon, coordsys=ToastCoordinateSystem.ASTRON x1 = min(min_x + halfsize + 1, 256) y1 = min(min_y + halfsize + 1, 256) - dist2_stamp = dist2[y0:y1,x0:x1] - lons_stamp = lons[y0:y1,x0:x1] - lats_stamp = lats[y0:y1,x0:x1] + dist2_stamp = dist2[y0:y1, x0:x1] + lons_stamp = lons[y0:y1, x0:x1] + lats_stamp = lats[y0:y1, x0:x1] flat_lons = lons_stamp.flatten() flat_lats = lats_stamp.flatten() - A = np.array([ - flat_lons * 0 + 1, - flat_lons, - flat_lats, - flat_lons**2, - flat_lons * flat_lats, - flat_lats**2, - ]).T + A = np.array( + [ + flat_lons * 0 + 1, + flat_lons, + flat_lats, + flat_lons ** 2, + flat_lons * flat_lats, + flat_lats ** 2, + ] + ).T ygrid, xgrid = np.indices(dist2_stamp.shape) x_coeff, _r, _rank, _s = np.linalg.lstsq(A, xgrid.flatten(), rcond=None) @@ -389,14 +400,16 @@ def toast_pixel_for_point(depth, lat, lon, coordsys=ToastCoordinateSystem.ASTRON # Evaluate the polynomial to get the refined pixel coordinates. - pt = np.array([ - 1, - lon, - lat, - lon**2, - lon * lat, - lat**2, - ]) + pt = np.array( + [ + 1, + lon, + lat, + lon ** 2, + lon * lat, + lat ** 2, + ] + ) x = np.dot(x_coeff, pt) y = np.dot(y_coeff, pt) return tile, x0 + x, y0 + y @@ -451,9 +464,9 @@ def _div4(tile): y *= 2 return [ - Tile(Pos(n=n, x=x, y=y ), (ul, to, ce, le), increasing), - Tile(Pos(n=n, x=x + 1, y=y ), (to, ur, ri, ce), increasing), - Tile(Pos(n=n, x=x, y=y + 1), (le, ce, bo, ll), increasing), + Tile(Pos(n=n, x=x, y=y), (ul, to, ce, le), increasing), + Tile(Pos(n=n, x=x + 1, y=y), (to, ur, ri, ce), increasing), + Tile(Pos(n=n, x=x, y=y + 1), (le, ce, bo, ll), increasing), Tile(Pos(n=n, x=x + 1, y=y + 1), (ce, ri, lr, bo), increasing), ] @@ -482,7 +495,7 @@ def create_single_tile(pos, coordsys=ToastCoordinateSystem.ASTRONOMICAL): """ if pos.n == 0: - raise ValueError('cannot create a Tile for the n=0 tile') + raise ValueError("cannot create a Tile for the n=0 tile") children = _create_level1_tiles(coordsys) cur_n = 0 @@ -499,7 +512,9 @@ def create_single_tile(pos, coordsys=ToastCoordinateSystem.ASTRONOMICAL): children = _div4(tile) -def generate_tiles(depth, bottom_only=True, coordsys=ToastCoordinateSystem.ASTRONOMICAL): +def generate_tiles( + depth, bottom_only=True, coordsys=ToastCoordinateSystem.ASTRONOMICAL +): """Generate a pyramid of TOAST tiles in deepest-first order. Parameters @@ -520,10 +535,14 @@ def generate_tiles(depth, bottom_only=True, coordsys=ToastCoordinateSystem.ASTRO The ``n = 0`` depth is not included. """ - return generate_tiles_filtered(depth, lambda t: True, bottom_only, coordsys=coordsys) + return generate_tiles_filtered( + depth, lambda t: True, bottom_only, coordsys=coordsys + ) -def generate_tiles_filtered(depth, filter, bottom_only=True, coordsys=ToastCoordinateSystem.ASTRONOMICAL): +def generate_tiles_filtered( + depth, filter, bottom_only=True, coordsys=ToastCoordinateSystem.ASTRONOMICAL +): """Generate a pyramid of TOAST tiles in deepest-first order, filtering out subtrees. Parameters @@ -552,7 +571,9 @@ def generate_tiles_filtered(depth, filter, bottom_only=True, coordsys=ToastCoord yield item -def count_tiles_matching_filter(depth, filter, bottom_only=True, coordsys=ToastCoordinateSystem.ASTRONOMICAL): +def count_tiles_matching_filter( + depth, filter, bottom_only=True, coordsys=ToastCoordinateSystem.ASTRONOMICAL +): """ Count the number of tiles matching a filter. @@ -570,7 +591,7 @@ def count_tiles_matching_filter(depth, filter, bottom_only=True, coordsys=ToastC :attr:`ToastCoordinateSystem.ASTRONOMICAL`. Returns - ------ + ------- The number of tiles matching the filter. Even if ``bottom_only`` is false, the ``n = 0`` tile is not counted. @@ -582,7 +603,9 @@ def count_tiles_matching_filter(depth, filter, bottom_only=True, coordsys=ToastC # With a generic filter function, brute force is our only option: n = 0 - for _tile in generate_tiles_filtered(depth, filter, bottom_only=bottom_only, coordsys=coordsys): + for _tile in generate_tiles_filtered( + depth, filter, bottom_only=bottom_only, coordsys=coordsys + ): n += 1 return n @@ -627,10 +650,13 @@ def sample_layer( """ from .par_util import resolve_parallelism + parallel = resolve_parallelism(parallel) if parallel > 1: - _sample_layer_parallel(pio, format, sampler, depth, coordsys, cli_progress, parallel) + _sample_layer_parallel( + pio, format, sampler, depth, coordsys, cli_progress, parallel + ) else: _sample_layer_serial(pio, format, sampler, depth, coordsys, cli_progress) @@ -647,15 +673,19 @@ def _sample_layer_serial(pio, format, sampler, depth, coordsys, cli_progress): print() -def _sample_layer_parallel(pio, format, sampler, depth, coordsys, cli_progress, parallel): +def _sample_layer_parallel( + pio, format, sampler, depth, coordsys, cli_progress, parallel +): import multiprocessing as mp done_event = mp.Event() - queue = mp.Queue(maxsize = 2 * parallel) + queue = mp.Queue(maxsize=2 * parallel) workers = [] for _ in range(parallel): - w = mp.Process(target=_mp_sample_worker, args=(queue, done_event, pio, sampler, format)) + w = mp.Process( + target=_mp_sample_worker, args=(queue, done_event, pio, sampler, format) + ) w.daemon = True w.start() workers.append(w) @@ -738,25 +768,38 @@ def sample_layer_filtered( """ from .par_util import resolve_parallelism + parallel = resolve_parallelism(parallel) if parallel > 1: - _sample_filtered_parallel(pio, tile_filter, sampler, depth, coordsys, cli_progress, parallel) + _sample_filtered_parallel( + pio, tile_filter, sampler, depth, coordsys, cli_progress, parallel + ) else: - _sample_filtered_serial(pio, tile_filter, sampler, depth, coordsys, cli_progress) + _sample_filtered_serial( + pio, tile_filter, sampler, depth, coordsys, cli_progress + ) def _sample_filtered_serial(pio, tile_filter, sampler, depth, coordsys, cli_progress): - n_todo = count_tiles_matching_filter(depth, tile_filter, bottom_only=True, coordsys=coordsys) + n_todo = count_tiles_matching_filter( + depth, tile_filter, bottom_only=True, coordsys=coordsys + ) with tqdm(total=n_todo, disable=not cli_progress) as progress: - for tile in generate_tiles_filtered(depth, tile_filter, bottom_only=True, coordsys=coordsys): + for tile in generate_tiles_filtered( + depth, tile_filter, bottom_only=True, coordsys=coordsys + ): lon, lat = toast_tile_get_coords(tile) sampled_data = sampler(lon, lat) img = Image.from_array(sampled_data) - with pio.update_image(tile.pos, masked_mode=img.mode, default='masked') as basis: - img.update_into_maskable_buffer(basis, slice(None), slice(None), slice(None), slice(None)) + with pio.update_image( + tile.pos, masked_mode=img.mode, default="masked" + ) as basis: + img.update_into_maskable_buffer( + basis, slice(None), slice(None), slice(None), slice(None) + ) progress.update(1) @@ -767,17 +810,23 @@ def _sample_filtered_serial(pio, tile_filter, sampler, depth, coordsys, cli_prog # chunks in parallel. -def _sample_filtered_parallel(pio, tile_filter, sampler, depth, coordsys, cli_progress, parallel): +def _sample_filtered_parallel( + pio, tile_filter, sampler, depth, coordsys, cli_progress, parallel +): import multiprocessing as mp - n_todo = count_tiles_matching_filter(depth, tile_filter, bottom_only=True, coordsys=coordsys) + n_todo = count_tiles_matching_filter( + depth, tile_filter, bottom_only=True, coordsys=coordsys + ) done_event = mp.Event() - queue = mp.Queue(maxsize = 2 * parallel) + queue = mp.Queue(maxsize=2 * parallel) workers = [] for _ in range(parallel): - w = mp.Process(target=_mp_sample_filtered, args=(queue, done_event, pio, sampler)) + w = mp.Process( + target=_mp_sample_filtered, args=(queue, done_event, pio, sampler) + ) w.daemon = True w.start() workers.append(w) @@ -785,7 +834,9 @@ def _sample_filtered_parallel(pio, tile_filter, sampler, depth, coordsys, cli_pr # Here we go: with tqdm(total=n_todo, disable=not cli_progress) as progress: - for tile in generate_tiles_filtered(depth, tile_filter, bottom_only=True, coordsys=coordsys): + for tile in generate_tiles_filtered( + depth, tile_filter, bottom_only=True, coordsys=coordsys + ): queue.put(tile) progress.update(1) @@ -828,5 +879,9 @@ def _mp_sample_filtered(queue, done_event, pio, sampler): sampled_data = sampler(lon, lat) img = Image.from_array(sampled_data) - with pio.update_image(tile.pos, masked_mode=img.mode, default='masked') as basis: - img.update_into_maskable_buffer(basis, slice(None), slice(None), slice(None), slice(None)) + with pio.update_image( + tile.pos, masked_mode=img.mode, default="masked" + ) as basis: + img.update_into_maskable_buffer( + basis, slice(None), slice(None), slice(None), slice(None) + ) From 0d112b207652cd9b9e11f3519c311855f694fc79 Mon Sep 17 00:00:00 2001 From: Peter Williams Date: Sat, 12 Feb 2022 12:27:15 -0500 Subject: [PATCH 03/11] toasty/tests: update for new wwt_data_formats constellation support --- setup.py | 4 ++-- toasty/tests/test_multi_tan.py | 1 + toasty/tests/test_study.py | 2 ++ 3 files changed, 5 insertions(+), 2 deletions(-) diff --git a/setup.py b/setup.py index 00671e6..92eae4f 100644 --- a/setup.py +++ b/setup.py @@ -1,5 +1,5 @@ # -*- mode: python; coding: utf-8 -*- -# Copyright 2013-2020 Chris Beaumont and the AAS WorldWide Telescope team +# Copyright 2013-2022 Chris Beaumont and the AAS WorldWide Telescope team # Licensed under the MIT License from __future__ import absolute_import, division, print_function @@ -78,7 +78,7 @@ def get_long_desc(): "reproject", "shapely", "tqdm>=4.0", - "wwt_data_formats>=0.12", + "wwt_data_formats>=0.13", ], extras_require={ "test": [ diff --git a/toasty/tests/test_multi_tan.py b/toasty/tests/test_multi_tan.py index 8b3711d..c933c0f 100644 --- a/toasty/tests/test_multi_tan.py +++ b/toasty/tests/test_multi_tan.py @@ -40,6 +40,7 @@ class TestMultiTan(object): Date: Tue, 22 Mar 2022 10:08:18 -0400 Subject: [PATCH 04/11] toasty/samplers.py: avoid HEALPix HDUs without data --- toasty/samplers.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/toasty/samplers.py b/toasty/samplers.py index acb1e50..1d3e91d 100644 --- a/toasty/samplers.py +++ b/toasty/samplers.py @@ -98,8 +98,13 @@ def _find_healpix_extension_index(pth): """ for i, hdu in enumerate(pth): - if hdu.header.get("PIXTYPE") == "HEALPIX": - return i + if hdu.data is None: + continue + + if hdu.header.get("PIXTYPE") != "HEALPIX": + continue + + return i else: raise IndexError("No HEALPIX extensions found in %s" % pth.filename()) From eb0c95029341cb81ef29935f89595a33426d9bda Mon Sep 17 00:00:00 2001 From: Peter Williams Date: Tue, 22 Mar 2022 10:17:18 -0400 Subject: [PATCH 05/11] Add --force-galactic mode for tile-healpix --- toasty/cli.py | 9 ++++++++- toasty/samplers.py | 12 ++++++++++-- 2 files changed, 18 insertions(+), 3 deletions(-) diff --git a/toasty/cli.py b/toasty/cli.py index cb3f1c2..d404fa5 100644 --- a/toasty/cli.py +++ b/toasty/cli.py @@ -366,6 +366,11 @@ def tile_allsky_impl(settings): def tile_healpix_getparser(parser): + parser.add_argument( + "--galactic", + action="store_true", + help="Force use of Galactic coordinate system, regardless of headers", + ) parser.add_argument( "--outdir", metavar="PATH", @@ -391,7 +396,9 @@ def tile_healpix_impl(settings): from .samplers import healpix_fits_file_sampler pio = PyramidIO(settings.outdir, default_format="npy") - sampler = healpix_fits_file_sampler(settings.fitspath) + sampler = healpix_fits_file_sampler( + settings.fitspath, force_galactic=settings.galactic + ) builder = Builder(pio) builder.toast_base(sampler, settings.depth) builder.write_index_rel_wtml() diff --git a/toasty/samplers.py b/toasty/samplers.py index 1d3e91d..0131fc6 100644 --- a/toasty/samplers.py +++ b/toasty/samplers.py @@ -78,7 +78,7 @@ def healpix_sampler(data, nest=False, coord="C", interpolation="nearest"): def vec2pix(l, b): if galactic: f = ICRS(l * u.rad, b * u.rad) - g = f.transform_to(Galactic) + g = f.transform_to(Galactic()) l, b = g.l.rad, g.b.rad theta = np.pi / 2 - b @@ -109,7 +109,9 @@ def _find_healpix_extension_index(pth): raise IndexError("No HEALPIX extensions found in %s" % pth.filename()) -def healpix_fits_file_sampler(path, extension=None, interpolation="nearest"): +def healpix_fits_file_sampler( + path, extension=None, interpolation="nearest", force_galactic=False +): """Create a sampler for HEALPix data read from a FITS file. Parameters @@ -123,6 +125,9 @@ def healpix_fits_file_sampler(path, extension=None, interpolation="nearest"): What interpolation scheme to use. WARNING: bilinear uses healpy's get_interp_val, which seems prone to segfaults + force_galactic : bool + If true, force use of Galactic coordinate system, regardless of + what the headers say. Returns ------- @@ -148,6 +153,9 @@ def healpix_fits_file_sampler(path, extension=None, interpolation="nearest"): nest = hdr.get("ORDERING") == "NESTED" coord = hdr.get("COORDSYS", "C") + if force_galactic: + coord = "G" + return healpix_sampler(data, nest, coord, interpolation) From 699608e03244304fc6a83222e3e0a8ca7bafd2e8 Mon Sep 17 00:00:00 2001 From: Peter Williams Date: Tue, 22 Mar 2022 10:35:19 -0400 Subject: [PATCH 06/11] toasty tile-healpix: support tiled FITS The default output format was npy, but we change it to FITS. As per usual we have to add some extra plumbing to flip image buffers vertically when saving to FITS. Also tidy up the CLI, including adding support for controlling the parallelization. --- toasty/cli.py | 16 ++++++++++++++-- toasty/toast.py | 28 ++++++++++++++++++++++++++++ 2 files changed, 42 insertions(+), 2 deletions(-) diff --git a/toasty/cli.py b/toasty/cli.py index d404fa5..cb8bdcd 100644 --- a/toasty/cli.py +++ b/toasty/cli.py @@ -377,6 +377,13 @@ def tile_healpix_getparser(parser): default=".", help="The root directory of the output tile pyramid", ) + parser.add_argument( + "--parallelism", + "-j", + metavar="COUNT", + type=int, + help="The parallelization level (default: use all CPUs if OS supports; specify `1` to force serial processing)", + ) parser.add_argument( "fitspath", metavar="PATH", @@ -395,12 +402,17 @@ def tile_healpix_impl(settings): from .pyramid import PyramidIO from .samplers import healpix_fits_file_sampler - pio = PyramidIO(settings.outdir, default_format="npy") + pio = PyramidIO(settings.outdir, default_format="fits") sampler = healpix_fits_file_sampler( settings.fitspath, force_galactic=settings.galactic ) builder = Builder(pio) - builder.toast_base(sampler, settings.depth) + builder.toast_base( + sampler, + settings.depth, + parallel=settings.parallelism, + cli_progress=True, + ) builder.write_index_rel_wtml() print( diff --git a/toasty/toast.py b/toasty/toast.py index 7d8d7a7..1fe2d2f 100644 --- a/toasty/toast.py +++ b/toasty/toast.py @@ -662,10 +662,17 @@ def sample_layer( def _sample_layer_serial(pio, format, sampler, depth, coordsys, cli_progress): + # The usual vertical flip that we may need if tiling into FITS: + invert_into_tiles = pio.get_default_vertical_parity_sign() == 1 + with tqdm(total=tiles_at_depth(depth), disable=not cli_progress) as progress: for tile in generate_tiles(depth, bottom_only=True, coordsys=coordsys): lon, lat = toast_tile_get_coords(tile) sampled_data = sampler(lon, lat) + + if invert_into_tiles: + sampled_data = sampled_data[::-1] + pio.write_image(tile.pos, Image.from_array(sampled_data), format=format) progress.update(1) @@ -716,6 +723,9 @@ def _mp_sample_worker(queue, done_event, pio, sampler, format): """ from queue import Empty + # The usual vertical flip that we may need if tiling into FITS: + invert_into_tiles = pio.get_default_vertical_parity_sign() == 1 + while True: try: tile = queue.get(True, timeout=1) @@ -726,6 +736,10 @@ def _mp_sample_worker(queue, done_event, pio, sampler, format): lon, lat = toast_tile_get_coords(tile) sampled_data = sampler(lon, lat) + + if invert_into_tiles: + sampled_data = sampled_data[::-1] + pio.write_image(tile.pos, Image.from_array(sampled_data), format=format) @@ -786,12 +800,19 @@ def _sample_filtered_serial(pio, tile_filter, sampler, depth, coordsys, cli_prog depth, tile_filter, bottom_only=True, coordsys=coordsys ) + # The usual vertical flip that we may need if tiling into FITS: + invert_into_tiles = pio.get_default_vertical_parity_sign() == 1 + with tqdm(total=n_todo, disable=not cli_progress) as progress: for tile in generate_tiles_filtered( depth, tile_filter, bottom_only=True, coordsys=coordsys ): lon, lat = toast_tile_get_coords(tile) sampled_data = sampler(lon, lat) + + if invert_into_tiles: + sampled_data = sampled_data[::-1] + img = Image.from_array(sampled_data) with pio.update_image( @@ -867,6 +888,9 @@ def _mp_sample_filtered(queue, done_event, pio, sampler): from queue import Empty + # The usual vertical flip that we may need if tiling into FITS: + invert_into_tiles = pio.get_default_vertical_parity_sign() == 1 + while True: try: tile = queue.get(True, timeout=1) @@ -877,6 +901,10 @@ def _mp_sample_filtered(queue, done_event, pio, sampler): lon, lat = toast_tile_get_coords(tile) sampled_data = sampler(lon, lat) + + if invert_into_tiles: + sampled_data = sampled_data[::-1] + img = Image.from_array(sampled_data) with pio.update_image( From 1d6f8e949a5ac6488f305b1adb46460eca1c2c88 Mon Sep 17 00:00:00 2001 From: Peter Williams Date: Fri, 1 Apr 2022 10:13:38 -0400 Subject: [PATCH 07/11] docs: document tile-healpix --- docs/cli.rst | 1 + docs/cli/tile-allsky.rst | 6 +++ docs/cli/tile-healpix.rst | 82 +++++++++++++++++++++++++++++++++++++++ docs/toasting.rst | 12 +++--- 4 files changed, 95 insertions(+), 6 deletions(-) create mode 100644 docs/cli/tile-healpix.rst diff --git a/docs/cli.rst b/docs/cli.rst index 1e31118..92162cc 100644 --- a/docs/cli.rst +++ b/docs/cli.rst @@ -16,6 +16,7 @@ CLI Reference cli/pipeline-publish cli/pipeline-refresh cli/tile-allsky + cli/tile-healpix cli/tile-study cli/tile-wwtl cli/transform-fx3-to-rgb diff --git a/docs/cli/tile-allsky.rst b/docs/cli/tile-allsky.rst index 7bf4c11..e370778 100644 --- a/docs/cli/tile-allsky.rst +++ b/docs/cli/tile-allsky.rst @@ -105,3 +105,9 @@ because ``fork()``-based multiprocessing is required. MacOS should support this, but there is currently (as of Python 3.8) `a bug`_ preventing that. .. _a bug: https://bugs.python.org/issue33725 + + +See Also +======== + +- :ref:`cli-tile-healpix` \ No newline at end of file diff --git a/docs/cli/tile-healpix.rst b/docs/cli/tile-healpix.rst new file mode 100644 index 0000000..d28965d --- /dev/null +++ b/docs/cli/tile-healpix.rst @@ -0,0 +1,82 @@ +.. _cli-tile-healpix: + +======================= +``toasty tile-healpix`` +======================= + +The ``tile-healpix`` command transforms all-sky FITS images in the HEALPix_ +format into WWT’s FITS TOAST format. + +.. _HEALPix: https://healpix.jpl.nasa.gov/ + +Usage +===== + +.. code-block:: shell + + toasty tile-healpix + [--galactic] + [--outdir DIR] + [--parallelism FACTOR] [-j FACTOR] + {FITS-PATH} + {TOAST-DEPTH} + +The ``FITS-PATH`` argument gives the filename of the input image, which should +be a single FITS file storing data in the HEALPix_ format. + +The ``TOAST-DEPTH`` argument specifies the resolution level of the TOAST +pixelization that will be generated. A depth of 0 means that the input will be +sampled onto a single 256×256 tile, a depth of 1 means that the input will be +sampled onto four tiles for a total resolution of 512×512, and so on. When +converting from HEALPix, the appropriate TOAST depth can be derived from the +HEALPix ``N_side`` parameter as approximately:: + + depth = log2(N_side) - 6.2 + +One should typically round up to the next larger integer to err on the side of +higher resolution, although the best policy will depend on the sampling of the +input HEALPix map. + +The ``--outdir DIR`` option specifies where the output data should be written. +If unspecified, the data root will be the current directory. + +The ``--galactic`` option forces the tiling to assume that the input map is in +Galactic coordinates, even if the FITS headers do not specify that. This has +been observed in some data sets in the wild. + +The ``--parallelism FACTOR`` argument (or ``-j FACTOR``) specifies the level of +parallism to use. On operating systems that support parallel processing, the +default is to use all CPUs. To disable parallel processing, explicitly specify a +factor of 1. + +Notes +===== + +This command requires that the healpy_ Python module is installed. + +.. _healpy: https://healpy.readthedocs.io/ + +This command will create FITS files for the highest-resolution tile layer, +corresponding to the ``DEPTH`` argument, and emit an ``index_rel.wtml`` file +containing projection information and template metadata. + +The WWT rendering engine can also display datasets in the hierarchical HiPS_ +FITS format, which addresses about the same use case as TOAST. If a HiPS FITS +dataset already exists, it is probably a wiser choice to use it rather than +creating a new TOAST pyramid. However, TOAST provides noticeably better +performance for the initial data display and avoiding positional shifts as the +user zooms in. + +.. _HiPS: https://aladin.unistra.fr/hips/ + +Currently, parallel processing is only supported on the Linux operating system, +because ``fork()``-based multiprocessing is required. MacOS should support this, +but there is currently (as of Python 3.8) `a bug`_ preventing that. + +.. _a bug: https://bugs.python.org/issue33725 + + +See Also +======== + +- :ref:`cli-tile-allsky` \ No newline at end of file diff --git a/docs/toasting.rst b/docs/toasting.rst index b516997..30c0713 100644 --- a/docs/toasting.rst +++ b/docs/toasting.rst @@ -35,12 +35,12 @@ workflow in this case is as follows: which the curved surface of the sphere is mapped onto a 2D image. For full-sphere imagery, only a few choices are ever used: probably the main thing to check is whether your image is in equatorial (RA/Dec) or galactic - (l/b) coordinates. Consult the :ref:`cli-tile-allsky` documentation for the - supported choices. If your image uses an unsupported projection, please `file - a request`_ with the developers. + (l/b) coordinates. Consult the :ref:`cli-tile-allsky` documentation (or + :ref:`cli-tile-healpix`) for the supported choices. If your image uses an + unsupported projection, please `file a request`_ with the developers. -4. Determine the *tiling depth* appropriate for your use case. The depth is a number - that specifies the highest resolution that your final map will attain. +4. Determine the *tiling depth* appropriate for your use case. The depth is a + number that specifies the highest resolution that your final map will attain. Consult the :ref:`cli-tile-allsky` documentation for the quantitative definition. The best choice will depend on your individual circumstances. But as a general guideline, you should probably choose the depth that yields a @@ -124,7 +124,7 @@ workflow in this case is as follows: 11. Finally, upload the complete contents of your ``tiled`` subdirectory to your web server. In this case, the upload location should be such that the url - ``http://myserver.org/datasetname/index.wtml``_ will yield the + ``http://myserver.org/datasetname/index.wtml`` will yield the ``index.wtml`` file created in the previous step. .. _file a request: https://github.com/WorldWideTelescope/toasty/issues/ From 355ca4d70e5135736aef60ef5a0f7a12a9adad4b Mon Sep 17 00:00:00 2001 From: Peter Williams Date: Tue, 12 Jul 2022 09:10:57 -0400 Subject: [PATCH 08/11] docs/conf.py: reformat with "black" --- docs/conf.py | 89 ++++++++++++++++++++++------------------------------ 1 file changed, 37 insertions(+), 52 deletions(-) diff --git a/docs/conf.py b/docs/conf.py index 098a937..fa0f951 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -1,78 +1,63 @@ # -*- coding: utf-8 -*- -project = 'toasty' -author = 'Chris Beaumont and the AAS WorldWide Telescope Team' -copyright = '2014-2020, ' + author +project = "toasty" +author = "Chris Beaumont and the AAS WorldWide Telescope Team" +copyright = "2014-2020, " + author -release = '0.dev0' # cranko project-version -version = '.'.join(release.split('.')[:2]) +release = "0.dev0" # cranko project-version +version = ".".join(release.split(".")[:2]) extensions = [ - 'sphinx.ext.autodoc', - 'sphinx.ext.doctest', - 'sphinx.ext.intersphinx', - 'sphinx.ext.mathjax', - 'sphinx.ext.viewcode', - 'sphinx_automodapi.automodapi', - 'sphinx_automodapi.smart_resolver', - 'numpydoc', + "sphinx.ext.autodoc", + "sphinx.ext.doctest", + "sphinx.ext.intersphinx", + "sphinx.ext.mathjax", + "sphinx.ext.viewcode", + "sphinx_automodapi.automodapi", + "sphinx_automodapi.smart_resolver", + "numpydoc", ] -templates_path = ['_templates'] -source_suffix = '.rst' -master_doc = 'index' +templates_path = ["_templates"] +source_suffix = ".rst" +master_doc = "index" language = None -exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store'] -pygments_style = 'sphinx' +exclude_patterns = ["_build", "Thumbs.db", ".DS_Store"] +pygments_style = "sphinx" todo_include_todos = False -html_theme = 'bootstrap-astropy' +html_theme = "bootstrap-astropy" html_theme_options = { - 'logotext1': 'toasty', - 'logotext2': '', - 'logotext3': ':docs', - 'astropy_project_menubar': False, + "logotext1": "toasty", + "logotext2": "", + "logotext3": ":docs", + "astropy_project_menubar": False, } -html_static_path = ['_static'] -htmlhelp_basename = 'toastydoc' +html_static_path = ["_static"] +htmlhelp_basename = "toastydoc" intersphinx_mapping = { - 'python': ( - 'https://docs.python.org/3/', - (None, 'http://data.astropy.org/intersphinx/python3.inv') + "python": ( + "https://docs.python.org/3/", + (None, "http://data.astropy.org/intersphinx/python3.inv"), ), - - 'astropy': ( - 'https://docs.astropy.org/en/stable/', - None - ), - - 'numpy': ( - 'https://numpy.org/doc/stable/', - (None, 'http://data.astropy.org/intersphinx/numpy.inv') - ), - - 'PIL': ( - 'https://pillow.readthedocs.io/en/stable/', - None - ), - - 'wwt_data_formats': ( - 'https://wwt-data-formats.readthedocs.io/en/stable/', - None + "astropy": ("https://docs.astropy.org/en/stable/", None), + "numpy": ( + "https://numpy.org/doc/stable/", + (None, "http://data.astropy.org/intersphinx/numpy.inv"), ), + "PIL": ("https://pillow.readthedocs.io/en/stable/", None), + "wwt_data_formats": ("https://wwt-data-formats.readthedocs.io/en/stable/", None), } numpydoc_show_class_members = False nitpicky = True -nitpick_ignore = [ - ('py:class', 'ipywidgets.widgets.domwidget.DOMWidget') -] +nitpick_ignore = [("py:class", "ipywidgets.widgets.domwidget.DOMWidget")] -default_role = 'obj' +default_role = "obj" -html_logo = 'images/logo.png' +html_logo = "images/logo.png" linkcheck_retries = 5 linkcheck_timeout = 10 From 0f72bb87e94e5cc95df5529fd881fdaf53e62321 Mon Sep 17 00:00:00 2001 From: Peter Williams Date: Tue, 12 Jul 2022 09:11:12 -0400 Subject: [PATCH 09/11] docs/conf.py: set language --- docs/conf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/conf.py b/docs/conf.py index fa0f951..5749c25 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -21,7 +21,7 @@ templates_path = ["_templates"] source_suffix = ".rst" master_doc = "index" -language = None +language = "en" exclude_patterns = ["_build", "Thumbs.db", ".DS_Store"] pygments_style = "sphinx" todo_include_todos = False From 1d41e091ab7bff957fba3c9aec7b40027e3dc8b3 Mon Sep 17 00:00:00 2001 From: Peter Williams Date: Tue, 12 Jul 2022 09:11:44 -0400 Subject: [PATCH 10/11] Add --avm-from option to "toasty tile-study" --- docs/cli/tile-study.rst | 7 +++++++ toasty/cli.py | 22 +++++++++++++++++++--- toasty/tests/test_study.py | 31 ++++++++++++++++++++++++------- 3 files changed, 50 insertions(+), 10 deletions(-) diff --git a/docs/cli/tile-study.rst b/docs/cli/tile-study.rst index 986a0c4..1e56872 100644 --- a/docs/cli/tile-study.rst +++ b/docs/cli/tile-study.rst @@ -18,6 +18,7 @@ Usage [--outdir DIR] [--name NAME] [--avm] + [--avm-from PATH] [--fits-wcs PATH] IMAGE-PATH @@ -49,6 +50,12 @@ requires the `pyavm`_ package to be installed on your system. .. _pyavm: https://astrofrog.github.io/pyavm/ +The ``--avm-from PATH`` argument is like ``--avm``, but loads the AVM +information from a *different* image file. This can be useful if you have +multiple files with different versions of the same image, and the one that you +want to tile has broken or missing AVM. It is up to you to ensure that +information contained in the AVM file corresponds to the main input image. + If the ``--fits-wcs`` argument is given, Toasty will attempt to load world-coordinate information from the headers of the named `FITS`_ file. It is up to you to ensure that information contained in that file corresponds to the diff --git a/toasty/cli.py b/toasty/cli.py index cb8bdcd..92285a5 100644 --- a/toasty/cli.py +++ b/toasty/cli.py @@ -501,6 +501,11 @@ def tile_study_getparser(parser): action="store_true", help="Expect the input image to have AVM positioning tags", ) + parser.add_argument( + "--avm-from", + metavar="PATH", + help="Set positioning based on AVM tags in a different image", + ) parser.add_argument( "--fits-wcs", metavar="PATH", @@ -540,7 +545,7 @@ def tile_study_impl(settings): # Now deal with the WCS - if settings.avm: + if settings.avm or settings.avm_from: # We don't *always* check for AVM because pyavm prints out a lot of junk # and may not be installed. try: @@ -551,10 +556,21 @@ def tile_study_impl(settings): try: with warnings.catch_warnings(): warnings.simplefilter("ignore") - avm = AVM.from_image(settings.imgpath) + + if settings.avm_from: + if settings.avm: + print( + "warning: `--avm` option superseded by `--avm-from`", + file=sys.stderr, + ) + avm_input = settings.avm_from + else: + avm_input = settings.imgpath + + avm = AVM.from_image(avm_input) except Exception: print( - f"error: failed to read AVM tags of input `{settings.imgpath}`", + f"error: failed to read AVM tags of input `{avm_input}`", file=sys.stderr, ) raise diff --git a/toasty/tests/test_study.py b/toasty/tests/test_study.py index 1f37391..34c981b 100644 --- a/toasty/tests/test_study.py +++ b/toasty/tests/test_study.py @@ -170,8 +170,7 @@ def test_sample_cli(self): """ - @pytest.mark.skipif("not HAS_AVM") - def test_avm(self): + def basic_avm_test_helper(self, args): from xml.etree import ElementTree as etree # Hack for macOS: XML textualization is ever-so-slightly different. @@ -183,7 +182,16 @@ def test_avm(self): expected = etree.fromstring(wtml) - cli.entrypoint( + cli.entrypoint(args) + + with open(self.work_path("index_rel.wtml"), "rt", encoding="utf8") as f: + observed = etree.fromstring(f.read()) + + assert_xml_elements_equal(observed, expected) + + @pytest.mark.skipif("not HAS_AVM") + def test_avm(self): + self.basic_avm_test_helper( [ "tile-study", "--avm", @@ -193,10 +201,19 @@ def test_avm(self): ] ) - with open(self.work_path("index_rel.wtml"), "rt", encoding="utf8") as f: - observed = etree.fromstring(f.read()) - - assert_xml_elements_equal(observed, expected) + @pytest.mark.skipif("not HAS_AVM") + def test_avm_from(self): + # Dumb smoke test here + self.basic_avm_test_helper( + [ + "tile-study", + "--avm-from", + test_path("geminiann11015a.jpg"), + "--outdir", + self.work_path(), + test_path("geminiann11015a.jpg"), + ] + ) # Blah: this is the same as AVM_WTML, minus the metadata fields that aren't # preserved in the FITS export. From 93396521eb1a31235572de7e4d15c764e1525730 Mon Sep 17 00:00:00 2001 From: Peter Williams Date: Tue, 12 Jul 2022 09:30:46 -0400 Subject: [PATCH 11/11] Fix up tests for latest JPEG2000/Numpy --- toasty/image.py | 35 +++++++++++++++++++++++++---------- toasty/tests/test_toast.py | 11 ++++++----- 2 files changed, 31 insertions(+), 15 deletions(-) diff --git a/toasty/image.py b/toasty/image.py index fb21258..8f8201a 100644 --- a/toasty/image.py +++ b/toasty/image.py @@ -1022,6 +1022,29 @@ def ensure_negative_parity(self): self.flip_parity() return self + def _as_writeable_array(self): + """ + Helper that does what it says. Should potentially become a public API at + some point? + """ + + arr = self.asarray() + + if not arr.flags.writeable: + try: + arr.setflags(write=True) + except ValueError: + # Source array is immutable, perhaps because it does not own its + # data (OWNDATA flag is false). In that case we have to copy: + arr = arr.copy() + self._array = arr + + # Ensure that we don't try to use the PIL representation anymore, since + # it will be out-of-date. + self._pil = None + + return arr + def fill_into_maskable_buffer(self, buffer, iy_idx, ix_idx, by_idx, bx_idx): """ Fill a maskable buffer with a rectangle of data from this image. @@ -1049,11 +1072,7 @@ def fill_into_maskable_buffer(self, buffer, iy_idx, ix_idx, by_idx, bx_idx): """ i = self.asarray() - b = buffer.asarray() - - # Ensure that we don't try to use the PIL representation of the buffer anymore, - # since it will be out-of-date. - buffer._pil = None + b = buffer._as_writeable_array() if self.mode == ImageMode.RGB: b.fill(0) @@ -1093,11 +1112,7 @@ def update_into_maskable_buffer(self, buffer, iy_idx, ix_idx, by_idx, bx_idx): """ i = self.asarray() - b = buffer.asarray() - - # Ensure that we don't try to use the PIL representation of the buffer anymore, - # since it will be out-of-date. - buffer._pil = None + b = buffer._as_writeable_array() sub_b = b[by_idx, bx_idx] sub_i = i[iy_idx, ix_idx] diff --git a/toasty/tests/test_toast.py b/toasty/tests/test_toast.py index e2664d2..b4074c0 100644 --- a/toasty/tests/test_toast.py +++ b/toasty/tests/test_toast.py @@ -305,12 +305,12 @@ def test_planet_zeroleft(self): test_path("Equirectangular_projection_SW-tweaked.jpg") ) hw = img.width // 2 - buf = img.asarray() + + buf = img._as_writeable_array() temp = buf[:, hw:].copy() buf[:, hw:] = buf[:, :hw] buf[:, :hw] = temp - img._array = buf # hack! - img._pil = None + p = self.work_path("zeroleft.jpg") img.save(p, format="jpg") self.inner_test( @@ -327,12 +327,13 @@ def test_planet_zeroright(self): test_path("Equirectangular_projection_SW-tweaked.jpg") ) hw = img.width // 2 - buf = img.asarray() + + buf = img._as_writeable_array() temp = buf[:, hw:].copy() buf[:, hw:] = buf[:, :hw] buf[:, :hw] = temp img._array = buf[:, ::-1] # hack! - img._pil = None + p = self.work_path("zeroright.jpg") img.save(p, format="jpg") self.inner_test(