Repository: astropy/astroquery

Issue:
Issue with spatial constraints in Alma.query
The following script:

```python
import matplotlib.pyplot as plt
from astroquery.alma import Alma
from astropy.coordinates import SkyCoord
archive = Alma.query({'spatial_resolution': '0..0.5', 'galactic':'0..360 -2..2'})
coords = SkyCoord(archive['gal_longitude'], archive['gal_latitude'], frame='galactic')
plt.scatter(coords.l.deg, coords.b.deg)
plt.savefig('distribution.png')
```

produces the following plot:

![distribution](https://user-images.githubusercontent.com/314716/133170651-f96e9b3d-c71d-464f-98b5-2a8d60a2e105.png)

I would have expected all sources to lie between -2 and 2 in latitude, so it seems there is a bug somewhere along the way?


Discussion:
keflavich: Could you repost the original query here?  I think this looks better.
andamian: before

```select * from ivoa.obscore WHERE (0.0<=spatial_resolution AND spatial_resolution<=0.5) AND (INTERSECTS(RANGE_S2D(268.37542277754466,264.47377034405804,-29.96387952993674,-27.88035232035121), s_region) = 1) AND science_observation='T' AND data_rights='Public'```

now

```select * from ivoa.obscore WHERE (0.0<=spatial_resolution AND spatial_resolution<=0.5) AND (INTERSECTS(RANGE_S2D(264.47377034405804,268.37542277754466,-29.96387952993674,-27.88035232035121), s_region) = 1) AND science_observation='T' AND data_rights='Public'```
andamian: The unit test fails now because `SkyCoord` automatically makes ra=360 into 0 and that messes up the min/max in formula. Is there a way to prevent that?
keflavich: SkyCoord has a `wrap_at` attribute; set that to 180deg  - e.g. `coord.wrap_at(180*u.deg)`.

Yes, this is acceptable, but we should note the limitations, and particularly that searches aligned with the Galactic plane are not supported.  I think that's the right interpretation of what you just said
andamian: @keflavich - the crude translation between galactic range to ICRS range is a crude approximation that probably includes extra areas. The accuracy depends on the size of the range and position. It is however better than the bug that we have right now.

A better alternative would be to approximate the galactic range with a polygon and use that in the ICRS intersect. To do that, we'd need to add more points to the polygon along the latitude parallels to keep the polygon vertices follow the parallels (and not on the spherical arc). More points means more accurate but also more complexity in executing the query. Is this worth from a practical point of view? Are queries like the one in the ticket common enough to justify the effort or the current approximation is sufficient? 
keflavich: Let's use this workaround for now.  I think the polygon approach is much better, but I'm not convinced demand is high enough right now to justify the extra work.
andamian: Just an update here. I've tested the polygon approach with all transformation on the client side and it seems to work for the most part. An option to do the geometry transformation between different reference systems on the service side with a query function has been just endorsed by IVOA and we might follow that lead for this issue.
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/astropy/astroquery/pull/3105?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=astropy) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 67.39%. Comparing base [(`d347559`)](https://app.codecov.io/gh/astropy/astroquery/commit/d34755956829c5ebcd6916790670319bcd661aa0?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=astropy) to head [(`88b433b`)](https://app.codecov.io/gh/astropy/astroquery/commit/88b433b95c724c66ec322ab464b71e12f628c219?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=astropy).
> Report is 56 commits behind head on main.

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##             main    #3105      +/-   ##
==========================================
- Coverage   67.42%   67.39%   -0.03%     
==========================================
  Files         233      233              
  Lines       18396    18421      +25     
==========================================
+ Hits        12403    12415      +12     
- Misses       5993     6006      +13     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/astropy/astroquery/pull/3105?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=astropy).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=astropy).

andamian: ![Fixed](https://github.com/user-attachments/assets/ee2898a2-0f42-4306-8d59-f1442da104d5)

I think this will do it. @keflavich - please review.

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index aba8224a..79a1cb27 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -22,6 +22,8 @@ alma
 
 - Added support for frequency_resolution in KHz [#3035]
 
+- Changed the way galactic ranges are used in queries [#3105]
+
 ehst
 ^^^^
 
diff --git a/astroquery/alma/tapsql.py b/astroquery/alma/tapsql.py
index b6357baf..52f90655 100644
--- a/astroquery/alma/tapsql.py
+++ b/astroquery/alma/tapsql.py
@@ -69,16 +69,21 @@ def _gen_pos_sql(field, value):
                     dec_min = -90
                 if dec_max is None:
                     dec_max = 90
-                min_pt = coord.SkyCoord(ra_min, dec_min, unit=u.deg,
-                                        frame=frame)
-                max_pt = coord.SkyCoord(ra_max, dec_max, unit=u.deg,
-                                        frame=frame)
-                result += \
-                    "(INTERSECTS(RANGE_S2D({},{},{},{}), s_region) = 1)".\
-                    format(min_pt.icrs.ra.to(u.deg).value,
-                           max_pt.icrs.ra.to(u.deg).value,
-                           min_pt.icrs.dec.to(u.deg).value,
-                           max_pt.icrs.dec.to(u.deg).value)
+                ra_min = coord.Angle(ra_min, unit=u.degree).deg
+                ra_max = coord.Angle(ra_max, unit=u.degree).deg
+                dec_min = coord.Angle(dec_min, unit=u.degree).deg
+                dec_max = coord.Angle(dec_max, unit=u.degree).deg
+                if frame == 'galactic':
+                    # intersect with s_region is too complicated. ALMA indicated that
+                    # the use of gal_longitude and gal_latitude is good enough
+                    # approximation in this less common use case
+                    result += ('gal_longitude>={} AND gal_longitude<={} AND '
+                               'gal_latitude>={} AND gal_latitude<={}').format(
+                        ra_min, ra_max, dec_min, dec_max)
+                else:
+                    result += \
+                        "(INTERSECTS(RANGE_S2D({},{},{},{}), s_region) = 1)".\
+                        format(ra_min, ra_max, dec_min, dec_max)
             else:
                 raise ValueError('Cannot interpret ra({}), dec({}'.
                                  format(ra, dec))


Test patch:
diff --git a/astroquery/alma/tests/test_alma.py b/astroquery/alma/tests/test_alma.py
index 19f3fd01..a4a11106 100644
--- a/astroquery/alma/tests/test_alma.py
+++ b/astroquery/alma/tests/test_alma.py
@@ -111,7 +111,7 @@ def test_gen_pos_sql():
         "-90.0,4.338888888888889), s_region) = 1)"
     assert _gen_sql({'ra_dec': '!(10..20) >60'}) == common_select + \
         "((INTERSECTS(RANGE_S2D(0.0,10.0,60.0,90.0), s_region) = 1) OR " \
-        "(INTERSECTS(RANGE_S2D(20.0,0.0,60.0,90.0), s_region) = 1))"
+        "(INTERSECTS(RANGE_S2D(20.0,360.0,60.0,90.0), s_region) = 1))"
     assert _gen_sql({'ra_dec': '0..20|40..60 <-50|>50'}) == common_select + \
         "((INTERSECTS(RANGE_S2D(0.0,20.0,-90.0,-50.0), s_region) = 1) OR " \
         "(INTERSECTS(RANGE_S2D(0.0,20.0,50.0,90.0), s_region) = 1) OR " \
@@ -123,18 +123,23 @@ def test_gen_pos_sql():
     assert _gen_sql({'galactic': '1 2, 3'}) == common_select + "(INTERSECTS(" \
         "CIRCLE('ICRS',{},{},3.0), s_region) = 1)".format(
         center.icrs.ra.to(u.deg).value, center.icrs.dec.to(u.deg).value)
-    min_point = coord.SkyCoord('12:13:14.0', '-00:01:02.1', unit=u.deg,
-                               frame='galactic')
-    max_point = coord.SkyCoord('12:14:14.0', '-00:00:02.1', unit=(u.deg, u.deg),
-                               frame='galactic')
+    gal_longitude = ('12:13:14.0', '12:14:14.0')
+    gal_latitude = ('-00:01:02.1', '-00:00:02.1')
+    min_pt = coord.SkyCoord(gal_longitude[0], gal_latitude[0], unit=u.deg)
+    long_min, lat_min = min_pt.ra.value, min_pt.dec.value
+    max_pt = coord.SkyCoord(gal_longitude[1], gal_latitude[1], unit=u.deg)
+    long_max, lat_max = max_pt.ra.value, max_pt.dec.value
     assert _gen_sql(
-        {'galactic': '12:13:14.0..12:14:14.0 -00:01:02.1..-00:00:02.1'}) == \
+        {'galactic': '{}..{} {}..{}'.format(
+            gal_longitude[0], gal_longitude[1], gal_latitude[0], gal_latitude[1])}) == \
         common_select +\
-        "(INTERSECTS(RANGE_S2D({},{},{},{}), s_region) = 1)".format(
-            min_point.icrs.ra.to(u.deg).value,
-            max_point.icrs.ra.to(u.deg).value,
-            min_point.icrs.dec.to(u.deg).value,
-            max_point.icrs.dec.to(u.deg).value)
+        'gal_longitude>={} AND gal_longitude<={} AND gal_latitude>={} AND gal_latitude<={}'.format(
+            long_min, long_max, lat_min, lat_max)
+
+    # test extremities
+    assert _gen_sql({'galactic': '0..360 -90..90'}) == \
+        common_select + ('gal_longitude>=0.0 AND gal_longitude<=360.0 AND '
+                         'gal_latitude>=-90.0 AND gal_latitude<=90.0')
 
     # combination of frames
     center = coord.SkyCoord(1, 2, unit=u.deg, frame='galactic')