Remove esutil dependency as it causes unpredictable errors on installation,...

Remove esutil dependency as it causes unpredictable errors on installation, and has only been used for spherical matching. Now we have a drop-in replacement for that.

doc/installation.rst

@@ -51,7 +51,7 @@ Install conda-forge dependencies
 
 .. prompt:: bash
 
-   conda install -c conda-forge esutil astroscrappy
+   conda install -c conda-forge astroscrappy
 
 Install conda astropy dependencies
 

doc/requirements_rtd.txt

@@ -16,7 +16,6 @@ astroquery>0.4.1
 sep
 photutils
 statsmodels
-esutil==0.6.0
 tqdm
 psycopg2
 regions

setup.py

@@ -16,7 +16,6 @@ requirements = [
     'sep',
     'photutils',
     'statsmodels',
-    'esutil',
     'tqdm',
     'regions'
 ]

stdpipe/astrometry.py

@@ -3,13 +3,12 @@ from __future__ import absolute_import, division, print_function, unicode_litera
 import os, tempfile, shutil, shlex, re, warnings
 import numpy as np
 
-from esutil import coords, htm
-
 from astropy.wcs import WCS
 from astropy.io import fits
 from astropy.wcs.utils import fit_wcs_from_points
-from astropy.coordinates import SkyCoord
+from astropy.coordinates import SkyCoord, search_around_sky
 from astropy.table import Table
+from astropy import units as u
 
 from scipy.stats import chi2
 
@@ -34,10 +33,10 @@ def get_frame_center(filename=None, header=None, wcs=None, width=None, height=No
         elif shape is not None:
             height,width = shape
 
-    [ra1],[dec1] = wcs.all_pix2world([0], [0], 1)
-    [ra0],[dec0] = wcs.all_pix2world([width/2], [height/2], 1)
+    ra1,dec1 = wcs.all_pix2world(0, 0, 0)
+    ra0,dec0 = wcs.all_pix2world(width/2, height/2, 0)
 
-    sr = coords.sphdist(ra0, dec0, ra1, dec1)[0]
+    sr = spherical_distance(ra0, dec0, ra1, dec1)
 
     return ra0, dec0, sr
 
@@ -71,6 +70,47 @@ def xyztoradec(xyz):
 
     return (np.rad2deg(ra), np.rad2deg(dec))
 
+def spherical_distance(ra1, dec1, ra2, dec2):
+    """Spherical distance.
+
+    :param ra1: First point or set of points RA
+    :param dec1: First point or set of points Dec
+    :param ra2: Second point or set of points RA
+    :param dec2: Second point or set of points Dec
+    :returns: Spherical distance in degrees
+
+    """
+
+    x = np.sin(np.deg2rad((ra1 - ra2)/2))
+    x *= x;
+    y = np.sin(np.deg2rad((dec1 - dec2)/2))
+    y *= y;
+
+    z = np.cos(np.deg2rad((dec1 + dec2)/2))
+    z *= z;
+
+    return np.rad2deg(2*np.arcsin(np.sqrt(x*(z - y) + y)))
+
+def spherical_match(ra1, dec1, ra2, dec2, sr=1/3600):
+    """Positional match on the sphere for two lists of coordinates.
+
+    Aimed to be a direct replacement for :func:`esutil.htm.HTM.match` method with :code:`maxmatch=0`.
+
+    :param ra1: First set of points RA
+    :param dec1: First set of points Dec
+    :param ra2: Second set of points RA
+    :param dec2: Second set of points Dec
+    :param sr: Maximal acceptable pair distance to be considered a match, in degrees
+    :returns: Two parallel sets of indices corresponding to matches from first and second lists, along with the pairwise distances in degrees
+
+    """
+
+    idx1,idx2,dist,_ = search_around_sky(SkyCoord(ra1, dec1, unit='deg'), SkyCoord(ra2, dec2, unit='deg'), sr*u.deg)
+
+    dist = dist.deg # convert to degrees
+
+    return idx1, idx2, dist
+
 def get_objects_center(obj, col_ra='ra', col_dec='dec'):
     """
     Returns the center RA, Dec, and radius in degrees for a cloud of objects on the sky.
@@ -79,7 +119,7 @@ def get_objects_center(obj, col_ra='ra', col_dec='dec'):
     xyz0 = np.mean(xyz, axis=1)
     ra0,dec0 = xyztoradec(xyz0)
 
-    sr0 = np.max(coords.sphdist(ra0, dec0, obj[col_ra], obj[col_dec]))
+    sr0 = np.max(spherical_distance(ra0, dec0, obj[col_ra], obj[col_dec]))
 
     return ra0, dec0, sr0
 
@@ -337,8 +377,7 @@ def refine_wcs(obj, cat, order=2, match=True, sr=3/3600, update=False,
 
     if match:
         # Perform simple nearest-neighbor matching within given radius
-        h = htm.HTM(10)
-        oidx,cidx,dist = h.match(obj['ra'], obj['dec'], cat[cat_col_ra], cat[cat_col_dec], sr, maxmatch=0)
+        oidx,cidx,dist = spherical_match(obj['ra'], obj['dec'], cat[cat_col_ra], cat[cat_col_dec], sr)
         _obj = obj[oidx]
         _cat = cat[cidx]
     else:

stdpipe/catalogs.py

@@ -16,8 +16,6 @@ from astropy.coordinates import SkyCoord
 from astropy import units as u
 from astropy.time import Time
 
-from esutil import htm
-
 from . import astrometry
 
 catalogs = {
@@ -199,8 +197,7 @@ def xmatch_skybot(obj, sr=10/3600, time=None, col_ra='ra', col_dec='dec', col_id
         return None
 
     # Cross-match objects
-    h = htm.HTM(10)
-    oidx,cidx,dist = h.match(obj[col_ra], obj[col_dec], xcat['RA'], xcat['DEC'], 10/3600)
+    oidx,cidx,dist = astrometry.spherical_match(obj[col_ra], obj[col_dec], xcat['RA'], xcat['DEC'], 10/3600)
 
     # Annotate the table with id from objects so that it is possible to identify the matches
     xcat[col_id] = MaskedColumn(len(xcat), dtype=np.dtype(obj[col_id][0]))

stdpipe/photometry.py

@@ -22,10 +22,10 @@ import photutils
 from photutils.utils import calc_total_error
 
 import statsmodels.api as sm
-from esutil import htm
-
 from scipy.optimize import minimize
 
+from . import astrometry
+
 try:
     import cv2
     # Much faster dilation
@@ -544,9 +544,7 @@ def match(obj_ra, obj_dec, obj_mag, obj_magerr, obj_flags, cat_ra, cat_dec, cat_
     # Simple wrapper around print for logging in verbose mode only
     log = (verbose if callable(verbose) else print) if verbose else lambda *args,**kwargs: None
 
-    h = htm.HTM(10)
-
-    oidx,cidx,dist = h.match(obj_ra, obj_dec, cat_ra, cat_dec, sr, maxmatch=0)
+    oidx,cidx,dist = astrometry.spherical_match(obj_ra, obj_dec, cat_ra, cat_dec, sr)
 
     log(len(dist), 'initial matches between', len(obj_ra), 'objects and', len(cat_ra), 'catalogue stars, sr =', sr*3600, 'arcsec')
     log('Median separation is %.2f arcsec' % (np.median(dist)*3600))

stdpipe/pipeline.py

@@ -15,8 +15,6 @@ from astropy.coordinates import SkyCoord
 from astropy.time import Time
 from astropy.table import Table
 
-from esutil import htm
-
 from . import photometry
 from . import astrometry
 from . import catalogs
@@ -165,8 +163,6 @@ def filter_transient_candidates(obj, sr=None, pixscale=None, time=None,
         # So let's create a copy of objects list where we may freely add extra columns without touching the original
         obj_in = obj_in.copy()
 
-    h = htm.HTM(10)
-
     log('Candidate filtering routine started with %d initial candidates and %.1f arcsec matching radius' % (len(obj), sr*3600))
     cand_idx = np.ones(len(obj), dtype=np.bool)
 
@@ -184,7 +180,7 @@ def filter_transient_candidates(obj, sr=None, pixscale=None, time=None,
     if cat is not None and remove == False:
         obj_in['candidate_refcat'] = False
     if cat is not None and np.any(cand_idx):
-        m = h.match(obj['ra'], obj['dec'], cat[cat_col_ra], cat[cat_col_dec], sr)
+        m = astrometry.spherical_match(obj['ra'], obj['dec'], cat[cat_col_ra], cat[cat_col_dec], sr)
         cand_idx[m[0]] = False
 
         if remove == False:

stdpipe/templates.py

@@ -16,8 +16,6 @@ from astropy.io import fits
 from astropy.stats import mad_std
 from astropy.table import Table
 
-from esutil import htm
-
 # from scipy.ndimage import binary_dilation
 from astropy.convolution import Tophat2DKernel, convolve, convolve_fft
 
@@ -296,8 +294,7 @@ def find_ps1_skycells(ra, dec, sr, band='r', ext='image', cell_radius=0.3, fullp
         __skycells = Table.read(utils.get_data_path('ps1skycells.txt'), format='ascii')
 
     # FIXME: here we may select the cells that are too far from actual footprint
-    h = htm.HTM(10)
-    _,idx,_ = h.match(ra, dec, __skycells['ra0'], __skycells['dec0'], sr + cell_radius, maxmatch=0)
+    _,idx,_ = astrometry.spherical_match(ra, dec, __skycells['ra0'], __skycells['dec0'], sr + cell_radius)
 
     if fullpath:
         # Get full path on the server