Enhanced results feature in ALMA

[andamian]

#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

print(Alma.query_object('M42', public=True, science=True)[0]['s_region'])
print(Alma(enhanced_results=True).query_object('M42', public=True, science=True)[0]['s_region'])

produces:

Polygon ICRS 83.831185 -5.264207 83.827356 -5.260845 83.830732 -5.257032 83.834561 -5.260394
WARNING: column frequency_support has a unit but is kept as a MaskedColumn as an attempt to convert it to Quantity failed with:
TypeError('The value must be a valid Python or Numpy numeric type.') [astropy.table.table]
Region: PolygonSkyRegion
vertices: <SkyCoord (ICRS): (ra, dec) in deg
    [(83.831185, -5.264207), (83.827356, -5.260845),
     (83.830732, -5.257032), (83.834561, -5.260394)]>

[pep8speaks]

Hello @andamian! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

There are currently no PEP 8 issues detected in this Pull Request. Cheers! :beers:

Comment last updated at 2024-06-10 21:01:26 UTC

[codecov[bot]]

Codecov Report

Attention: Patch coverage is 93.42105% with 5 lines in your changes missing coverage. Please review.

Project coverage is 67.24%. Comparing base (462c2da) to head (c62f9de). Report is 16 commits behind head on main.

:exclamation: Current head c62f9de differs from pull request most recent head 434c4ec

Please upload reports for the commit 434c4ec to get more accurate results.

Files	Patch %	Lines
astroquery/alma/core.py	93.33%	5 Missing :warning:

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     

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[keflavich]

Yes, that looks like a good solution.

[andamian]

@keflavich besides Polygon and Circle, ALMA has compose shapes such as Union (of Polygons). I've solved that by returning regions.Regions (essentially a list of shapes) - I hope that's the correct approach. While running the parser over the entire ALMA collection I came across the Not operator. What kind of regions is this, for example and how can this be translated to regions.Regions?

Polygon ICRS 197.872406 -1.363399 197.864495 -1.361505 197.859768 -1.353969 197.857210 -1.345208 197.857042 -1.339734 197.858840 -1.337233 197.855764 -1.334732 197.855568 -1.331622 197.862544 -1.327086 197.864642 -1.323653 197.869721 -1.324977 197.871969 -1.320892 197.876753 -1.323038 197.886090 -1.318814 197.891460 -1.318334 197.894748 -1.324836 197.894074 -1.337986 197.888867 -1.341784 197.888477 -1.347778 197.883592 -1.359877 Not (Polygon 197.865541 -1.335028 197.866812 -1.336468 197.871866 -1.334542 197.878114 -1.335570 197.881730 -1.340990 197.884624 -1.332046 197.876248 -1.330377) Not (Polygon 197.871726 -1.341833 197.870114 -1.342896 197.874558 -1.343510 197.873973 -1.341334) Not (Polygon 197.875258 -1.349535 197.869529 -1.350188 197.866855 -1.348762 197.865919 -1.350032 197.869432 -1.356172 197.872312 -1.354448 197.875506 -1.355953 197.878824 -1.354446 197.880479 -1.348022)

[keflavich]

Yes, just stitching them together in a list is the way to go.

Regions do have built-in exclude/include logic: the .meta['include'] attribute is 1 for include, 0 for exclude. ds9 regions serialize this as a leading - sign or a trailing # include = 0. CRTF serializes as leading -.

[keflavich]

You can also consider making CompositeSkyRegions, e.g.: regions.CompoundSkyRegion(reg1, reg2, operator.or) ...not entirely sure what operator is and not...

[andamian]

Thanks. I've somehow missed the "Combining Regions" section. I'll use CompoundSkyRegion for the ALMA union, but I still don't know what the above regions represents. We meet with ALMA tomorrow and I'll request clarifications.

[bsipocz]

Side note: what about supporting MOC?

[andamian]

@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 A|B (OR), A&B (AND) and A^B(XOR).

[keflavich]

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.

[andamian]

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 https://github.com/astropy/astroquery/pull/2855#issuecomment-1777575895, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABOTQLWDBHNYO6YMRCWHKKLYA7SQVAVCNFSM6AAAAAA6GJ3ALSVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTONZXGU3TKOBZGU . You are receiving this because you were mentioned.Message ID: @.***>

[andamian]

@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 (UNION and NOT)? They result in PolygonSkyRegion and CompoundSkyRegion). How can one tell if they are the correct shapes? The test examples are real one.

[keflavich]

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)

[andamian]

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)

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

[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')).celestial

that WCS will work

[andamian]

@keflavich - making little progress with the wcs and ALMA compound regions , however while playing with it I've run into https://github.com/astropy/regions/pull/470. The corresponding PR only deals with the missing centre attribute, however the issue is that other shape operations such as UNION are still not supported. It would be nice to use something that's working if we are to propose this feature to the ALMA user community.

[keflavich]

I agree that plotting is an important feature - I think we can push to get https://github.com/astropy/regions/pull/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?

[andamian]

This case for example (https://almascience.eso.org/aq/?observationsProjectCode=2021.1.00547.S) with this shape: Polygon ICRS 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 Not (Polygon 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)

This is the a - b case. If I do it with a xor b it works because b is contained in a, but if I do it with the generic (a or b) xor b it fails on the or operation.

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.

[keflavich]

For my sake, I'm including graphics here: is the first one (https://almascience.eso.org/aq/?observationsProjectCode=2021.1.00547.S)

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).

but if I do it with the generic (a or b) xor b it fails on the or operation.

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.

[andamian]

The problem is here: https://github.com/astropy/regions/blob/1c88fa31e3685068b2ee42822e4efa22f49649bb/regions/core/compound.py#L147

Union is not expected in a compound region.

[keflavich]

Right; I think what we should be doing here is either: (A) plot all regions as individual artists and, maybe, color the not-include regions as different from the include (B) create an include mask and use a contour to display it Both of these feel somewhat out-of-scope for astroquery and are better redirected to regions, but I recommend we include a demo of how to do (B) in this PR and call it done. That's the only approach that is straightforward.

[keflavich]

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) ^ poly2

One 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]);

but done as comp.plot(...) instead of manually

[andamian]

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 PyVO release to get quantities as well as regions. Does the testing look OK? Anything else missing? Didn't want to add any info about plotting in the documentation although those footprints would catch the eye.

[andamian]

@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.

[andamian]

@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 ObsCore standard and propagation to ALMA services - i.e. years away), this should work for our ALMA users.