How to snap object local rotation to its nearest global axis with python
Snap closest local to global rotation.

Example of use, object with some arbitrary rotation, can see the local axes from the widget and the global axis from navigation widget, rotates in this case such that Z axes match
The columns of the 3x3 rotation matrix of the world matrix of an object is its unit vector global axis alignment.
Find the minimum angle between any and the global axes.
Create a rotation matrix from the difference between the vectors with minimal axis and transform the original.
import bpy
from mathutils import Matrix
from bpy import context

ob = context.object
mw = ob.matrix_world.copy()
R = mw.to_3x3().col
I = Matrix.Identity(3)

diffs = [locl.angle(globl) for locl, globl in zip(R, I)]
idx = diffs.index(min(diffs))

M  = (Matrix.Translation(mw.translation) @ 
      R[idx].rotation_difference(I[idx]).to_matrix().to_4x4() @
      Matrix.Translation(-mw.translation)) 
ob.matrix_world = M @ ob.matrix_world