Input (lon, lat) in degrees instead of (theta, phi) in radians.
    Returns HEALPix map at the desired resolution
500541
def healpixMap(nside, lon, lat, fill_value=0., nest=False):
    

    lon_median, lat_median = np.median(lon), np.median(lat)
    max_angsep = np.max(ugali.utils.projector.angsep(lon, lat, lon_median, lat_median))
    
    pix = angToPix(nside, lon, lat, nest=nest)
    if max_angsep < 10:
        # More efficient histograming for small regions of sky
        m = np.tile(fill_value, healpy.nside2npix(nside))
        pix_subset = ugali.utils.healpix.angToDisc(nside, lon_median, lat_median, max_angsep, nest=nest)
        bins = np.arange(np.min(pix_subset), np.max(pix_subset) + 1)
        m_subset = np.histogram(pix, bins=bins - 0.5)[0].astype(float)
        m[bins[0:-1]] = m_subset
    else:
        m = np.histogram(pix, np.arange(hp.nside2npix(nside) + 1))[0].astype(float)
    if fill_value != 0.:
        m[m == 0.] = fill_value
    return m
572333
def lon_lat(value):
    
    lon, lat = value.split()
    return longitude(lon), latitude(lat)
1180275
def pfLon(jd, lat, lon):
    
    sun = swe.sweObjectLon(const.SUN, jd)
    moon = swe.sweObjectLon(const.MOON, jd)
    asc = swe.sweHousesLon(jd, lat, lon,
                           const.HOUSES_DEFAULT)[1][0]
    
    if isDiurnal(jd, lat, lon):
        return angle.norm(asc + moon - sun)
    else:
        return angle.norm(asc + sun - moon)
819716
def get_sky_coords(self):
        
        if self._ipix is None:
            theta, phi = hp.pix2ang(
                self._nside, list(range(self._npix)), self._nest)
        else:
            theta, phi = hp.pix2ang(self._nside, self._ipix, self._nest)

        lat = np.degrees((np.pi / 2) - theta)
        lon = np.degrees(phi)
        return np.vstack([lon, lat]).T
847936
function getDistanceInKilometers(p1, p2) {
  // Credit: adapts code for Haversine formula and degrees-to-radian conversion found at:
  // http://www.movable-type.co.uk/scripts/latlong.html

  var R = 6371; // earth's radius in kilometers
  var dLat = degreesToRadians(p2.lat - p1.lat);
  var dLon = degreesToRadians(p2.lng - p1.lng); 
  var lat1 = degreesToRadians(p1.lat);
  var lat2 = degreesToRadians(p2.lat);

  var a = 
    Math.sin(dLat/2) * Math.sin(dLat/2)
    + Math.sin(dLon/2) * Math.sin(dLon/2)
    * Math.cos(lat1) * Math.cos(lat2);

  var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a)); 
  var distance = R * c; // kilomter
  return distance;
}
161689
def pixel_coords(self, latlon, reverse=False):
        '''
        '''
        state = self.state
        if reverse:
            (x,y) = latlon
            return self.coordinates(x,y)
        (lat,lon) = (latlon[0], latlon[1])
        return state.mt.coord_to_pixel(state.lat, state.lon, state.width, state.ground_width, lat, lon)
102417
def ToLatLng(self):
    
    rad_lat = math.atan2(self.z, math.sqrt(self.x * self.x + self.y * self.y))
    rad_lng = math.atan2(self.y, self.x)
    return (rad_lat * 180.0 / math.pi, rad_lng * 180.0 / math.pi)
1081224
def format_coord(self, x, y):
        """"""
        p, b = stereonet_math.geographic2plunge_bearing(x, y)
        s, d = stereonet_math.geographic2pole(x, y)
        pb = u'P/B={:0.0f}\u00b0/{:03.0f}\u00b0'.format(p[0], b[0])
        sd = u'S/D={:03.0f}\u00b0/{:0.0f}\u00b0'.format(s[0], d[0])
        return u'{}, {}'.format(pb, sd)
316793
def convert(lat, lon):
    """"""

    x = numpy.cos(lon) * numpy.cos(lat)
    y = numpy.sin(lon) * numpy.cos(lat)
    z = numpy.sin(lat)

    # Invariable plane: values in arcseconds.
    epsilon = 5713.86
    omega = 387390.8

    coseps = numpy.cos(epsilon * secrad)
    sineps = numpy.sin(epsilon * secrad)
    cosom = numpy.cos(omega * secrad)
    sinom = numpy.sin(omega * secrad)

    xi = x * cosom + y * sinom
    yi = coseps * (-sinom * x + cosom * y) + sineps * z
    zi = - sineps * (-sinom * x + cosom * y) + coseps * z

    lat = numpy.arcsin(zi)
    lon = numpy.arctan2(yi, xi)
    return (lat, lon)
14167
def _convert_latitude(self, latitude):
        """"""
        return int((180 - (180 / pi * log(tan(
            pi / 4 + latitude * pi / 360)))) * (2 ** self._zoom) * self._size / 360)
1338924
def pixscale(self):
        
        try:
            (x, y) = self['NAXIS1'] / 2.0, self['NAXIS2'] / 2.0
            p1 = SkyCoord(*self.wcs.xy2sky(x, y) * units.degree)
            p2 = SkyCoord(*self.wcs.xy2sky(x + 1, y + 1) * units.degree)
            return round(p1.separation(p2).to(units.arcsecond).value / math.sqrt(2), 3)
        except Exception as ex:
            logging.debug("Failed to compute PIXSCALE using WCS: {}".format(ex))

        return float(self['PIXSCAL'])
1335915
def ang2pix(nside, theta, phi):
    
    """
    nside, theta, phi = numpy.lib.stride_tricks.broadcast_arrays(nside, theta, phi)

    def equatorial(nside, tt, z):
        t1 = nside * (0.5 + tt)
        t2 = nside * z * 0.75
        jp = (t1 - t2).astype('i8')
        jm = (t1 + t2).astype('i8')
        ir = nside + 1 + jp - jm # in {1, 2n + 1}
        kshift = 1 - (ir & 1) # kshift=1 if ir even, 0 odd

        ip = (jp + jm - nside + kshift + 1) // 2 # in {0, 4n - 1}

        ip = ip % (4 * nside)
        return nside * (nside - 1) * 2 + (ir - 1) * 4 * nside + ip

    def polecaps(nside, tt, z, s):
        tp = tt - numpy.floor(tt)
        za = numpy.abs(z)
        tmp = nside * s / ((1 + za) / 3) ** 0.5
        mp = za > 0.99
        tmp[mp] = nside[mp] * (3 *(1-za[mp])) ** 0.5
        jp = (tp * tmp).astype('i8')
        jm = ((1 - tp) * tmp).astype('i8')
        ir = jp + jm + 1
        ip = (tt * ir).astype('i8')
        ip = ip % (4 * ir)

        r1 = 2 * ir * (ir - 1)
        r2 = 2 * ir * (ir + 1)

        r = numpy.empty_like(r1)

        r[z > 0] = r1[z > 0] + ip[z > 0]
        r[z < 0] = 12 * nside[z < 0] * nside[z < 0] - r2[z < 0] + ip[z < 0]
        return r

    z = numpy.cos(theta)
    s = numpy.sin(theta)

    tt = (phi / (0.5 * numpy.pi) ) % 4 # in [0, 4]

    result = numpy.zeros(z.shape, dtype='i8')
    mask = (z < 2. / 3) & (z > -2. / 3)

    result[mask] = equatorial(nside[mask], tt[mask], z[mask])
    result[~mask] = polecaps(nside[~mask], tt[~mask], z[~mask], s[~mask])
    return result
517792
def gnomonicImageToSphere(x, y):
    
    # Convert angle to [-180, 180] interval
    x = x - 360.*(x>180)
    x = np.asarray(x)
    y = np.asarray(y)
    lon = np.degrees(np.arctan2(y, x))
    r_theta = np.sqrt(x**2 + y**2)
    lat = np.degrees(np.arctan(180. / (np.pi * r_theta)))
    return lon, lat
553881
def projection_to_raster_coords(self, lat, lon):
        
        r_px_py = np.array([1, lon, lat])
        tg = inv(np.array([[1, 0, 0], self.geotransform[0:3], self.geotransform[3:6]]))
        return np.inner(tg, r_px_py)[1:]
1293560
def make_dom_map(pmt_directions, values, nside=512, d=0.2, smoothing=0.1):
    
    import healpy as hp
    discs = [hp.query_disc(nside, dir, 0.2) for dir in pmt_directions]
    npix = hp.nside2npix(nside)
    pixels = np.zeros(npix)
    for disc, value in zip(discs, values):
        for d in disc:
            pixels[d] = value
    if smoothing > 0:
        return hp.sphtfunc.smoothing(pixels, fwhm=smoothing, iter=1)
    return pixels
35836
def eqCoords(lon, lat):
    
    # Convert to radians
    _lambda = math.radians(lon)
    _beta = math.radians(lat)
    _epson = math.radians(23.44)  # The earth's inclination
    
    # Declination in radians
    decl = math.asin(math.sin(_epson) * math.sin(_lambda) * math.cos(_beta) + \
           math.cos(_epson) * math.sin(_beta))
    
    # Equatorial Distance in radians
    ED = math.acos(math.cos(_lambda) * math.cos(_beta) / math.cos(decl))
    
    # RA in radians
    ra = ED if lon < 180 else math.radians(360) - ED
    
    # Correctness of RA if longitude is close to 0º or 180º in a radius of 5º
    if (abs(angle.closestdistance(lon, 0)) < 5 or 
        abs(angle.closestdistance(lon, 180)) < 5):
            a = math.sin(ra) * math.cos(decl)
            b = math.cos(_epson) * math.sin(_lambda) * math.cos(_beta) - \
                math.sin(_epson) * math.sin(_beta)
            if (math.fabs(a-b) > 0.0003):
                ra = math.radians(360) - ra
                
    return (math.degrees(ra), math.degrees(decl))
383384
function (radians) {
            var lon = radians % this.TWO_PI;

            return lon > Math.PI ? lon - this.TWO_PI : lon < -Math.PI ? this.TWO_PI + lon : lon;
        }