Use `spherical_geometry` for accurate positioning

[martinjohndyer]

The SkyGrid.get_tile(coords) function is very questionable, and often returns the wrong answer.

For instance, with the GOTO-8 grid the position 114.48700 -25.94968 is clearly within T0314, but the basic get_tile function only checks for the nearest centre, meaning it returns the wrong responce:

In [1]: from gototile.grid import SkyGrid
   ...: grid = SkyGrid.from_name('GOTO-8')

In [2]: from astropy.coordinates import SkyCoord
   ...: coords = SkyCoord('114.48700 -25.94968', unit='deg')

In [3]: grid.get_tile(coords)
Out[3]: 'T0272'

In [4]: grid.get_tile(coords, overlap=True)
Out[4]: ['T0314']

Using overlap=True forces a different method, based on a base HealPIX grid to determine which points are within each tiles. That was correct in this instance, but can also fail. For example, if the point was within a pixel that's included in a neighbouring tile due to the inclusive=True in mhealpy's query_polygon. For instance the position 114.0 -26.2 slightly below that target is located in a pixel that's within both tiles. If we had to chose between them we'd still pick the closest centre, which clearly is still wrong.

In [6]: coords = SkyCoord('114 -26.2', unit='deg')

In [7]: grid.get_tile(coords)
Out[7]: 'T0272'

In [8]: grid.get_tile(coords, overlap=True)
Out[8]: ['T0272', 'T0314']

A far better method would be to use the spherical_geometry package https://github.com/spacetelescope/spherical_geometry, as I already point out in a comment.

goto-tile/gototile/grid.py

Lines 764 to 765 in c8ad3f3

	# TODO: spherical_geometry has better functions for things like overlaps, areas and
	# finding if a point is within an area

This would be perfectly simple to implement, for instance:

In [10]: from spherical_geometry import polygon, vector

In [11]: for tilename in ['T0314', 'T0272']:
    ...:     vertices = grid.vertices[grid.tilenames.index(tilename)]
    ...:     poly = polygon.SphericalPolygon.from_radec(vertices.ra, vertices.dec)
    ...:     target_vector = vector.radec_to_vector(coords.ra.deg, coords.dec.deg)
    ...:     print(tilename, poly.contains_point(target_vector))
T0314 True
T0272 False

Really the backend of the grid should be rewritten to use that package, with the tile polygons being precalculated with the grid. It would also give us better functions such as getting the area covered by each tile, or a test for overlaps.

Of course, you could still get cases where a point is is multiple tiles. Then if requested with overlap=False (the default) the best would be to return the case where the point was closer to the tile centre, which is what that function was supposed to do.

[martinjohndyer]

This is still planned, I'd really like it to be the big functional update once #101 is done with.

For now, here's a minimal example of finding which grid tile a particular coordinate is within. This should be a lot more reliable than the current grid.get_tile.

from astropy.coordinates import SkyCoord
from gototile.grid import SkyGrid
from spherical_geometry.polygon import SphericalPolygon
from spherical_geometry.vector import radec_to_vector

target = SkyCoord('114.48700 -25.94968', unit='deg')
target_vector = radec_to_vector(target.ra.deg, target.dec.deg)

grid = SkyGrid.from_name('GOTO')
tile_polys = [SphericalPolygon.from_radec(v.ra, v.dec) for v in grid.vertices]

for i in range(grid.ntiles):
    if tile_polys[i].contains_point(target_vector):
        print(grid.tilenames[i])

Ultimately of course the grid should define polygons on creation, and any functions like get_tile would convert to a vector internally.

[martinjohndyer]

Here's some still way too slow code to check which tiles overlap. No way we could run this at grid creation, even with some optimisations, but there could be some method on the grid to get overlapping tiles given a name.

import numpy as np
from gototile.grid import SkyGrid
from spherical_geometry.polygon import SphericalPolygon

grid = SkyGrid.from_name("GOTO")
tile_polys = [SphericalPolygon.from_radec(v.ra, v.dec) for v in grid.vertices]

overlaps = []
checked_overlaps = {}  # Dictionary to store already-checked overlaps

for i in range(grid.ntiles):
    poly = tile_polys[i]
    overlap_mask = [False] * grid.ntiles  # Initialize with False
    for j in range(i + 1, grid.ntiles):  # Only check pairs where j > i
        # Check if this pair has already been processed
        if (i, j) in checked_overlaps:
            overlap = checked_overlaps[(i, j)]
        else:
            overlap = poly.intersects_poly(tile_polys[j])
            # Store the result for both (i, j) and (j, i)
            checked_overlaps[(i, j)] = overlap
            checked_overlaps[(j, i)] = overlap
        overlap_mask[j] = overlap
    # Convert overlap_mask to include symmetric results
    for j in range(i):
        overlap_mask[j] = overlaps[j][i]
    print(grid.tilenames[i], np.array(grid.tilenames)[overlap_mask])
    overlaps.append(overlap_mask)