'SCRIPT_NUMBER_ID' column disappeared in Simbad.query_region() after v0.4.8 update

[vsquicciarini]

Hello there,
it appears that the v0.4.8 of astroquery deleted the column 'SCRIPT_NUMBER_ID' when using the function Simbad.query_region(). The column, introduced here (#2047) was crucial to connect the unequivocally inputs and outputs of a query, by assigning each row a number corresponding to the i-th input region.

Apparently the column was renamed as 'object_number_id' in Simbad.query_objects() (#2954). However, no such column exists for Simbad.query_region().

Would it be possible to restore this column (even with a different name)?

Thank you in advance for your the consideration. :)

[keflavich]

@ManonMarchand I expect you'll have something to say on this?

(@vsquicciarini you should expect to wait ~a week to hear back on this at least, as the relevant folks are on holiday right now)

[ManonMarchand]

Indeed, there is now a single script even with multiple radii.

The behaviour could be mimicked just as in query_objects. What would be sensible, keep the historical name script_number_id or something closer to reality like center_number or center_number_id?

Edit: Might not be possible actually, cannot use an upload like in query_objects. Thinking about how to make it work, but a workaround is to cross match the output with the entry. Something using astropy's search_around_sky maybe. In that case, that'd be a documentation tip to explain how to interpret the results of query region:

from astropy.coordinates import SkyCoord
from astropy.table import Table, hstack
import astropy.units as u
from astroquery.simbad import Simbad
# let's do a query for two centers
simbad = Simbad()
centers = Table({"ra": [0, 1], "dec": [0, 1], "cone_number": [0, 1]})
coo_centers = SkyCoord(centers["ra"], centers["dec"], unit="deg")
regions = simbad.query_region(coo_centers, radius="0.1d")
coo_result = SkyCoord(regions["ra"], regions["dec"], unit="deg")
# here, we do the cross match
idx1, idx2, _, _ = coo_centers.search_around_sky(coo_result, 0.1*u.deg)
hstack([regions[idx1], centers[["cone_number"]][idx2]])

<Table length=51>
          main_id                     ra         ... cone_number
                                     deg         ...
           object                  float64       ...    int32
---------------------------- ------------------- ... -----------
                    PB  5669       0.04984129467 ...           0
              QSO B2357-003A     359.99988315297 ...           0
    SDSS J000351.87+010511.4            0.966156 ...           1
    SDSS J000421.84+005743.6       1.09105216744 ...           1
            SDSS-II SN 14480 0.01745833333333333 ...           0
Gaia DR3 2738341959609844992  0.9648817681758335 ...           1
            UCAC4 450-138351  359.91676860000007 ...           0
                    Z 382-15       0.03259767004 ...           0
    SDSS J000342.95+010109.2  0.9290047483600001 ...           1
                         ...                 ... ...         ...
     2MASX J00034546+0102404  0.9394874682400001 ...           1
               SDSSCGB 49921  1.0616666666666668 ...           1
            IRAS F23572-0021   359.9485833333334 ...           0
            SDSS-II SN 15730  359.92387916666667 ...           0
    SDSS J235951.01+000453.9  359.96257883642255 ...           0
Gaia DR2 2738345636100534016  0.9106980628416665 ...           1
  [LVO2003] J000355.2+005519  0.9800000000000001 ...           1
            UCAC4 455-000077  0.9549125032562502 ...           1
            UCAC4 455-000082  1.0003119072266666 ...           1
               MCG+00-01-014       0.01849340605 ...           0

[vsquicciarini]

@ManonMarchand thanks a lot for your help (and apologies for the late answer). The proposed solution would work, provided that the queried regions, as defined from coordinates and the search radius, do not overlap. Unfortunately, if some regions overlap, the cone_number would pick for each entry the closest region center and not the original region ID. An example to clarify this point:

Here I have used on purpose the old astroquery v0.4.7 to compare 'SCRIPT_NUMBER_ID' and 'cone_number'.
Do you have thoughts on this point?

[ManonMarchand]

No worries 🙂

The cross match is not not giving you the closest match, it's giving you dupplicates since the points are in both regions. You have the same dupplicate in the case of script number ID.

Here in this example each source that belongs in the two regions appears four times. Two from the script ID and two from the cone number. If you apply only one of the two methods you should have the same outputs.

Maybe a way to make it easier to see is to sort by cone_number and compare to the version you've got in the screenshot.

Edit: there is something strange about the quasar though. It's showing up in the two cones for script_id and only in one with astropy's cross-match. Maybe a precision issue to dig into between the two methods?

Edit2: It is not a precision issue, you just used 2*u.arcmin separation in the cone_id case while the angle you chose for query_radius is 0.1d. Hence the difference in the output when the separations are respectively 0.034d and 0.033d between the quasar and the two centers 😆
All good here, the two methods give the same output 🙂

[vsquicciarini]

Hi again, thanks for the clarification! Indeed, if I define a pseudo-index (main_id+cone_number) so as to remove the duplication given by the cone number, and then I sort by cone number, I get the same result. So the issue appears to be solved :)
Do you think the column has any chance to be reinserted in the next version of the package using this workaround? I now know thanks to you how to write a snippet of code for all my scripts, but everyone else using the old column will have to independently find out a solution.
Thanks again for your help!

[ManonMarchand]

Hi!
Just to keep you informed, I tried adding this to the method, but search_around_sky does not work when the radii are different. I'm looking for a solution that would work in every case