Calculating and exporting global rotation of faces created from splitting an isosphere
A script to do this


Shows origin and local orientation of one split off face


Adds an icosphere, edge splits, separates and returns to object mode
Selected objects are the result of separate modifier.
For each face object get its lone face and calculate the global center and normal.  See how far the normal is rotated from default quaternion.


Use this to create the global matrix aligned to normal, origin at center.  The object space inverse is applied to the mesh with transfrom (makes origin (0, 0, 0).

The global transform is given to the object by way of world matrix. 


Script.

import bpy
from mathutils import Quaternion

context = bpy.context
bpy.ops.mesh.primitive_ico_sphere_add(enter_editmode=True)
bpy.ops.mesh.edge_split()
bpy.ops.mesh.separate(type='LOOSE')

bpy.ops.object.mode_set()
z = Quaternion()
for ob in context.selected_objects:
    print(ob)
    me = ob.data
    mw = ob.matrix_world
    mwi = mw.inverted()
    f = ob.data.polygons[0]
    c = mw @ f.center
    n = mw @ f.normal

    q = z.rotation_difference(n.to_track_quat())
    M =  q.to_matrix().to_4x4()
    M.translation = c
    me.transform(mwi @ M.inverted())
    ob.matrix_world = M


Related. Another example of making a matrix based on known orthogonals.

https://blender.stackexchange.com/a/177331/15543