radio-astro-tools / spectral-cube

Library for reading and analyzing astrophysical spectral data cubes
http://spectral-cube.rtfd.org
BSD 3-Clause "New" or "Revised" License
98 stars 65 forks source link

Add cube mosaicking #829

Closed SpacialTree closed 2 years ago

SpacialTree commented 2 years ago

Created functions to mosaic cubes together.

codecov-commenter commented 2 years ago

Codecov Report

Base: 77.88% // Head: 77.36% // Decreases project coverage by -0.51% :warning:

Coverage data is based on head (d50faf2) compared to base (a37fbaa). Patch coverage: 6.81% of modified lines in pull request are covered.

Additional details and impacted files ```diff @@ Coverage Diff @@ ## master #829 +/- ## ========================================== - Coverage 77.88% 77.36% -0.52% ========================================== Files 24 24 Lines 5982 6026 +44 ========================================== + Hits 4659 4662 +3 - Misses 1323 1364 +41 ``` | [Impacted Files](https://codecov.io/gh/radio-astro-tools/spectral-cube/pull/829?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=radio-astro-tools) | Coverage Δ | | |---|---|---| | [spectral\_cube/spectral\_cube.py](https://codecov.io/gh/radio-astro-tools/spectral-cube/pull/829/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=radio-astro-tools#diff-c3BlY3RyYWxfY3ViZS9zcGVjdHJhbF9jdWJlLnB5) | `76.49% <0.00%> (-0.11%)` | :arrow_down: | | [spectral\_cube/cube\_utils.py](https://codecov.io/gh/radio-astro-tools/spectral-cube/pull/829/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=radio-astro-tools#diff-c3BlY3RyYWxfY3ViZS9jdWJlX3V0aWxzLnB5) | `74.38% <7.14%> (-8.69%)` | :arrow_down: | Help us with your feedback. Take ten seconds to tell us [how you rate us](https://about.codecov.io/nps?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=radio-astro-tools). Have a feature suggestion? [Share it here.](https://app.codecov.io/gh/feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=radio-astro-tools)

:umbrella: View full report at Codecov.
:loudspeaker: Do you have feedback about the report comment? Let us know in this issue.

SpacialTree commented 2 years ago

@keflavich Ready for further review!

keflavich commented 2 years ago

It's time to start looking at the error messages. For example, this is a syntax error: https://github.com/radio-astro-tools/spectral-cube/runs/7883143970?check_suite_focus=true#step:5:162

keflavich commented 2 years ago

We need to have the test skipped if reproject is not installed

keflavich commented 2 years ago

OK we're really close now, but the WCSes don't match up https://github.com/radio-astro-tools/spectral-cube/runs/7996917393?check_suite_focus=true#step:5:596

SpacialTree commented 2 years ago

There is an issue where the edge pixels of the image are turned into NaNs when reprojected. For very small cubes, such as those used in tests, this may result in the entire cube being turned into NaNs.

edges

Code to reproduce:

from spectral_cube import SpectralCube
import numpy as np
import matplotlib.pyplot as plt

# Read in a spectral cube
cubename = '/orange/adamginsburg/cmz/g5/G5/sum/newcombination/g5.H2CO.spw23_test.fits' # Example data
cube = SpectralCube.read(cubename)

# Make 2D example so it will run faster, same result of edges being lost when using a cube
mom0 = cube.moment0() # Moment 0 of original cube
repro_mom0 = mom0.reproject(mom0.header) # Reproject Moment 0 to same header

# Get the masks of the Moment 0's (there is no get_mask_array for Projection objects)
mom0_mask = np.isfinite(mom0)
repro_mom0_mask = np.isfinite(repro_mom0)

# Find where the masks are different with xor
edges = np.logical_xor(mom0_mask, repro_mom0_mask)

# Plot an image of the edges
plt.imshow(edges)
keflavich commented 2 years ago

Good catch. This isn't really a bug as much as an annoying feature. I think we can work around it by specifying order='nearest-neighbor' in the call to reproject.

SpacialTree commented 2 years ago

Looks like that fixed the issue for moment 0's, but not for cubes.

Try this for a cube:

from spectral_cube import SpectralCube
import numpy as np
import matplotlib.pyplot as plt

# Read in a spectral cube
cubename = '/orange/adamginsburg/cmz/g5/G5/sum/newcombination/g5.H2CO.spw23_test.fits' # Example data
cube = SpectralCube.read(cubename)

cutcube = cube[:1]
repro_cutcube = cutcube.reproject(cutcube.header)

cutcube_mask = np.isfinite(cutcube)
repro_cutcube_mask = np.isfinite(repro_cutcube)

edges = np.logical_xor(cutcube_mask, repro_cutcube_mask)

plt.imshow(edges[0])
keflavich commented 2 years ago

Still looks ok if you pass reproject(..., order='nearest-neighbor'), no?

SpacialTree commented 2 years ago

Looks like I missed the second reproject in the code, should be fixed now!

keflavich commented 2 years ago

It looks like tests are failing only on windows, so this is now good to go. However, I want to see if I can fix the windows tests first (I think main is also failing on windows now)

keflavich commented 2 years ago

I manage to fix the tests on windows, but now something has gone very weird with this PR - there is some documentation stuff showing up that has no right to be.

keflavich commented 2 years ago

We're seeing the 'all-nans' errors on Windows only, I think, but it is coming from this PR.

keflavich commented 2 years ago

OK, I've pinned this down to a really obscene error:

(Pdb) wcs_out.spectral.pixel_to_world(np.arange(6))
<SpectralCoord
   (doppler_rest=0.0 m
    doppler_convention=optical)
  [-321214.698632  , -319926.48366321, -318638.26869442, -317350.05372563,
   -316061.83875684, -314773.62378805] m / s>
(Pdb) wcs_in.spectral.world_to_pixel(wcs_out.spectral.pixel_to_world(np.arange(6)))
array([nan, nan, nan, nan, nan, nan])

which arises because of a deep, stupid, painful bug in WCS:

(Pdb) wcs_in.wcs.restfrq
1420405718.41
(Pdb) wcs_in.wcs.restwav
0.0

the solution is apparently to hack in a restwav...

keflavich commented 2 years ago

https://github.com/astropy/astropy/issues/13774 is the issue that caused our nan problems.

keflavich commented 2 years ago

Well dangit, that's the same error.

keflavich commented 2 years ago

@astrofrog @e-koch we're getting "file not closed" errors in addition to other unexpected errors. Any idea why? It's in test_reproject. https://github.com/radio-astro-tools/spectral-cube/actions/runs/3176766813/jobs/5176439021#step:5:235

I can't reproduce the other errors, which is disconcerting. Looks like dev versions are the problem, so I'll try again with newer astropy etc.

keflavich commented 2 years ago

Nevermind the previous message - the "Errors" were caused by the legitimate failures. All green now!

@e-koch want a final look? I'm good to merge now