BUG: Square search queries actually produce rectangular regions of sky

[cam92473]

I've tested 3 of astroquery's Gaia catalogue query methods: Gaia.launch_job_async, Gaia.query_object_async and Gaia.cone_search_async. For Gaia.launch_job_async I used two ADQL queries, one for a circular region and one for a square region (rectangular region with both sidelengths equal). The circular queries return sources in the correct shape, that is, the space the sources occupy is elliptical in right ascension and declination, but circular when the right ascension and declination are projected (say using the Gnomonic projection). However, the square queries return sources that occupy a square region in right ascension and declination, and when these sources are projected, the region is rectangular.

Here is the test code for Gaia.launch_job_async with a square query:

from astroquery.gaia import Gaia
import astropy.units as u
from astropy.coordinates import SkyCoord
import matplotlib.pyplot as plt
import numpy as np
from math import sin, cos, pi

center_ra = 80
center_dec = -70
search_length = .01

#JOBSTR BOX
jobstr = "SELECT TOP 1000 gaia_source.ra,gaia_source.dec FROM gaiaedr3.gaia_source\n"
jobstr += "WHERE 1=CONTAINS(POINT('ICRS', gaiaedr3.gaia_source.ra,gaiaedr3.gaia_source.dec),"
jobstr += "BOX('ICRS',{},{},{},{}))\n".format(center_ra,center_dec,search_length*2,search_length*2)
job = Gaia.launch_job_async(jobstr, dump_to_file=False, verbose=False)
results = job.get_results()
print(results)

ras = np.array(results['ra'])
decs = np.array(results['dec'])

xs = np.zeros(152)
ys = np.zeros(152)

center_ra_rad = center_ra*pi/180
center_dec_rad = center_dec*pi/180

for i in range(152):
    a = ras[i]*pi/180
    d = decs[i]*pi/180
    X = cos(d)*sin(a-center_ra_rad) / (cos(center_dec_rad)*cos(d)*cos(a-center_ra_rad)+sin(center_dec_rad)*sin(d))
    Y = (cos(center_dec_rad)*sin(d) - cos(d)*sin(center_dec_rad)*cos(a-center_ra_rad)) / (cos(center_dec_rad)*cos(d)*cos(a-center_ra_rad) + sin(center_dec_rad)*sin(d))
    xs[i] = X
    ys[i] = Y

fig = plt.figure(1,figsize=(16,8))
ax1 = fig.add_subplot(121)
ax1.scatter(ras,decs,marker='.')
ax1.axis('equal')
ax2 = fig.add_subplot(122)
ax2.scatter(xs,ys,marker='.')
ax2.axis('equal')

plt.show()
plt.close()

Here is the test code for Gaia.launch_job_async with a circular query:

from astroquery.gaia import Gaia
import astropy.units as u
from astropy.coordinates import SkyCoord
import matplotlib.pyplot as plt
import numpy as np
from math import sin, cos, pi

center_ra = 80
center_dec = -70
search_length = .01


#JOBSTR CIRCLE
jobstr = "SELECT TOP 1000 gaia_source.ra,gaia_source.dec FROM gaiaedr3.gaia_source\n"
jobstr += "WHERE 1=CONTAINS(POINT('ICRS', gaiaedr3.gaia_source.ra,gaiaedr3.gaia_source.dec),"
jobstr += "CIRCLE('ICRS',{},{},{}))\n".format(center_ra,center_dec,search_length)
job = Gaia.launch_job_async(jobstr, dump_to_file=False, verbose=False)
results = job.get_results()
print(results)

ras = np.array(results['ra'])
decs = np.array(results['dec'])
 
xs = np.zeros(352)
ys = np.zeros(352)

center_ra_rad = center_ra*pi/180
center_dec_rad = center_dec*pi/180

for i in range(352):
    a = ras[i]*pi/180
    d = decs[i]*pi/180
    X = cos(d)*sin(a-center_ra_rad) / (cos(center_dec_rad)*cos(d)*cos(a-center_ra_rad)+sin(center_dec_rad)*sin(d))
    Y = (cos(center_dec_rad)*sin(d) - cos(d)*sin(center_dec_rad)*cos(a-center_ra_rad)) / (cos(center_dec_rad)*cos(d)*cos(a-center_ra_rad) + sin(center_dec_rad)*sin(d))
    xs[i] = X
    ys[i] = Y

fig = plt.figure(2,figsize=(16,8))
ax1 = fig.add_subplot(121)
ax1.scatter(ras,decs,marker='.')
ax1.axis('equal')
ax2 = fig.add_subplot(122)
ax2.scatter(xs,ys,marker='.')
ax2.axis('equal')

plt.show()
plt.close()

Here is the test code for Gaia.query_object_async (square):

from astroquery.gaia import Gaia
import astropy.units as u
from astropy.coordinates import SkyCoord
import matplotlib.pyplot as plt
import numpy as np
from math import sin, cos, pi

center_ra = 80
center_dec = -70
search_length = .01

#QUERY OBJECT SQUARE
Gaia.ROW_LIMIT = 1000
Gaia.MAIN_GAIA_TABLE = 'gaiaedr3.gaia_source'
coord = SkyCoord(ra=center_ra,dec=center_dec,unit=(u.degree,u.degree),frame='icrs')
width = u.Quantity(search_length*5,u.deg)
height = u.Quantity(search_length*5,u.deg)
r = Gaia.query_object_async(coordinate=coord, width=width, height=height, columns=['ra','dec'])
print(r)

ras = np.array(r['ra'])
decs = np.array(r['dec'])

xs = np.zeros(1000)
ys = np.zeros(1000)

center_ra_rad = center_ra*pi/180
center_dec_rad = center_dec*pi/180

for i in range(500):
    a = ras[i]*pi/180
    d = decs[i]*pi/180
    X = cos(d)*sin(a-center_ra_rad) / (cos(center_dec_rad)*cos(d)*cos(a-center_ra_rad)+sin(center_dec_rad)*sin(d))
    Y = (cos(center_dec_rad)*sin(d) - cos(d)*sin(center_dec_rad)*cos(a-center_ra_rad)) / (cos(center_dec_rad)*cos(d)*cos(a-center_ra_rad) + sin(center_dec_rad)*sin(d))
    xs[i] = X
    ys[i] = Y

fig = plt.figure(3,figsize=(16,8))
ax1 = fig.add_subplot(121)
ax1.scatter(ras,decs,marker='.')
ax1.axis('equal')
ax2 = fig.add_subplot(122)
ax2.scatter(xs,ys,marker='.')
ax2.axis('equal')

plt.show()
plt.close()

And lastly here is the test code for Gaia.cone_search_async:

from astroquery.gaia import Gaia
import astropy.units as u
from astropy.coordinates import SkyCoord
import matplotlib.pyplot as plt
import numpy as np
from math import sin, cos, pi

center_ra = 80
center_dec = -70
search_length = .01

#CONE SEARCH
Gaia.ROW_LIMIT = 1000
Gaia.MAIN_GAIA_TABLE = 'gaiaedr3.gaia_source'
coord = SkyCoord(ra=center_ra,dec=center_dec,unit=(u.degree,u.degree),frame='icrs')
radius = u.Quantity(search_length,u.deg)
j = Gaia.cone_search_async(coord, radius, columns=['ra','dec'])
r = j.get_results()
print(r)
    
ras = np.array(r['ra'])
decs = np.array(r['dec'])

xs = np.zeros(352)
ys = np.zeros(352)

center_ra_rad = center_ra*pi/180
center_dec_rad = center_dec*pi/180

for i in range(352):
    a = ras[i]*pi/180
    d = decs[i]*pi/180
    X = cos(d)*sin(a-center_ra_rad) / (cos(center_dec_rad)*cos(d)*cos(a-center_ra_rad)+sin(center_dec_rad)*sin(d))
    Y = (cos(center_dec_rad)*sin(d) - cos(d)*sin(center_dec_rad)*cos(a-center_ra_rad)) / (cos(center_dec_rad)*cos(d)*cos(a-center_ra_rad) + sin(center_dec_rad)*sin(d))
    xs[i] = X
    ys[i] = Y

fig = plt.figure(4,figsize=(16,8))
ax1 = fig.add_subplot(121)
ax1.scatter(ras,decs,marker='.')
ax1.axis('equal')
ax2 = fig.add_subplot(122)
ax2.scatter(xs,ys,marker='.')
ax2.axis('equal')

plt.show()
plt.close()

Thanks.

[eerovaher]

Doesn't the box search in the Gaia archive web interface work the same way? If the behavior of astroquery matches that of the web interface then there is no bug, though there might be an issue of insufficient documentation.

[cam92473]

Good point. Astroquery does match the behaviour of the web interface, so there is no bug. I'll check with those at the Gaia archive to see what they say about the issue.