While I am aware of vaguely similar questions, I seem to be stuck on this.

Buckminster Fuller introduced a spherical mapping of the world onto an icosahedron – it’s known as the Dymaxion Map

A common way of identifying the cartesian coordinates of an icosahedron is by using the coordinates, where ???? is the golden ratio: (1+√5)/2) or cos(π/5.0)

(0,±1,±1????),(±1,±1????,0),(±1????,0,±1)

Expanding this out gives me the location of the 12 vertices of a regular icosahedron with side length 2:

    ???? = math.cos(π / 5.0)
    ico_vertices = [
        (0, -1, -????), (0, -1, +????), (0, +1, -????), (0, +1, +????),
        (-1, -????, 0), (-1, +????, 0), (+1, -????, 0), (+1, +????, 0),
        (-????, 0, -1), (-????, 0, +1), (+????, 0, -1), (+????, 0, +1)
    ]

Needless to say these need to be normalised.

    iv = np.array(ico_vertices)
    iv_n = ((iv[:, None] ** 2).sum(2) ** 0.5).reshape(-1, 1)
    ico = iv / iv_n  #This is the starting set of vertices.

Here is an image of the golden ratio ???? coordinates projected onto Google Maps.

Fuller’s icosahedron, mapped onto a spherical projection of the globe, defined using these 12 vertices: (the xyz values are already normalised to a unit sphere)

{
    "vertices":[
        {"name":"CHINA",                 "ll": [[39.10000000, "N"],[122.30000000,"E"]], "xyz":[-0.41468222, 0.65596241,  0.63067581]},
        {"name":"NORWAY",                "ll": [[64.70000000, "N"],[10.53619898 ,"E"]], "xyz":[ 0.42015243, 0.07814525,  0.90408255]},
        {"name":"ARABIAN SEA",           "ll": [[10.44734504, "N"],[58.15770555 ,"E"]], "xyz":[ 0.51883673, 0.83542038,  0.18133184]},
        {"name":"LIBERIA",               "ll": [[2.30088201 , "N"],[5.24539058  ,"W"]], "xyz":[ 0.99500944, -0.09134780, 0.04014717]},
        {"name":"PUERTO RICO",           "ll": [[23.71792533, "N"],[67.13232659 ,"W"]], "xyz":[ 0.35578140, -0.84358000, 0.40223423]},
        {"name":"ALASKA",                "ll": [[50.10320164, "N"],[143.47849033,"W"]], "xyz":[-0.51545596, -0.38171689, 0.76720099]},
        {"name":"BUENOS AIRES",          "ll": [[39.10000000, "S"],[57.70000000 ,"W"]], "xyz":[ 0.41468222, -0.65596241,-0.63067581]},
        {"name":"ANTARCTICA",            "ll": [[64.70000000, "S"],[169.46380102,"W"]], "xyz":[-0.42015243, -0.07814525,-0.90408255]},
        {"name":"PITCAIRN ISLAND",       "ll": [[10.44734504, "S"],[121.84229445,"W"]], "xyz":[-0.51883673, -0.83542038,-0.18133184]},
        {"name":"GILBERT ISLAND",        "ll": [[2.30088201 , "S"],[174.75460942,"E"]], "xyz":[-0.99500944, 0.09134780, -0.04014717]},
        {"name":"AUSTRALIA",             "ll": [[23.71792533, "S"],[112.86767341,"E"]], "xyz":[-0.35578140, 0.84358000, -0.40223423]},
        {"name":"PRINCE EDWARD ISLAND",  "ll": [[50.10320164, "S"],[36.52150967 ,"E"]], "xyz":[ 0.51545596, 0.38171689, -0.76720099]}
    ]
}

Here is an image of the dymaxion coordinates projected onto Google Maps.

The dymaxion coordinates are (via a json load) loaded into a numpy array ‘dym_in’.
The order of the two definitions is not the same – so the mapping is (this may be wrong).

    i2d = [6, 4, 10, 0, 8, 9, 3, 2, 7, 5, 11, 1]
    # ico[i] is equivalent to dym_in[i2d[i]]
    dym = np.array([dym_in[m] for m in i2d ])

So now I have 12 normalised vertices in ‘ico’ and 12 dymaxion map vertices in ‘dym’, which are ordered such that ico[x] => dym[x].
I want to find the rotation (or approximate rotation) matrix that transforms ico to dym.

I say approximate, because the coordinates in given dym may not exactly mathematically define an icosahedron. I do not know because I do not know how to derive the transform!

What I know for sure is that the geoid is not relevant here – the Dymaxion starts from a spherical earth projection.

Likewise, I freely admit there may be bugs in my assumptions above.

What I want is to be able to derive the rotational matrix of any set of 12 icosahedral points from the initial golden-ratio starting set – bearing in mind that there are several 12(?) rotations to choose from, of course.