>>> import bmesh
>>> bm = bmesh.from_edit_mesh(C.object.data)
>>> bm.verts.ensure_lookup_table()
>>> verts = [v for v in bm.verts]
>>> print([v.index for v in verts])
[0, 1, 2, 3, 4]
>>> print([v.co.xy.to_tuple() for v in verts])
[(1.0, 1.5), (0.0, -2.0), (-1.0, -2.0), (-1.0, 1.5), (1.0, 0.0)]
>>> angles = [atan2(*v.co.yx) for v in verts]
>>> angles
[0.982793723247329, -1.5707963267948966, -2.0344439357957027, 2.158798930342464, 0.0]
>>> combined = zip(angles, verts)
>>> sorted_by_angle = list(sorted(combined))
>>> print([t[1].index + 1 for t in sorted_by_angle])
[3, 2, 5, 1, 4]
In 3D space you probably want to estimate the normal of the face (if all vertices lay on the same plane, that's just cross product of 2 vectors between vertices, otherwise perhaps you want to average the plane first, then calculate normal) and rotate the coordinate system so the normal points up...