Enhanced results feature in ALMA #2855
Conversation
andamian
changed the title
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #2855 +/- ##
==========================================
- Coverage 67.35% 67.24% -0.11%
==========================================
Files 236 231 -5
Lines 18320 18255 -65
==========================================
- Hits 12339 12276 -63
+ Misses 5981 5979 -2 ☔ View full report in Codecov by Sentry. |
|
Yes, that looks like a good solution. |
|
@keflavich besides |
|
Yes, just stitching them together in a list is the way to go. Regions do have built-in exclude/include logic: the |
|
You can also consider making |
|
Thanks. I've somehow missed the "Combining Regions" section. I'll use |
|
Side note: what about supporting MOC? |
|
@keflavich - I'm wondering if those "Not Polygon" are holes in the bigger Polygon. If they are fully "contained" inside the big Polygon I can apply XOR to get the Compound sky region. But if they are not, I don't know how to get A-B when regions seem to offer only |
|
Any logic can be constructed from those operators, right? So you could do |
|
Indeed. Haven’t had my morning coffee yet. Thanks
…On Tue, Oct 24, 2023 at 9:15 AM Adam Ginsburg ***@***.***> wrote:
Any logic can be constructed from those operators, right? So you could do (a
or b) xor b if you want a and not b.
—
Reply to this email directly, view it on GitHub
<#2855 (comment)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/ABOTQLWDBHNYO6YMRCWHKKLYA7SQVAVCNFSM6AAAAAA6GJ3ALSVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTONZXGU3TKOBZGU>
.
You are receiving this because you were mentioned.Message ID:
***@***.***>
|
|
@keflavich - I've added support for all the shapes that are currently in the ALMA archive. I know this because I've run the parser over each one of them. Do you have any suggestions on how to test the more complex shapes ( |
|
Make a fake image that overlaps with the region.... e.g., using a lot of hasty off-the-top-of-my-head pseudocode: fake = np.ones([500,500])
ww = WCS()
ww.wcs.crpix=[250,250]
ww.wcs.crval=[center_ra, center_dec]
ww.wcs.cdelt=[1/3600,1/3600.] # or something appropriate for the footprint
ww.wcs.ctype=['RA---TAN', 'DEC--TAN']
preg = reg.to_pixel(ww)
masked = preg.to_mask().multiply(fake)
pl.imshow(masked) |
astroquery/alma/core.py
Outdated
| @@ -20,6 +20,7 @@ | |||
| from astropy import units as u | |||
| from astropy.time import Time | |||
| from astropy.coordinates import SkyCoord | |||
| import regions | |||
There was a problem hiding this comment.
regions is optional dependency, so do the same try/except import the other modules using it do (mocserver and xmatch, or even here in alma)
Where can I find the appropriate values for the wcs? Are they anywhere in the metadata or just in the FITS files? And if the latter, which of the FITS files? There are a lot of files corresponding to each of those observations. Thanks @keflavich |
|
All of the FITS files should have the same celestial footprint, so: from astropy.wcs import WCS
from astropy.io import fits
ww = WCS(fits.getheader('any_alma_fits_file.fits')).celestialthat WCS will work |
|
@keflavich - making little progress with the |
|
I agree that plotting is an important feature - I think we can push to get astropy/regions#470 completed. But I don't understand "the issue is that other shape operations such as UNION are still not supported": the PR provides the workaround for union'd shapes that is to show all of them individually. I think the general case of union of regions is not easily solvable. Almost all use cases should involve building masks - even for display, making a mask and showing contours may be a better approach. We could include examples of how to do that in documentation. How many observations in the archive are these corner-case multi-polygons? |
|
This case for example (https://almascience.eso.org/aq/?observationsProjectCode=2021.1.00547.S) with this shape: This is the Your comment about mask makes me think of MOCs. This could be the solution going forward but requires changes in the spec. It would be easier for push for such a change if we show that it works with contours. In ALMA there are other examples of multi polygon shapes but we don't know if there are valid or not. More recent example: https://almascience.eso.org/aq/?observationsProjectCode=2022.1.01742.S. We need to deal with them somehow. |
|
For my sake, I'm including graphics here: 
is the second (https://almascience.eso.org/aq/?observationsProjectCode=2022.1.01742.S). I think this one doesn't make sense - it's not a polygon, it's 3 individual single pointings, so it's weird to parametrize this way. It's because the target is a moving object wrt to the fixed sky frame (Titan).
Is there an Exception here? I think I missed this above; if there's a bug in regions, we should just fix that, but I'm not sure if that's where we are yet. |
|
The problem is here: https://github.com/astropy/regions/blob/1c88fa31e3685068b2ee42822e4efa22f49649bb/regions/core/compound.py#L147 Union is not expected in a compound region. |
|
Right; I think what we should be doing here is either: |
|
So I expect the region to be returned will be equivalent to this: import regions
from astropy.coordinates import SkyCoord
from astropy import units as u
from astropy.wcs import WCS
x='''53.095155 -27.862094 53.079461 -27.868350 53.070554 -27.864641 53.071634 -27.854716 53.063795 -27.848756 53.065662 -27.840191 53.057830 -27.834233 53.059697 -27.825669 53.051864 -27.819711 53.053728 -27.811149 53.045898 -27.805188 53.047766 -27.796624 53.040165 -27.791131 53.040727 -27.783321 53.168803 -27.740304 53.176556 -27.741168 53.180746 -27.747006 53.178683 -27.754479 53.186636 -27.761019 53.183259 -27.768949 53.192133 -27.774067 53.189238 -27.783465 53.198118 -27.788584 53.196264 -27.797148 53.204102 -27.803100 53.202244 -27.811668 53.209857 -27.817148 53.209307 -27.824959'''
ra, dec = list(map(float, x.split()[::2])), list(map(float, x.split()[1::2]))
coords = SkyCoord(ra, dec, frame='icrs', unit=(u.deg, u.deg))
poly1 = regions.PolygonSkyRegion(coords)
y = '53.160470 -27.801462 53.170993 -27.798574 53.173032 -27.789459 53.183944 -27.785173 53.175851 -27.781416 53.177962 -27.770657 53.168833 -27.764876 53.148836 -27.770820 53.146803 -27.779935 53.136929 -27.783386 53.144761 -27.789340 53.143267 -27.798785 53.153398 -27.795816'
ra2, dec2 = list(map(float, y.split()[::2])), list(map(float, y.split()[1::2]))
coords2 = SkyCoord(ra2, dec2, frame='icrs', unit=(u.deg, u.deg))
poly2 = regions.PolygonSkyRegion(coords2)
comp = (poly1 | poly2) ^ poly2One can then use this with, e.g.: ww = WCS(naxis=2) # it would be best to grab this from an ALMA fits file
ww.wcs.crval=53.1, -27.8
ww.wcs.cdelt[:] = 1/3600
ww.wcs.crpix = [512,512]
ww.wcs.ctype=['RA---TAN', 'DEC--TAN']
ww._naxis=[1024,1024]
pcomp = comp.to_pixel(ww)
mask = pcomp.to_mask()
ax = pl.subplot(projection=ww)
ax.imshow(mask)to get I think we'd ideally like to see something like this supported by regions, though: ax = pl.subplot(projection=ww)
ppoly1 = poly1.to_pixel(ww)
ppoly1.plot(ax=ax)
ppoly2 = poly2.to_pixel(ww)
ppoly2.plot(ax=ax, facecolor=(1,0,0,0.1), hatch='x', edgecolor='r')
ax.axis([0, 1024, 0, 1024]);
|
|
I think this should do it but please let me know if you think it's not the case. We need to wait for a |
docs/alma/alma.rst
Outdated
| The above footprint could be transformed into a pixel region and have the mask | ||
| extracted or plotted. Refer to the Astropy affiliated ''regions'' package | ||
| for more details. |
There was a problem hiding this comment.
We should try to point to a more specific example here; maybe we can build an astropy tutorial showing how to do this.
There was a problem hiding this comment.
I was thinking of doing that except that Sphinx always intimidates me :-). It will try to execute the code in CI so we have to make sure it can do this. I have this that I used for testing (and it's commented out in test_alma.py:
import matplotlib.pyplot as plt
artist = pix_region.as_artist()
axes = plt.subplot(projection=ww)
axes.set_aspect('equal')
axes.add_artist(artist)
axes.set_xlim([x_min, x_max])
axes.set_ylim([y_min, y_max])
pix_region.plot()
plt.show()
Should I give it a try? I agree that a graph will grab attention in a help document.
Alternatively, I can ask ALMA to do a notebook on this which we can then link. (see #2734)
There was a problem hiding this comment.
Yes, give it a try - we can help debug.
It would be great if ALMA could also point people in this direction with a notebook, though.
There was a problem hiding this comment.
OK with a simply Polygon or Circle? I cannot display more interesting Compound regions like the one I have in the documentation because of this condition:
https://github.com/astropy/regions/blob/1c88fa31e3685068b2ee42822e4efa22f49649bb/regions/core/compound.py#L147
There was a problem hiding this comment.
Go ahead with a simple circle or polygon, but the example workaround I gave for those more complex cases is something we should also highlight somewhere - but maybe not in the default docs.
|
@keflavich - found a polygon. The image is a bit larger and I tried to scale it down but it doesn't work for some reason. Also, I don't know if there's a simpler way (less LOC) to achieve the same result. |
setup.cfg
Outdated
| @@ -135,7 +135,7 @@ install_requires= | |||
| beautifulsoup4>=4.8 | |||
| html5lib>=0.999 | |||
| keyring>=15.0 | |||
| pyvo>=1.1 | |||
| pyvo>=1.5 | |||
There was a problem hiding this comment.
pyvo 1.5 isn't available yet - this is way too big a jump forward in requirements
There was a problem hiding this comment.
I'm pushing to have that PyVO released soon as it's about time to do. We need it for quantities in the enhanced table.
There was a problem hiding this comment.
this is a really big leap in required versions, especially if it's only required for alma.
(I suppose we'll need some of the SIA2 fixes in the 1.4.x series for irsa, too, but still 1.5 feels too strict a requirement)
Maybe, for now, do a version check in the module instead?
There was a problem hiding this comment.
Yeah, I feel strongly about this: a big version jump should be in-module only, not an install-level requirement
There was a problem hiding this comment.
there will be a lot of updates in pyvo that we'll rely on here (around download and handling cloud data), so at some point we'll need to bite the bullet
There was a problem hiding this comment.
I can easily isolate the 1.5 dependency in alma but it's just annoying that some users will have to install it manually (it's just the quantity feature).
1.4 has been out for a while and probably already used by the astroquery community. We haven't bumped the dependency here because the code did not depend on any new 1.4 API.
We obviously need to release 1.5 first and then can wait for a bit before merging this draft.
There was a problem hiding this comment.
Let's have a separate PR for updating the pyvo version requirement, and we can just make this PR depend on that one
There was a problem hiding this comment.
(we shouldn't be having discussions like this in a side comment on an unrelated PR... hard for people to find if it becomes important)
There was a problem hiding this comment.
I'll create an issue after 1.5 is out and we'll take it from there.
There was a problem hiding this comment.
let's aim to have a tagged 0.4.7 out, and then update dependencies. I want to drop some ancient python and astropy support as we don't need to drag those along indefinitely. I'm still on the fence of mandating to have a latest and greatest of something, so we may likely sit on requiring 1.5 pyvo for a while (but we can certainly upgrade the install instructions in a way that it would help users to upgrade the dependencies for them)
andamian
force-pushed
the
CADC-12747
branch
6 times, most recently
from
71344cb to
d841ca5
Compare
|
@keflavich, @bsipocz - I think this is finally ready. Footprints are a longstanding issue that is still being discussed in the IVOA circles. Until a general, protocol-supported solution is agreed upon (followed by an update of |
andamian
force-pushed
the
CADC-12747
branch
4 times, most recently
from
d23f14f to
4f87c18
Compare
There was a problem hiding this comment.
I have few minor comment about
- documentation (would like to have those before merging),
- namespace (we can/should discuss, and if it's a go ahead, can be done in a follow-up)
- nitpicking about import style, either way is fine by me, it was just strange to see the
import astroquery
docs/alma/alma.rst
Outdated
|
|
||
|
|
||
| .. doctest-remote-data:: | ||
| .. plot:: |
There was a problem hiding this comment.
Oh, I'm sorry, I didn't mention this before. The plot directive restarts the context, so it needs to be a completely standalone code snippet.
However, as with all the directives, empty lines and inserted narratives ends it, thus you run into the failures.
There was a problem hiding this comment.
It needs 2 directives (dockets-remote-data and plot) and that is not currently supported upstream.
#2687
@keflavich - this is a possible solution to this issue.
It turns out that because STC-S is not an IVOA recommendation and it's going to be replaced in the next version of the ObsCore spec, the fix could not be done upstream. It is instead a temporary ALMA solution that will be converted to the permanent one once the new spec is approved and the changes are propagated in the ALMA archive. There are a few more details to sort out, but at this stage I'm interested if this is an acceptable solution from the user perspective. I think that the enhanced results is preferred but I didn't make it default to ensure backwards compatibility. We could probably go further and return the position as coordinates as well if that's better than ra/dec and anything else that can be converted to Python objects.
Please let me know what you think.
Currently
produces: