lsst.geom.geometry Namespace Reference

Classes
class	_SubList
class	PartitionMap
class	SphericalBox
class	SphericalBoxPartitionMap
class	SphericalCircle
class	SphericalConvexPolygon
class	SphericalEllipse
class	SphericalRegion

Functions
def	isinf (x)
def	dot (v1, v2)
def	cross (v1, v2)
def	invScale (v, s)
def	normalize (v)
def	sphericalCoords (args)
def	cartesianUnitVector (args)
def	sphericalAngularSep (p1, p2)
def	clampPhi (phi)
def	reduceTheta (theta)
def	northEast (v)
def	alpha (r, centerPhi, phi)
def	maxAlpha (r, centerPhi)
def	cartesianAngularSep (v1, v2)
def	minEdgeSep (p, n, v1, v2)
def	minPhiEdgeSep (p, phi, minTheta, maxTheta)
def	minThetaEdgeSep (p, theta, minPhi, maxPhi)
def	centroid (vertices)
def	between (p, n, v1, v2)
def	segments (phiMin, phiMax, width)
def	median (array)
def	hemispherical (points)
def	convex (vertices)
def	convexHull (points)

Variables
float	ARCSEC_PER_DEG = 3600.0
float	DEG_PER_ARCSEC = 1.0 / 3600.0
float	ANGLE_EPSILON = 0.001 * DEG_PER_ARCSEC
float	POLE_EPSILON = 1.0 * DEG_PER_ARCSEC
float	COS_MAX = 1.0 - 1.0e-15
int	CROSS_N2MIN = 2e-15
	SIN_MIN = math.sqrt(CROSS_N2MIN)
	INF = float('inf')
	NEG_INF = float('-inf')
	MAX_FLOAT = float_info.max
	MIN_FLOAT = float_info.min
	EPSILON = float_info.epsilon
	isinf = math.isinf

Function Documentation

alpha()

def lsst.geom.geometry.alpha	(		r,
			centerPhi,
			phi
	)

Returns the longitude angle alpha of the intersections (alpha, phi),
(-alpha, phi) of the circle centered on (0, centerPhi) with radius r and
the plane z = sin(phi). If there is no intersection, None is returned.

Examples:

forEachPixel.cc.

Definition at line 222 of file geometry.py.

def alpha(r, centerPhi, phi):
    """Returns the longitude angle alpha of the intersections (alpha, phi),
    (-alpha, phi) of the circle centered on (0, centerPhi) with radius r and
    the plane z = sin(phi). If there is no intersection, None is returned.
    """
    if phi < centerPhi - r or phi > centerPhi + r:
        return None
    a = abs(centerPhi - phi)
    if a <= r - (90.0 - phi) or a <= r - (90.0 + phi):
        return None
    p = math.radians(phi)
    cp = math.radians(centerPhi)
    x = math.cos(math.radians(r)) - math.sin(cp) * math.sin(cp)
    u = math.cos(cp) * math.cos(p)
    y = math.sqrt(abs(u * u - x * x))
    return math.degrees(abs(math.atan2(y, x)))

between()

def lsst.geom.geometry.between	(		p,
			n,
			v1,
			v2
	)

Tests whether p lies on the shortest great circle arc from cartesian
unit vectors v1 to v2, assuming that p is a unit vector on the plane
defined by the origin, v1 and v2.

Definition at line 329 of file geometry.py.

def between(p, n, v1, v2):
    """Tests whether p lies on the shortest great circle arc from cartesian
    unit vectors v1 to v2, assuming that p is a unit vector on the plane
    defined by the origin, v1 and v2.
    """
    p1 = cross(n, v1)
    p2 = cross(n, v2)
    return dot(p1, p) >= 0.0 and dot(p2, p) <= 0.0

cartesianAngularSep()

def lsst.geom.geometry.cartesianAngularSep	(		v1,
			v2
	)

Returns the angular separation in degrees between points v1 and v2,
which must both be specified as cartesian 3-vectors.

Definition at line 260 of file geometry.py.

def cartesianAngularSep(v1, v2):
    """Returns the angular separation in degrees between points v1 and v2,
    which must both be specified as cartesian 3-vectors.
    """
    cs = dot(v1, v2)
    n = cross(v1, v2)
    ss = math.sqrt(dot(n, n))
    if cs == 0.0 and ss == 0.0:
        return 0.0
    return math.degrees(math.atan2(ss, cs))

cartesianUnitVector()

def lsst.geom.geometry.cartesianUnitVector

(

args

)

Returns a unit cartesian 3-vector corresponding to the input
coordinates, which can be spherical or 3D cartesian. The 2 (spherical)
or 3 (cartesian 3-vector) input coordinates can either be passed
individually or as a tuple/list, and can be of any type convertible
to a float.

Definition at line 141 of file geometry.py.

def cartesianUnitVector(*args):
    """Returns a unit cartesian 3-vector corresponding to the input
    coordinates, which can be spherical or 3D cartesian. The 2 (spherical)
    or 3 (cartesian 3-vector) input coordinates can either be passed
    individually or as a tuple/list, and can be of any type convertible
    to a float.
    """
    if len(args) == 1:
        args = args[0]
    if len(args) == 2:
        theta = math.radians(float(args[0]))
        phi = math.radians(float(args[1]))
        cosPhi = math.cos(phi)
        return (math.cos(theta) * cosPhi,
                math.sin(theta) * cosPhi,
                math.sin(phi))
    elif len(args) == 3:
        x = float(args[0])
        y = float(args[1])
        z = float(args[2])
        n = math.sqrt(x * x + y * y + z * z)
        if n == 0.0:
            raise RuntimeError('Cannot normalize a 3-vector with 0 magnitude')
        return (x / n, y / n, z / n)
    raise TypeError('Expecting 2 or 3 coordinates, or a tuple/list ' +
                    'containing 2 or 3 coordinates')

centroid()

def lsst.geom.geometry.centroid

(

vertices

)

Computes the centroid of a set of vertices (which must be passed in
as a list of cartesian unit vectors) on the unit sphere.

Definition at line 316 of file geometry.py.

def centroid(vertices):
    """Computes the centroid of a set of vertices (which must be passed in
    as a list of cartesian unit vectors) on the unit sphere.
    """
    x, y, z = 0.0, 0.0, 0.0
    # Note: could use more accurate summation routines
    for v in vertices:
        x += v[0]
        y += v[1]
        z += v[2]
    return normalize((x, y, z))

clampPhi()

def lsst.geom.geometry.clampPhi

(

phi

)

Clamps the input latitude angle phi to [-90.0, 90.0] deg.

Definition at line 184 of file geometry.py.

def clampPhi(phi):
    """Clamps the input latitude angle phi to [-90.0, 90.0] deg.
    """
    if phi < -90.0:
        return -90.0
    elif phi >= 90.0:
        return 90.0
    return phi

convex()

def lsst.geom.geometry.convex

(

vertices

)

Tests whether an ordered list of vertices (which must be specified
as cartesian unit vectors) form a spherical convex polygon and determines
their winding order. Returns a 2-tuple as follows:

(True, True):   The vertex list forms a spherical convex polygon and is in
                counter-clockwise order when viewed from outside the unit
                sphere in a right handed coordinate system.
(True, False):  The vertex list forms a spherical convex polygon and is in
                in clockwise order.
(False, msg):   The vertex list does not form a spherical convex polygon -
                msg is a string describing why the test failed.

The algorithm completes in O(n) time, where n is the number of
input vertices.

Definition at line 1782 of file geometry.py.

def convex(vertices):
    """Tests whether an ordered list of vertices (which must be specified
    as cartesian unit vectors) form a spherical convex polygon and determines
    their winding order. Returns a 2-tuple as follows:

    (True, True):   The vertex list forms a spherical convex polygon and is in
                    counter-clockwise order when viewed from outside the unit
                    sphere in a right handed coordinate system.
    (True, False):  The vertex list forms a spherical convex polygon and is in
                    in clockwise order.
    (False, msg):   The vertex list does not form a spherical convex polygon -
                    msg is a string describing why the test failed.

    The algorithm completes in O(n) time, where n is the number of
    input vertices.
    """
    if len(vertices) < 3:
        return (False, '3 or more vertices must be specified')
    if not hemispherical(vertices):
        return (False, 'vertices are not hemispherical')
    center = centroid(vertices)
    windingAngle = 0.0
    clockwise = False
    counterClockwise = False
    p1 = cross(center, vertices[-1])
    n2 = dot(p1, p1)
    if abs(n2) < CROSS_N2MIN:
        return (False, 'centroid of vertices is too close to a vertex')
    for i in range(len(vertices)):
        beg = vertices[i - 2]
        mid = vertices[i - 1]
        end = vertices[i]
        plane = cross(mid, end)
        n2 = dot(plane, plane)
        if dot(mid, end) >= COS_MAX or n2 < CROSS_N2MIN:
            return (False, 'vertex list contains [near-]duplicates')
        plane = invScale(plane, math.sqrt(n2))
        d = dot(plane, beg)
        if d > SIN_MIN:
            if clockwise:
                return (False, 'vertices wind around centroid in both ' +
                        'clockwise and counter-clockwise order')
            counterClockwise = True
        elif d < -SIN_MIN:
            if counterClockwise:
                return (False, 'vertices wind around centroid in both ' +
                        'clockwise and counter-clockwise order')
            clockwise = True
        # center must be inside polygon formed by vertices if they form a
        # convex polygon - an equivalent check is to test that polygon
        # vertices always wind around center in the same direction.
        d = dot(plane, center)
        if d < SIN_MIN and counterClockwise or d > -SIN_MIN and clockwise:
            return (False, 'centroid of vertices is not conclusively ' +
                    'inside all edges')
        # sum up winding angle for edge (mid, end)
        p2 = cross(center, end)
        n2 = dot(p2, p2)
        if abs(n2) < CROSS_N2MIN:
            return (False, 'centroid of vertices is too close to a vertex')
        windingAngle += cartesianAngularSep(p1, p2)
        p1 = p2
    # For convex polygons, the closest multiple of 360 to
    # windingAngle is 1.
    m = math.floor(windingAngle / 360.0)
    if m == 0.0 and windingAngle > 180.0 or m == 1.0 and windingAngle < 540.0:
        return (True, counterClockwise)
    return (False, 'vertices do not completely wind around centroid, or ' +
            'wind around it multiple times')

convexHull()

def lsst.geom.geometry.convexHull

(

points

)

Computes the convex hull (a spherical convex polygon) of an unordered
list of points on the unit sphere, which must be passed in as cartesian
unit vectors. The algorithm takes O(n log n) time, where n is the number
of points.

Definition at line 1853 of file geometry.py.

def convexHull(points):
    """Computes the convex hull (a spherical convex polygon) of an unordered
    list of points on the unit sphere, which must be passed in as cartesian
    unit vectors. The algorithm takes O(n log n) time, where n is the number
    of points.
    """
    if len(points) < 3:
        return None
    if not hemispherical(points):
        return None
    center = centroid(points)
    maxSep = 0.0
    extremum = 0
    # Vertex furthest from the center is on the hull
    for i in range(len(points)):
        sep = cartesianAngularSep(points[i], center)
        if sep > maxSep:
            extremum = i
            maxSep = sep
    anchor = points[extremum]
    refPlane = cross(center, anchor)
    n2 = dot(refPlane, refPlane)
    if n2 < CROSS_N2MIN:
        # vertex and center are too close
        return None
    refPlane = invScale(refPlane, math.sqrt(n2))
    # Order points by winding angle from the first (extreme) vertex
    av = [(0.0, anchor)]
    for i in range(0, len(points)):
        if i == extremum:
            continue
        v = points[i]
        plane = cross(center, v)
        n2 = dot(plane, plane)
        if n2 >= CROSS_N2MIN:
            plane = invScale(plane, math.sqrt(n2))
            p = cross(refPlane, plane)
            sa = math.sqrt(dot(p, p))
            if dot(p, center) < 0.0:
                sa = -sa
            angle = math.atan2(sa, dot(refPlane, plane))
            if angle < 0.0:
                angle += 2.0 * math.pi
            av.append((angle, v))
    if len(av) < 3:
        return None
    av.sort(key=lambda t: t[0])  # stable, so av[0] still contains anchor
    # Loop over vertices using a Graham scan adapted for spherical geometry.
    # Chan's algorithm would be an improvement, but seems likely to be slower
    # for small vertex counts (the expected case).
    verts = [anchor]
    edges = [None]
    edge = None
    i = 1
    while i < len(av):
        v = av[i][1]
        p = cross(anchor, v)
        n2 = dot(p, p)
        if dot(anchor, v) < COS_MAX and n2 >= CROSS_N2MIN:
            if edge is None:
                # Compute first edge
                edge = invScale(p, math.sqrt(n2))
                verts.append(v)
                edges.append(edge)
                anchor = v
            else:
                d = dot(v, edge)
                if d > SIN_MIN:
                    # v is inside the edge defined by the last
                    # 2 vertices on the hull
                    edge = invScale(p, math.sqrt(n2))
                    verts.append(v)
                    edges.append(edge)
                    anchor = v
                elif d < -SIN_MIN:
                    # backtrack - the most recently added hull vertex
                    # is not actually on the hull.
                    verts.pop()
                    edges.pop()
                    anchor = verts[-1]
                    edge = edges[-1]
                    # reprocess v to decide whether to add it to the hull
                    # or whether further backtracking is necessary
                    continue
                # v is coplanar with edge, skip it
        i += 1
    # Handle backtracking necessary for last edge
    if len(verts) < 3:
        return None
    v = verts[0]
    while True:
        p = cross(anchor, v)
        n2 = dot(p, p)
        if dot(anchor, v) < COS_MAX and n2 >= CROSS_N2MIN:
            if dot(v, edge) > SIN_MIN:
                edges[0] = invScale(p, math.sqrt(n2))
                break
        verts.pop()
        edges.pop()
        anchor = verts[-1]
        edge = edges[-1]
        if len(verts) < 3:
            return None
    return SphericalConvexPolygon(verts, edges)


# -- Generic unit sphere partitioning scheme base class ----

# TODO: make this an ABC once python 2.6 becomes the default

cross()

def lsst.geom.geometry.cross	(		v1,
			v2
	)

Returns the cross product of cartesian 3-vectors v1 and v2.

Definition at line 85 of file geometry.py.

def cross(v1, v2):
    """Returns the cross product of cartesian 3-vectors v1 and v2.
    """
    return (v1[1] * v2[2] - v1[2] * v2[1],
            v1[2] * v2[0] - v1[0] * v2[2],
            v1[0] * v2[1] - v1[1] * v2[0])

dot()

def lsst.geom.geometry.dot	(		v1,
			v2
	)

Returns the dot product of cartesian 3-vectors v1 and v2.

Definition at line 79 of file geometry.py.

def dot(v1, v2):
    """Returns the dot product of cartesian 3-vectors v1 and v2.
    """
    return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]

hemispherical()

def lsst.geom.geometry.hemispherical

(

points

)

Tests whether a set of points is hemispherical, i.e. whether a plane
exists such that all points are strictly on one side of the plane. The
algorithm used is Megiddo's algorithm for linear programming in R2 and
has run-time O(n), where n is the number of points. Points must be passed
in as a list of cartesian unit vectors.

Definition at line 1740 of file geometry.py.

def hemispherical(points):
    """Tests whether a set of points is hemispherical, i.e. whether a plane
    exists such that all points are strictly on one side of the plane. The
    algorithm used is Megiddo's algorithm for linear programming in R2 and
    has run-time O(n), where n is the number of points. Points must be passed
    in as a list of cartesian unit vectors.
    """
    # Check whether the set of linear equations
    # x * v[0] + y * v[1] + z * v[2] > 0.0 (for v in points)
    # has a solution (x, y, z). If (x,y,z) is a solution (is feasible),
    # so is C*(x,y,z), C > 0. Therefore we can fix z to 1, -1 and
    # perform 2D feasibility tests.
    if _feasible2D(points, 1.0):
        return True
    if _feasible2D(points, -1.0):
        return True
    # At this point a feasible solution must have z = 0. Fix y to 1, -1 and
    # perform 1D feasibility tests.
    if _feasible1D(points, 1.0):
        return True
    if _feasible1D(points, -1.0):
        return True
    # At this point a feasible solution must have y = z = 0. If all v[0]
    # are non-zero and have the same sign, then there is a feasible solution.
    havePos = False
    haveNeg = False
    for v in points:
        if v[0] > 0.0:
            if haveNeg:
                return False
            havePos = True
        elif v[0] < 0.0:
            if havePos:
                return False
            haveNeg = True
        else:
            return False
    return True


# -- Convexity test and convex hull algorithm ----

invScale()

def lsst.geom.geometry.invScale	(		v,
			s
	)

Returns a copy of the cartesian 3-vector v scaled by 1 / s.

Definition at line 93 of file geometry.py.

def invScale(v, s):
    """Returns a copy of the cartesian 3-vector v scaled by 1 / s.
    """
    return (v[0] / s, v[1] / s, v[2] / s)

isinf()

def lsst.geom.geometry.isinf

(

x

)

Definition at line 73 of file geometry.py.

    def isinf(x):
        return x == INF or x == NEG_INF

# -- Utility methods ----

maxAlpha()

def lsst.geom.geometry.maxAlpha	(		r,
			centerPhi
	)

Computes alpha, the extent in longitude angle [-alpha, alpha] of the
circle with radius r and center (0, centerPhi) on the unit sphere.
Both r and centerPhi are assumed to be in units of degrees.
centerPhi is clamped to lie in the range [-90,90] and r must
lie in the range [0, 90].

Definition at line 240 of file geometry.py.

def maxAlpha(r, centerPhi):
    """Computes alpha, the extent in longitude angle [-alpha, alpha] of the
    circle with radius r and center (0, centerPhi) on the unit sphere.
    Both r and centerPhi are assumed to be in units of degrees.
    centerPhi is clamped to lie in the range [-90,90] and r must
    lie in the range [0, 90].
    """
    assert r >= 0.0 and r <= 90.0
    if r == 0.0:
        return 0.0
    centerPhi = clampPhi(centerPhi)
    if abs(centerPhi) + r > 90.0 - POLE_EPSILON:
        return 180.0
    r = math.radians(r)
    c = math.radians(centerPhi)
    y = math.sin(r)
    x = math.sqrt(abs(math.cos(c - r) * math.cos(c + r)))
    return math.degrees(abs(math.atan(y / x)))

median()

def lsst.geom.geometry.median

(

array

)

Finds the median element of the given array in linear time.

Examples:

imageStatistics.cc.

Definition at line 1544 of file geometry.py.

def median(array):
    """Finds the median element of the given array in linear time.
    """
    left = 0
    right = len(array)
    if right == 0:
        return None
    k = right >> 1
    while True:
        i = _medianOfMedians(array, left, right)
        i = _partition(array, left, right, i)
        if k == i:
            return array[k]
        elif k < i:
            right = i
        else:
            left = i + 1


# -- Testing whether a set of points is hemispherical ----
#
# This test is used by the convexity test and convex hull algorithm. It is
# implemented using Megiddo's algorithm for linear programming in R2, see:
#
# Megiddo, N. 1982. Linear-time algorithms for linear programming in R3 and related problems.
# In Proceedings of the 23rd Annual Symposium on Foundations of Computer Science (November 03 - 05, 1982).
# SFCS. IEEE Computer Society, Washington, DC, 329-338. DOI=
# http://dx.doi.org/10.1109/SFCS.1982.74

minEdgeSep()

def lsst.geom.geometry.minEdgeSep	(		p,
			n,
			v1,
			v2
	)

Returns the minimum angular separation in degrees between p
and points on the great circle edge with plane normal n and
vertices v1, v2. All inputs must be unit cartesian 3-vectors.

Definition at line 272 of file geometry.py.

def minEdgeSep(p, n, v1, v2):
    """Returns the minimum angular separation in degrees between p
    and points on the great circle edge with plane normal n and
    vertices v1, v2. All inputs must be unit cartesian 3-vectors.
    """
    p1 = cross(n, v1)
    p2 = cross(n, v2)
    if dot(p1, p) >= 0.0 and dot(p2, p) <= 0.0:
        return abs(90.0 - cartesianAngularSep(p, n))
    else:
        return min(cartesianAngularSep(p, v1), cartesianAngularSep(p, v2))

minPhiEdgeSep()

def lsst.geom.geometry.minPhiEdgeSep	(		p,
			phi,
			minTheta,
			maxTheta
	)

Returns the minimum angular separation in degrees between p
and points on the small circle edge with constant latitude angle
phi and vertices (minTheta, phi), (maxTheta, phi). p must be in
spherical coordinates.

Definition at line 285 of file geometry.py.

def minPhiEdgeSep(p, phi, minTheta, maxTheta):
    """Returns the minimum angular separation in degrees between p
    and points on the small circle edge with constant latitude angle
    phi and vertices (minTheta, phi), (maxTheta, phi). p must be in
    spherical coordinates.
    """
    if minTheta > maxTheta:
        if p[0] >= minTheta or p[0] <= maxTheta:
            return abs(p[1] - phi)
    else:
        if p[0] >= minTheta and p[0] <= maxTheta:
            return abs(p[1] - phi)
    return min(sphericalAngularSep(p, (minTheta, phi)),
               sphericalAngularSep(p, (maxTheta, phi)))

minThetaEdgeSep()

def lsst.geom.geometry.minThetaEdgeSep	(		p,
			theta,
			minPhi,
			maxPhi
	)

Returns the minimum angular separation in degrees between p
and points on the great circle edge with constant longitude angle
theta and vertices (theta, minPhi), (theta, maxPhi). p must be a
unit cartesian 3-vector.

Definition at line 301 of file geometry.py.

def minThetaEdgeSep(p, theta, minPhi, maxPhi):
    """Returns the minimum angular separation in degrees between p
    and points on the great circle edge with constant longitude angle
    theta and vertices (theta, minPhi), (theta, maxPhi). p must be a
    unit cartesian 3-vector.
    """
    v1 = cartesianUnitVector(theta, minPhi)
    v2 = cartesianUnitVector(theta, maxPhi)
    n = cross(v1, v2)
    d2 = dot(n, n)
    if d2 == 0.0:
        return min(cartesianAngularSep(p, v1), cartesianAngularSep(p, v2))
    return minEdgeSep(p, invScale(n, math.sqrt(d2)), v1, v2)

normalize()

def lsst.geom.geometry.normalize

(

v

)

Returns a normalized copy of the cartesian 3-vector v.

Definition at line 99 of file geometry.py.

def normalize(v):
    """Returns a normalized copy of the cartesian 3-vector v.
    """
    n = math.sqrt(dot(v, v))
    if n == 0.0:
        raise RuntimeError('Cannot normalize a 3-vector with 0 magnitude')
    return (v[0] / n, v[1] / n, v[2] / n)

northEast()

def lsst.geom.geometry.northEast

(

v

)

Returns unit N,E basis vectors for a point v, which must be a
cartesian 3-vector.

Definition at line 205 of file geometry.py.

def northEast(v):
    """Returns unit N,E basis vectors for a point v, which must be a
    cartesian 3-vector.
    """
    north = (-v[0] * v[2],
             -v[1] * v[2],
             v[0] * v[0] + v[1] * v[1])
    if north == (0.0, 0.0, 0.0):
        # pick an arbitrary orthogonal basis with z = 0
        north = (-1.0, 0.0, 0.0)
        east = (0.0, 1.0, 0.0)
    else:
        north = normalize(north)
        east = normalize(cross(north, v))
    return north, east

reduceTheta()

def lsst.geom.geometry.reduceTheta

(

theta

)

Range reduces the given longitude angle to lie in the range
[0.0, 360.0).

Definition at line 194 of file geometry.py.

def reduceTheta(theta):
    """Range reduces the given longitude angle to lie in the range
    [0.0, 360.0).
    """
    theta = math.fmod(theta, 360.0)
    if theta < 0.0:
        return theta + 360.0
    else:
        return theta

segments()

def lsst.geom.geometry.segments	(		phiMin,
			phiMax,
			width
	)

Computes the number of segments to divide the given latitude angle
range [phiMin, phiMax] (degrees) into. Two points within the range
separated by at least one segment are guaranteed to have angular
separation of at least width degrees.

Definition at line 339 of file geometry.py.

def segments(phiMin, phiMax, width):
    """Computes the number of segments to divide the given latitude angle
    range [phiMin, phiMax] (degrees) into. Two points within the range
    separated by at least one segment are guaranteed to have angular
    separation of at least width degrees.
    """
    p = max(abs(phiMin), abs(phiMax))
    if p > 90.0 - 1.0 * DEG_PER_ARCSEC:
        return 1
    if width >= 180.0:
        return 1
    elif width < 1.0 * DEG_PER_ARCSEC:
        width = 1.0 * DEG_PER_ARCSEC
    p = math.radians(p)
    cw = math.cos(math.radians(width))
    sp = math.sin(p)
    cp = math.cos(p)
    x = cw - sp * sp
    u = cp * cp
    y = math.sqrt(abs(u * u - x * x))
    return int(math.floor((2.0 * math.pi) / abs(math.atan2(y, x))))


# -- Regions on the unit sphere ----

# TODO: make this an ABC once python 2.6 becomes the default

sphericalAngularSep()

def lsst.geom.geometry.sphericalAngularSep	(		p1,
			p2
	)

Returns the angular separation in degrees between points p1 and p2,
which must both be specified in spherical coordinates. The implementation
uses the halversine distance formula.

Definition at line 169 of file geometry.py.

def sphericalAngularSep(p1, p2):
    """Returns the angular separation in degrees between points p1 and p2,
    which must both be specified in spherical coordinates. The implementation
    uses the halversine distance formula.
    """
    sdt = math.sin(math.radians(p1[0] - p2[0]) * 0.5)
    sdp = math.sin(math.radians(p1[1] - p2[1]) * 0.5)
    cc = math.cos(math.radians(p1[1])) * math.cos(math.radians(p2[1]))
    s = math.sqrt(sdp * sdp + cc * sdt * sdt)
    if s > 1.0:
        return 180.0
    else:
        return 2.0 * math.degrees(math.asin(s))

sphericalCoords()

def lsst.geom.geometry.sphericalCoords

(

args

)

Returns spherical coordinates in degrees for the input coordinates,
which can be spherical or 3D cartesian. The 2 (spherical) or 3
(cartesian 3-vector) inputs can be passed either individually
or as a tuple/list, and can be of any type convertible to a float.

Definition at line 108 of file geometry.py.

def sphericalCoords(*args):
    """Returns spherical coordinates in degrees for the input coordinates,
    which can be spherical or 3D cartesian. The 2 (spherical) or 3
    (cartesian 3-vector) inputs can be passed either individually
    or as a tuple/list, and can be of any type convertible to a float.
    """
    if len(args) == 1:
        args = args[0]
    if len(args) == 2:
        t = (float(args[0]), float(args[1]))
        if t[1] < -90.0 or t[1] > 90.0:
            raise RuntimeError('Latitude angle is out of bounds')
        return t
    elif len(args) == 3:
        x = float(args[0])
        y = float(args[1])
        z = float(args[2])
        d2 = x * x + y * y
        if d2 == 0.0:
            theta = 0.0
        else:
            theta = math.degrees(math.atan2(y, x))
            if theta < 0.0:
                theta += 360.0
        if z == 0.0:
            phi = 0.0
        else:
            phi = math.degrees(math.atan2(z, math.sqrt(d2)))
        return (theta, phi)
    raise TypeError('Expecting 2 or 3 coordinates, or a tuple/list ' +
                    'containing 2 or 3 coordinates')

Variable Documentation

ANGLE_EPSILON

float lsst.geom.geometry.ANGLE_EPSILON = 0.001 * DEG_PER_ARCSEC

Definition at line 37 of file geometry.py.

ARCSEC_PER_DEG

float lsst.geom.geometry.ARCSEC_PER_DEG = 3600.0

Definition at line 33 of file geometry.py.

COS_MAX

float lsst.geom.geometry.COS_MAX = 1.0 - 1.0e-15

Definition at line 44 of file geometry.py.

CROSS_N2MIN

int lsst.geom.geometry.CROSS_N2MIN = 2e-15

Definition at line 48 of file geometry.py.

DEG_PER_ARCSEC

float lsst.geom.geometry.DEG_PER_ARCSEC = 1.0 / 3600.0

Definition at line 34 of file geometry.py.

EPSILON

float lsst.geom.geometry.EPSILON = float_info.epsilon

Definition at line 62 of file geometry.py.

INF

lsst.geom.geometry.INF = float('inf')

Definition at line 55 of file geometry.py.

isinf

lsst.geom.geometry.isinf = math.isinf

Definition at line 71 of file geometry.py.

MAX_FLOAT

float lsst.geom.geometry.MAX_FLOAT = float_info.max

Definition at line 60 of file geometry.py.

MIN_FLOAT

float lsst.geom.geometry.MIN_FLOAT = float_info.min

Definition at line 61 of file geometry.py.

NEG_INF

lsst.geom.geometry.NEG_INF = float('-inf')

Definition at line 56 of file geometry.py.

POLE_EPSILON

float lsst.geom.geometry.POLE_EPSILON = 1.0 * DEG_PER_ARCSEC

Definition at line 40 of file geometry.py.

SIN_MIN

lsst.geom.geometry.SIN_MIN = math.sqrt(CROSS_N2MIN)

Definition at line 52 of file geometry.py.