Why are my icosahedron triangle subdivisions not equilateral (PostGIS)?

Question

I am creating a geodesic polyhedron using PostGIS. I am dividing the base icosahedron points (see Appendix A), into the first level of Class I subdivisions.

The setup

For simplicity let's take just the first triangle formed by these points:

CREATE TABLE pix (
    id bigserial
  , name text
  , geog geography(POINT,4326) -- use WGS84
);

insert into pix values (DEFAULT, 'China', st_point(122.3, 39.1));
insert into pix values (DEFAULT, 'Norway', st_point(10.53619898, 64.7));
insert into pix values (DEFAULT, 'Arabian sea', st_point(58.15770555, 10.44734504));

Which creates the base icosahedron triangle

The numbers of the triangle sides correspond to the row numbers shown in the table generated by this query

select p1.name
  , p2.name
  , st_distance(p1.geog, p2.geog) as dist
from pix p1
  cross join pix p2
where p1.id < p2.id

The points are roughly equidistant. There are small variances because the seed data above seems to have been generated assuming a sphere, not spheroid, which st_distance is returning - use_spheroid false in st_distance will yield much closer distances, but note that using false in all the presented queries does not change the problem described below.

First level of Class I subdivision

insert into pix (name, geog)
select p1.name || '-' || p2.name
  , st_project(p1.geog, st_distance(p1.geog, p2.geog) / 2, st_azimuth(p1.geog, p2.geog)) as geog
from pix p1
  cross join pix p2
where p1.id < p2.id

The above generates the mid-points of each of the pairs of points on the original triangle

And these have the following distances, generated by this query

select p1.name
  , p2.name
  , st_distance(p1.geog, p2.geog) as dist
from pix p1
  cross join pix p2
where p1.id < p2.id
  and p2.id > 3
order by dist

The problem

The distances of the segments 7, 8 and 9 (around 4000kms) are much longer than the distances of segments 1-6 (around 3500kms).

Why is this algorithm not creating something more close to equliateral triangles? Is my use of st_azimuth and st_project wrong? Or did I miss some more basic maths about projections of points onto spheres?

Answer 1

Yes, this is due to the tiled icosahedron to-sphere projections. I think of this process as:

Take a normal icosahedron, all corners lie on the same sphere, and therefore, the distance is also the same, when projected on a sphere. However, if we take our non-projected icosahedron and subdivide all triangles, the generated points will not be on that sphere. When these move outwards, they will be bigger, like the redtriangles in this image.

Also, these constructions are used in architecture. Using the Geodesic dome calculator, (image from their site), we find for an edge length B of 3.5, the corresponding length for A is 3.9:

So indeed it was the mathematics that you missed.