diff --git a/CHANGELOG.md b/CHANGELOG.md index 6f9beb8..92208a9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,16 @@ +# toasty 0.17.0 (2022-07-12) + +- Add an `--avm-from` option to `toasty tile-study`, which would have been + useful with the JWST imagery released yesterday (#82, @pkgw). +- Support tiled FITS in `toasty tile-healpix` (#79, #80, @pkgw), and add a + `--force-galactic` option. +- Avoid HEALPix HDUs without data (#79, @pkgw) +- Add some diagnostics for spatial AVM information in `toasty check-avm` (#78, + @pkgw) +- Update tests for new `wwt_data_formats` constellation support and other + changing dependencies. + + # toasty 0.16.1 (2022-01-27) - Toasty is now more forgiving with FITS and WCS shapes when tiling FITS data 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/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/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/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/docs/conf.py b/docs/conf.py index 098a937..5749c25 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' -language = None -exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store'] -pygments_style = 'sphinx' +templates_path = ["_templates"] +source_suffix = ".rst" +master_doc = "index" +language = "en" +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 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/ diff --git a/setup.py b/setup.py index 31a72a6..667a682 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 @@ -39,7 +39,7 @@ def get_long_desc(): setup_args = dict( name="toasty", # cranko project-name - version="0.16.1", # cranko project-version + version="0.17.0", # cranko project-version description="Generate TOAST image tile pyramids from existing image data", long_description=get_long_desc(), long_description_content_type="text/markdown", @@ -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/cli.py b/toasty/cli.py index d78d857..92285a5 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 @@ -318,12 +366,24 @@ 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", 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", @@ -342,10 +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") - sampler = healpix_fits_file_sampler(settings.fitspath) + 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( @@ -434,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", @@ -473,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: @@ -484,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/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/samplers.py b/toasty/samplers.py index acb1e50..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 @@ -98,13 +98,20 @@ 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()) -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 @@ -118,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 ------- @@ -143,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) 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" 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): """ - @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. @@ -182,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", @@ -192,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. @@ -204,6 +222,7 @@ def test_avm(self): 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) 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) @@ -647,15 +680,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) @@ -686,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) @@ -696,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) @@ -738,25 +782,45 @@ 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 + ) + + # 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): + 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(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 +831,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 +855,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) @@ -816,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) @@ -826,7 +901,15 @@ 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(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) + )