Assume you have a table of locations with Latitude and Longitude for each one. In my case the table is "station" and the primary key is "LocID".

First we create a view to help with the 3D geometry (6378 = Earth's radius in km):

CREATE VIEW gpsGlb AS
    SELECT 
        LocID
        ,6378 * COS(RADIANS(Latitude)) * COS(RADIANS(Longitude)) AS x
        ,6378 * COS(RADIANS(Latitude)) * SIN(RADIANS(Longitude)) AS y
        ,6378 * SIN(RADIANS(Latitude)) AS z
    FROM station;

Here I query for distances to all my locations that are NOT LocID = 405 (rounded miles in my case):

SELECT 
    LocID
    ,ROUND((2 * 6378 * ASIN(d / 2 / 6378)) * 0.621371192) AS dist_mi
FROM
    (SELECT
        SQRT(dx * dx + dy * dy + dz * dz) AS d
        ,LocID
     FROM
        (SELECT
            p1.x - p2.x AS dx
            ,p1.y - p2.y AS dy
            ,p1.z - p2.z AS dz
            ,p2.LocID
        FROM gpsGlb p1
        JOIN gpsGlb p2 ON (p1.LocID = 405 AND p2.LocID != 405)
       ) t1
    ) t2
ORDER BY dist_mi

Here I get the initial bearing to the locations. The "boxed" calculation will come in handy later.

SELECT
    LocID
    ,(360 + DEGREES(ATAN2(y, x))) % 360 AS initBearing_deg
    ,ROUND(((360 + DEGREES(ATAN2(y, x))) % 360) / 22.5) * 22.5 
     AS initBearingBoxed_deg
FROM
    (SELECT
        SIN(RADIANS(s2.Longitude - s1.Longitude)) * COS(RADIANS(s2.Latitude)) 
        AS y
        ,COS(RADIANS(s1.Latitude)) * SIN(RADIANS(s2.Latitude))
            - SIN(RADIANS(s1.Latitude)) * COS(RADIANS(s2.Latitude))
               * COS(RADIANS(s2.Longitude - s1.Longitude)) 
        AS x
        ,s2.LocID
    FROM station s1
    JOIN station s2 ON (s1.LocID = 405 AND s2.LocID != 405)
    ) q1

Here's the combined query plus boxed degrees converted to 'NNE', etc. I've also added a limit for the distance in the qq1 subquery.

SELECT
    qq2.LocID
    ,dist_mi
    ,CASE initBearingBoxed_deg
        WHEN 22.5 THEN 'NNE'   WHEN 45 THEN 'NE'
        WHEN 67.5 THEN 'ENE'   WHEN 90 THEN 'E'
        WHEN 112.5 THEN 'ESE'  WHEN 135 THEN 'SE'
        WHEN 157.5 THEN 'SSE'  WHEN 180 THEN 'S'
        WHEN 202.5 THEN 'SSW'  WHEN 225 THEN 'SW'
        WHEN 247.5 THEN 'WSW'  WHEN 270 THEN 'W'
        WHEN 292.5 THEN 'WNW'  WHEN 315 THEN 'NW'
        WHEN 337.5 THEN 'NNW'  ELSE 'N'
     END AS bearing
FROM (
    SELECT 
        LocID
        ,ROUND((2 * 6378 * ASIN(d / 2 / 6378)) * 0.621371192) AS dist_mi
    FROM
        (SELECT
            SQRT(dx * dx + dy * dy + dz * dz) AS d
            ,LocID
         FROM
            (SELECT
                p1.x - p2.x AS dx
                ,p1.y - p2.y AS dy
                ,p1.z - p2.z AS dz
                ,p2.LocID
            FROM gpsGlb p1
            JOIN gpsGlb p2 ON (p1.LocID = 405 AND p2.LocID != 405)
           ) t1
        ) t2
    ) qq1
JOIN (
    SELECT
        LocID
        ,(360 + DEGREES(ATAN2(y, x))) % 360 AS initBearing_deg
        ,(360 + ROUND((DEGREES(ATAN2(y, x))) / 22.5) * 22.5) % 360 
         AS initBearingBoxed_deg
    FROM
        (SELECT
            SIN(RADIANS(s2.Longitude - s1.Longitude)) * COS(RADIANS(s2.Latitude)) 
             AS y
            ,COS(RADIANS(s1.Latitude)) * SIN(RADIANS(s2.Latitude))
                - SIN(RADIANS(s1.Latitude)) * COS(RADIANS(s2.Latitude))
                   * COS(RADIANS(s2.Longitude - s1.Longitude)) 
             AS x
            ,s2.LocID
        FROM station s1
        JOIN station s2 ON (s1.LocID = 405 AND s2.LocID != 405)
        ) q1
    ) qq2 ON (qq1.LocID = qq2.LocID
              AND qq1.dist_mi <= 60)
ORDER BY dist_mi

// Multiply a point by a transformation matrix.

/*!
 * Multiply a point by a transformation matrix.
 *
 * Applies the given transformation matrix to the given point.  With some
 * transformation matrices, a vector may also be transformed.
 *
 * \param c Result. (just a float[3])
 * \param m Transformation matrix. (just a float[4][4])
 * \param a Input point. (just a float[3])
 */
void
lib3ds_vector_transform(Lib3dsVector c, Lib3dsMatrix m, Lib3dsVector a) {
    c[0] = m[0][0] * a[0] + m[1][0] * a[1] + m[2][0] * a[2] + m[3][0];
    c[1] = m[0][1] * a[0] + m[1][1] * a[1] + m[2][1] * a[2] + m[3][1];
    c[2] = m[0][2] * a[0] + m[1][2] * a[1] + m[2][2] * a[2] + m[3][2];
}

// Multiply two matrices

/*!
* Multiplies a matrix by a second one (m = m * n).
*/
void lib3ds_matrix_mult(Lib3dsMatrix m, Lib3dsMatrix n) {
    Lib3dsMatrix tmp;
    int i, j, k;
    float ab;

    memcpy(tmp, m, sizeof(Lib3dsMatrix));
    for (j = 0; j < 4; j++) {
        for (i = 0; i < 4; i++) {
            ab = 0.0f;
            for (k = 0; k < 4; k++)
                ab += tmp[k][i] * n[j][k];
            m[j][i] = ab;
        }
    }
}

// Apply a rotation about an arbitrary axis to a matrix.

/*!
* Apply a rotation about an arbitrary axis to a matrix.
* m is just a float[4][4], q is just a float[4]
*/
void lib3ds_matrix_rotate(Lib3dsMatrix m, Lib3dsQuat q) {
    float s, xs, ys, zs, wx, wy, wz, xx, xy, xz, yy, yz, zz, l;
    Lib3dsMatrix R;

    l = q[0] * q[0] + q[1] * q[1] + q[2] * q[2] + q[3] * q[3];
    s = 2.0f / l;

    xs = q[0] * s;
    ys = q[1] * s;
    zs = q[2] * s;
    wx = q[3] * xs;
    wy = q[3] * ys;
    wz = q[3] * zs;
    xx = q[0] * xs;
    xy = q[0] * ys;
    xz = q[0] * zs;
    yy = q[1] * ys;
    yz = q[1] * zs;
    zz = q[2] * zs;

    R[0][0] = 1.0f - (yy + zz);
    R[1][0] = xy - wz;
    R[2][0] = xz + wy;
    R[0][1] = xy + wz;
    R[1][1] = 1.0f - (xx + zz);
    R[2][1] = yz - wx;
    R[0][2] = xz - wy;
    R[1][2] = yz + wx;
    R[2][2] = 1.0f - (xx + yy);
    R[3][0] = R[3][1] = R[3][2] = R[0][3] = R[1][3] = R[2][3] = 0.0f;
    R[3][3] = 1.0f;

    lib3ds_matrix_mult(m, R);
}

// Build a quaternion from a rotation about an arbitrary axis

/*!
* Compute a quaternion from axis and angle.
*
* \param c Computed quaternion (this is just a float[4])
* \param axis Rotation axis (this is just a float[3])
* \param angle Angle of rotation, radians.
*/
void lib3ds_quat_axis_angle(Lib3dsQuat c, Lib3dsVector axis, float angle)
{
    double omega, s;
    double l = sqrt(axis[0] * axis[0] + axis[1] * axis[1] + axis[2] * axis[2]);

    omega = -0.5 * angle;
    s = sin(omega) / l;

    c[0] = (float)s * axis[0];
    c[1] = (float)s * axis[1];
    c[2] = (float)s * axis[2];
    c[3] = (float)cos(omega);
}

Given a negative or huge angle, it returns a number between 0 and 359 :b

//+ Jonas Raoni Soares Silva
//@ http://jsfromhell.com

function realAngle(a){
	return (a %= 360) < 0 ? a + 360 : a;
}

usage

alert(realAngle(359));
alert(realAngle(-45));
alert(realAngle(1234));