There's probably a cleaner solution, but for the time being, you can use:

import bpy
from math import radians
from mathutils import Vector, Matrix

obj = bpy.context.active_object
mesh = obj.data


def get_layers():
    """
    return: tuple of 20 booleans. True if layer is active layer,False if not
    Allows the mesh to be created on the active layer. There is only 1 active layer
    """
    active_layer = bpy.context.scene.active_layer
    layers = [False] * 20
    layers[active_layer] = True
    return tuple(layers)


def create_light_at_face():
    """
    layers: tuple. 20 booleans where exaclty one boolen is True
    return: None
    creates an area light the has the same orientation as the normal at the center of each face
    """
    mat_world = obj.matrix_world
    up = Vector((0,0,1))
    
    for poly in mesh.polygons:
        co = mat_world * Vector(poly.center)
        
        forward = poly.normal.copy()
        forward.rotate(mat_world)
        right = forward.cross(up).normalized() # Vector.length closer to 1.0
        up = right.cross(forward)
    
        rot = Matrix((right, up, -forward)).transposed().normalized().to_4x4()
        mat = Matrix.Translation(co) * rot
    
        bpy.ops.object.lamp_add(type='AREA', view_align=False, location=co, layers=get_layers())
        # scale lamp up a little fo we can see the area
        bpy.context.object.data.size = 3
        bpy.context.object.matrix_world = mat

create_light_at_face()

Also tried eulers, but did not work reliably (some lamps pointing inwards).