How to calculate the angles between axis and bounding box normals?
The orientation of an object is obj.rotation_euler or obj.rotation_quaternion (and whichever one does not match the object's rotation_mode can be invalid or inaccurate).

The orientation can also be recovered from obj.matrix_world.decompose()[1], but that's a little bit of overkill.

For a slightly deeper dive down the rabbit hole, check out Relationship between global vertex coordinates and mesh object&#39;s matrix attributes .