Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Handle aperture photometry when images linked by WCS for advanced cases #2154

Merged
merged 8 commits into from
Aug 14, 2023

Conversation

pllim
Copy link
Contributor

@pllim pllim commented Apr 14, 2023

Description

This pull request is to investigate and hopefully fix the bug of aperture photometry not able to provide correct result in certain cases, e.g.:

  • When images have different pixel scales.
  • When Subset no longer retains original shape after being projected onto a non-reference image.

Blocked by

TODO

Change log entry

  • Is a change log needed? If yes, is it added to CHANGES.rst? If you want to avoid merge conflicts,
    list the proposed change log here for review and add to CHANGES.rst before merge. If no, maintainer
    should add a no-changelog-entry-needed label.

Checklist for package maintainer(s)

This checklist is meant to remind the package maintainer(s) who will review this pull request of some common things to look for. This list is not exhaustive.

  • Are two approvals required? Branch protection rule does not check for the second approval. If a second approval is not necessary, please apply the trivial label.
  • Do the proposed changes actually accomplish desired goals? Also manually run the affected example notebooks, if necessary.
  • Do the proposed changes follow the STScI Style Guides?
  • Are tests added/updated as required? If so, do they follow the STScI Style Guides?
  • Are docs added/updated as required? If so, do they follow the STScI Style Guides?
  • Did the CI pass? If not, are the failures related?
  • Is a milestone set? Set this to bugfix milestone if this is a bug fix and needs to be released ASAP; otherwise, set this to the next major release milestone.
  • After merge, any internal documentations need updating (e.g., JIRA, Innerspace)? 🐱 🐱 🐱

@pllim pllim added this to the 3.5 milestone Apr 14, 2023
@github-actions github-actions bot added the documentation Explanation of code and concepts label Apr 14, 2023
"axs[1].set_ylim(70, 80)\n",
"axs[2].imshow(rect_grp.subsets[2].to_mask(), origin='lower')\n",
"axs[2].set_xlim(105, 125)\n",
"axs[2].set_ylim(70, 80);"
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This would give:

Screenshot 2023-04-14 173003

This corresponds to the following images loaded into Imviz and linked by WCS:

Screenshot 2023-04-14 173442

This comment was marked as outdated.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@larrybradley , to answer your question, if I draw a rectangle on the third image (down-sampled and rotated), the rectangle comes back warped as if I drew it on the reference image, which is what I suspected would happen.

Screenshot 2023-04-17 172540

@pllim
Copy link
Contributor Author

pllim commented Apr 17, 2023

@larrybradley and I discussed several options. The most scientifically correct way is still using sky regions (as suggested by @bmorris3). However, this does disable the plugin for any data without WCS.

Hence, Larry also requested that we put in extra logic to still fall back to the old way of doing things (pixel regions) if image does not have WCS and is the reference image for the Subset. But is that the same as a viewer reference image (have to check what becomes the viewer reference image when you only display a non-reference image in the viewer; answer is still the original reference image)?

While we're at it, I guess technically we can always falls back to using pixel regions if we detect the link is in pixel space whether the data of interest has WCS or not. Keep in mind that even when Imviz is asked to link by WCS and the reference data has WCS, any image without WCS still links to the reference data by pixel.

And since glue-astronomy does not handle sky regions (see glue-viz/glue-astronomy#89), we would have to rely on our own region_translators.py.

@pllim
Copy link
Contributor Author

pllim commented Apr 18, 2023

For completeness, Larry said photutils would ignore distortion anyway (that is, if you pass in a circular aperture and WCS with distortion, it would still be circular when extracting the data array, not a distorted circle shape) so we should not worry about distortion here.

@pllim
Copy link
Contributor Author

pllim commented Apr 21, 2023

For future me: With this, I don't need GWCS.

And since I will add all these as static test data files, hopefully we never have to rerun these again.

To downsample by 2

w2_hdr['CRPIX1'] = 250 / 2
w2_hdr['CRPIX2'] = 150 / 2
w2_hdr['NAXIS1'] = 250
w2_hdr['NAXIS2'] = 150
w2_hdr['CDELT1'] = 2
w2_hdr['CDELT2'] = 2

Goes with:

data_sm = block_reduce(data, 2)

Where:

data = make_100gaussians_image()
w = make_wcs(data.shape)

To rotate

https://github.com/musevlt/mpdaf/blob/321cec4e95dc70c3e7eaf7402a2e3a24ebfb5b8e/lib/mpdaf/obj/coords.py#L1550

def rotate(wcs_in, theta):
    wcs = deepcopy(wcs_in)
    theta = np.deg2rad(theta)
    sinq = np.sin(theta)
    cosq = np.cos(theta)
    mrot = np.array([[cosq, -sinq],
                     [sinq, cosq]])

    if wcs.wcs.has_cd():    # CD matrix
        newcd = np.dot(mrot, wcs.wcs.cd)
        wcs.wcs.cd = newcd
        wcs.wcs.set()
    elif wcs.wcs.has_pc():      # PC matrix + CDELT
        newpc = np.dot(mrot, wcs.wcs.get_pc())
        wcs.wcs.pc = newpc
        wcs.wcs.set()
    else:
        raise TypeError("Unsupported wcs type (need CD or PC matrix)")
    
    return wcs

Goes with:

data_rot, _ = reproject_interp(image_2, w3, shape_out=data_sm.shape)

@pllim pllim added bug Something isn't working imviz labels Apr 21, 2023
@pllim

This comment was marked as resolved.

@haticekaratay haticekaratay modified the milestones: 3.5, 3.6 May 25, 2023
@javerbukh javerbukh modified the milestones: 3.6, 3.7 Jul 28, 2023
@pllim pllim force-pushed the aper-subset-not-ref branch 3 times, most recently from f6a38a9 to cfb3ff1 Compare August 8, 2023 22:04
@codecov
Copy link

codecov bot commented Aug 10, 2023

Codecov Report

Patch coverage: 100.00% and project coverage change: +0.02% 🎉

Comparison is base (0c97b76) 90.70% compared to head (3ab15d3) 90.72%.

Additional details and impacted files
@@            Coverage Diff             @@
##             main    #2154      +/-   ##
==========================================
+ Coverage   90.70%   90.72%   +0.02%     
==========================================
  Files         157      157              
  Lines       17969    18015      +46     
==========================================
+ Hits        16299    16345      +46     
  Misses       1670     1670              
Files Changed Coverage Δ
...igs/default/plugins/subset_plugin/subset_plugin.py 93.64% <100.00%> (+0.10%) ⬆️
jdaviz/configs/imviz/helper.py 97.35% <100.00%> (+0.21%) ⬆️
...imviz/plugins/aper_phot_simple/aper_phot_simple.py 92.87% <100.00%> (+0.21%) ⬆️
jdaviz/configs/imviz/plugins/viewers.py 89.65% <100.00%> (-0.86%) ⬇️
...daviz/configs/imviz/tests/test_simple_aper_phot.py 100.00% <100.00%> (ø)
jdaviz/core/region_translators.py 98.92% <100.00%> (-0.09%) ⬇️

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

@pllim pllim marked this pull request as ready for review August 10, 2023 03:54
Copy link
Contributor

@bmorris3 bmorris3 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Neat!

jdaviz/configs/imviz/helper.py Outdated Show resolved Hide resolved
Comment on lines +282 to 289
# TODO: Brett Morris might want to look at this for
# https://github.com/spacetelescope/jdaviz/pull/2179
# ref_label = self.state.reference_data ???
#
# The original links were created against data_collection[0], not necessarily
# against the current viewer reference_data
ref_label = self.session.application.data_collection[0].label
Copy link
Contributor

@bmorris3 bmorris3 Aug 10, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If I understand your comment correctly, you're thinking about (near-)future-proofing this by looking up the current reference data's label? If so, I'll suggest that change here.

jdaviz/configs/imviz/plugins/viewers.py Show resolved Hide resolved
jdaviz/core/region_translators.py Outdated Show resolved Hide resolved
jdaviz/core/region_translators.py Outdated Show resolved Hide resolved
Copy link
Contributor

@bmorris3 bmorris3 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for the updates 👌🏻

Comment on lines +81 to +87
<v-text-field v-if="item.name === 'Parent'"
:label="item.name"
:value="item.value"
style="padding-top: 0px; margin-top: 0px"
:readonly="true"
hint="Subset was defined with respect to this reference data (read-only)"
></v-text-field>
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we'll eventually want to move this out of the subregions once reparenting is implemented (or perhaps remove entirely), but since its read-only, it doesn't hurt to have for now and could be useful for debugging.

Comment on lines +219 to +222
"""Find the type of ``glue`` linking between the given
data labels. A link is bi-directional. If there are
more than 2 data in the collection, one of the given
labels should be the reference data or look-up will fail.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@bmorris3 - does the rotation work remove some of this flexibility (by forbidding mixed linking)? If we might need to change or remove this, should we make it private in the meantime?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this method is useful with or without the rotation work. It's likely the case that it will be redundant with the link type attribute after image rotation is merged and WCS-linked in Imviz, but I think it's handy for debugging.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you think it should be public or private in the meantime?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't see why it should be private. The information is so useful in the meantime that I think it should be public. If it comes a day we don't need it anymore, then we can deprecate it.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree that public is fine.

@pllim
Copy link
Contributor Author

pllim commented Aug 14, 2023

There is a conflict now, I will have to rebase and let the CI run again.

pllim and others added 8 commits August 14, 2023 12:05
advanced aperture photometry.

Add tests that will fail without fix

WIP: How to actually fix this mess

Remove warning from doc

[ci skip] [rtd skip]
BUG: Make sure recentering works in Subset Tools

WIP: What about the actual bug

[ci skip] [rtd skip]
to images linked by WCS if sky projection is different
and remove outdated comment
and remove debug comment
Co-authored-by: Brett M. Morris <[email protected]>
@pllim
Copy link
Contributor Author

pllim commented Aug 14, 2023

Looks like a new deprecation warning from scipy over the weekend. It is unrelated.

`numpy.distutils` is deprecated since NumPy 1.23.0, as a result

@pllim pllim merged commit 7dcc714 into spacetelescope:main Aug 14, 2023
13 of 14 checks passed
@pllim pllim deleted the aper-subset-not-ref branch August 14, 2023 17:04
@pllim
Copy link
Contributor Author

pllim commented Aug 14, 2023

Thanks, all!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working documentation Explanation of code and concepts imviz Ready for final review
Projects
None yet
5 participants