I am working on a Python script to:

  1. compute the vertex coordinates of a geodesic sphere/icosahedron,
  2. project the triangles onto a sphere, then
  3. find which spherical triangle contains an arbitrary point.

I would ultimately like to provide my script with a list of sample point coordinates, and the script would indicate which triangles contain points. These characters would then be used for presence-absence data in biological studies.

I'm mostly done, but struggling with step 3. Is there a fast algorithm to identify whether a point is contained within a spherical equilateral triangle? In two dimensions for general polygons, you can count how many times a ray crosses a perimeter, but this seems challenging to do in three dimensions since the perimeter is curved.

It seems like there should be a clever way to exploit the symmetry of spherical equilateral triangles, but I'm just not seeing it. I don't work with GIS or spherical trigonometry, so just let me know if my question is improperly phrased.

  • 2
    Welcome to GIS SE! However, looking at your question, it does not look like there is a clearly GIS component to it so I am wondering whether it might not be better posted to somewhere like our sister site Mathematics. – PolyGeo Nov 12 '13 at 20:39

Using the spherical coordinates, create the icosahedral tessellation and project it with a suitable projection for your work, then apply any GIS to solve this point in polygon problem.

If you want to write your own code, note that each pair of vertices determines a reflection plane (containing the origin) and that points can be classified according to which side of each of the 15 distinct reflection planes they lie on. This classification is actually a finer one than requested, because each of the 20 faces of the icosahedron is further divided into six congruent triangles, giving 120 regions (equal to a fundamental domain for the action of the reflections together with its 119 images).

For the record, a set of normal vectors to these planes (in geocentric Cartesian coordinates) is

X   Y   Z
1-f -f  1
f   -1  f-1
-2  0   0
f-1 -f  1
-f  -1  f-1
-1  f-1 f
1-f f   1
0   -2  0
-1  1-f f
-f  1   f-1
1   1-f f
0   0   2
f-1 f   1
f   1   f-1
1   f-1 f

where f = (1 + sqrt(5))/2 = 1.61803398874989... is the Golden Ratio. (These all point towards the Equator or northwards.) To use them, convert any point (on the sphere or off) into geocentric Cartesian coordinates (x,y,z) and form the 15 dot products with the 15 normals (which essentially multiplies the 15 by 3 matrix by the 3 by 1 vector (x,y,z)'). Classify those products as negative or non-negative, yielding a 15-element Boolean vector. Considering that as a 15-bit binary integer, there are only 120 possible values, one for each of the 120 little triangles. These values encode the desired information about the face (and which part of the face) where the point (x,y,z) is located.

To illustrate, I associated 120 random colors with these codes to display 10,000 random points on the unit sphere:


Both the icosahedron and its dual, the dodecahedron, are evident.

|improve this answer|||||
  • This is very helpful and clear. It also seems if I would like to bin a point into higher frequency tessellations, once I know which major triangle contains the point, I can recursively use your approach to narrow down which minor triangle contains the point, so I don't have to perform very large matrix-vector multiplies. – mlandis Nov 12 '13 at 22:26
  • You might be interested, then, in the related answer at gis.stackexchange.com/questions/53674/… which demonstrates that recursive approach. – whuber Nov 12 '13 at 23:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.