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 -
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 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_project wrong? Or did I miss some more basic maths about projections of points onto spheres?