Nearest Neighbor Search in Spherical Coordinates
Due to accuracy complications as well as loss of information I cannot convert them to Cartesian coordinates.


I don't understand this restriction. Converting to Cartesian coordinates uses essentially the same trigonometric functions that you would need to find the distance between two points in spherical coordinates: the only fundamental difference is that working directly in spherical coordinates you can take the cosine of the difference of azimuthal angles rather than taking sine and cosine of the individual azimuthal angles.

In any case, you should be able to bound the error in converting to Cartesian coordinates and take that error bound into account to create a shortlist of candidates using k-d trees; then recalculate the distances for the points in the shortlist using the spherical coordinates directly.