Merge pull request #61 from pkgw/fits-study
Fix up tile-study with FITS files
This commit is contained in:
Peter Williams
GitHub
Коммит
3d23dd7ec7
|
|
@ -104,6 +104,8 @@ jobs:
|
|||
pyavm \
|
||||
pytest-cov \
|
||||
pyyaml \
|
||||
reproject \
|
||||
shapely \
|
||||
tqdm
|
||||
pip install wwt_data_formats
|
||||
pip install -e .
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -71,7 +71,74 @@ class Builder(object):
|
|||
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()
|
||||
|
|
|
|||
|
|
@ -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()
|
||||
|
||||
|
|
|
|||
|
|
@ -168,8 +168,8 @@ class MultiTanProcessor(object):
|
|||
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()
|
||||
|
|
|
|||
|
|
@ -320,6 +320,12 @@ class StudyTiling(object):
|
|||
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 @@ class StudyTiling(object):
|
|||
|
||||
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)
|
||||
|
|
|
|||
|
|
@ -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 collection
|
|||
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 = """
|
||||
<Folder Group="Explorer" Name="TestName">
|
||||
<Place Angle="0" DataSetType="Sky" Dec="0.74380165289257" Name="TestName"
|
||||
Opacity="100" RA="14.419753086419734" Rotation="0.0"
|
||||
ZoomLevel="0.1433599999999144">
|
||||
<Folder Browseable="True" Group="Explorer" Name="Toasty" MSRCommunityId="0"
|
||||
MSRComponentId="0" Permission="0" Searchable="True" Type="Sky">
|
||||
<Place Angle="0.0" AngularSize="0.0" DataSetType="Sky" Dec="0.7438249862258411"
|
||||
Distance="0.0" DomeAlt="0.0" DomeAz="0.0" Lat="0.0" Lng="0.0"
|
||||
Magnitude="0.0" MSRCommunityId="0" MSRComponentId="0" Name="Toasty"
|
||||
Opacity="100.0" Permission="0" RA="14.41975153073335" Rotation="0.0"
|
||||
Thumbnail="thumb.jpg" ZoomLevel="0.2437119999998555" >
|
||||
<ForegroundImageSet>
|
||||
<ImageSet BandPass="Gamma" BaseDegreesPerTile="0.023893333333319066"
|
||||
BaseTileLevel="0" BottomsUp="False" CenterX="216.296296296296"
|
||||
CenterY="0.74380165289257" DataSetType="Sky" FileType=".png"
|
||||
Name="TestName" OffsetX="4.66666666666388e-05"
|
||||
OffsetY="4.66666666666388e-05" Projection="Tan" Rotation="0.0"
|
||||
Sparse="True" TileLevels="1" Url="UP{1}/{3}/{3}_{2}.png"
|
||||
WidthFactor="2">
|
||||
<Credits>CT</Credits>
|
||||
<CreditsUrl>CU</CreditsUrl>
|
||||
<ThumbnailUrl>TU</ThumbnailUrl>
|
||||
<Description>DT</Description>
|
||||
<ImageSet BandPass="Visible" BaseDegreesPerTile="0.023893333333319167"
|
||||
BaseTileLevel="0" BottomsUp="False" CenterX="216.2962962963"
|
||||
CenterY="0.74380165289257" DataSetType="Sky" ElevationModel="False"
|
||||
FileType=".fits" Generic="False" MeanRadius="0.0" MSRCommunityId="0"
|
||||
MSRComponentId="0" Name="Toasty" OffsetX="2.33333333333195e-05"
|
||||
OffsetY="2.33333333333195e-05" Permission="0" Projection="Tan"
|
||||
Rotation="-0.0" Sparse="True" StockSet="False" TileLevels="1"
|
||||
Url="{1}/{3}/{3}_{2}.fits" WidthFactor="2">
|
||||
<ThumbnailUrl>thumb.jpg</ThumbnailUrl>
|
||||
</ImageSet>
|
||||
</ForegroundImageSet>
|
||||
</Place>
|
||||
</Folder>"""
|
||||
|
||||
# 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.
|
||||
# Gross workaround for platform differences in the XML output.
|
||||
|
||||
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"')
|
||||
)
|
||||
if sys.platform == 'darwin':
|
||||
WTML = WTML.replace('Dec="0.7438249862258411"', 'Dec="0.743824986225841"')
|
||||
|
||||
# Back to the non-gross stuff.
|
||||
|
||||
|
|
@ -67,6 +74,7 @@ class TestMultiTan(object):
|
|||
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 @@ class TestMultiTan(object):
|
|||
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('<ThumbnailUrl>thumb.jpg</ThumbnailUrl>', '')
|
||||
)
|
||||
|
||||
args = [
|
||||
'tile-multi-tan',
|
||||
|
|
@ -94,3 +146,77 @@ class TestMultiTan(object):
|
|||
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)
|
||||
|
|
|
|||
Загрузка…
Ссылка в новой задаче