Position an object's origin point to the face below (python)?
This likely needs more testing but should do what you want. This is currently set up for explicitly moving objects on a 'work layer' to the selected faces of a target object in the z direction. Currently if multiple faces are capable of being snapped to the max movement limit should help to limit but be sure to test before saving your file.

import bpy
from bpy.types import Operator
from bpy.types import Panel
from mathutils import Vector
from mathutils import geometry

bl_info = {
    "name": "Blank1",
    "author": "Your Name Here",
    "version": (1, 0),
    "blender": (2, 79, 0),
    "location": "View3D &gt; Tools panel 'T' &gt; Blank1 &gt; Blank2",
    "description": "",
    "warning": "",
    "wiki_url": "",
    "category": "",
    }


def makeprop():
    bpy.types.Scene.my_tolerance = bpy.props.FloatProperty(
        name="bpy.types.Scene.my_tolerance",
        description="",
        default=5.0,
        )
    bpy.types.Scene.my_worklayer = bpy.props.IntProperty(
        name="bpy.types.Scene.my_worklayer",
        description="",
        default=0,
        min=0,
        max=19,
        )


def delprop():
    del bpy.types.Scene.my_tolerance
    del bpy.types.Scene.my_worklayer


def get_index_selected():
    obj = bpy.context.edit_object
    me = obj.data
    myfaces = []
    if obj.type not in ['MESH']:
            return myfaces
    for f in me.polygons:
        if f.select:
            myfaces.append(f.index)
    return myfaces


def GetTargetVert(target, object, tol):
    me = target.data
    bpy.ops.object.editmode_toggle()
    bpy.ops.object.editmode_toggle()
    my_target_faces = get_index_selected()
    for face in my_target_faces:
        x_max = -1000000.0
        x_min = 1000000.0
        y_max = -1000000.0
        y_min = 1000000.0
        face_verts = []
        for v in me.polygons[face].vertices:
            face_verts.append(v)
        bpy.ops.object.editmode_toggle()
        bpy.ops.object.editmode_toggle()
        W_mat = target.matrix_world
        my_line_a = object.location
        my_line_b = object.location + Vector((0, 0, -5))
        vert_index = me.polygons[face].vertices[0]
        my_plane_co = W_mat * me.vertices[vert_index].co
        my_plane_no = W_mat * me.polygons[face].normal
        for v in face_verts:
            l_co = me.vertices[v].co
            loc = W_mat * l_co
            if loc[0] &gt; x_max:
                x_max = loc[0]
            if loc[0] &lt; x_min:
                x_min = loc[0]
            if loc[1] &gt; y_max:
                y_max = loc[1]
            if loc[1] &lt; y_min:
                y_min = loc[1]
        print(
            'x_max %.2f x_min %.2f y_max %.2f y_min %.2f' % (
                x_max, x_min, y_max, y_min))
        print(
            'obj_x %.2f obj_y %.2f' % (object.location[0], object.location[1]))
        if ((x_min &lt; object.location[0] &lt; x_max) and
            (y_min &lt; object.location[1] &lt; y_max)):
                collision = geometry.intersect_line_plane(
                    my_line_a, my_line_b, my_plane_co, my_plane_no)
        else:
            collision = None
#            print('not in x,y range')
        if collision is None:
            print('no collision available between %s Face[%i] and %s'
                  % (target.name, face, object.name))
        elif abs(collision[2]-object.location[2]) &lt; tol:
            print('collision found between %s Face[%i] and %s'
                  % (target.name, face, object.name))
            print(collision)
            object.location = collision
        else:
            print('collision out of range between %s Face[%i] and %s'
                  % (target.name, face, object.name))
#            print(collision)


def align_objs(context):
    print("-"*20 + "align objects button pressed" + "-"*20)
    scene = context.scene
    target = scene.objects.active
    tol = scene.my_tolerance
    my_layer_objs = [
        ob for ob in scene.objects if ob.layers[scene.my_worklayer]]
    for obj in my_layer_objs:
        if (not obj == target) and (obj.type == 'MESH'):
            GetTargetVert(target, obj, tol)


class AlignOperator(bpy.types.Operator):
    """Tooltip"""
    bl_idname = "object.align_operator"
    bl_label = "Align Operator"

    @classmethod
    def poll(cls, context):
        return (context.mode == 'EDIT_MESH')

    def execute(self, context):
        align_objs(context)
        return {'FINISHED'}


class VIEW3D_PT_Solve_Intersect(Panel):
    bl_space_type = 'VIEW_3D'
    bl_region_type = 'TOOLS'
    bl_category = "Intersect Face"
    bl_label = "Intersect"

    def draw(self, context):

        scene = context.scene
        layout = self.layout
        row = layout.row()
        col = layout.column(align=True)
        box = col.box()
        box.label(text='Move target &amp; collision objects to work layer')
        box.label(text='Active Layer : %s' % scene.active_layer)
        box.prop(scene, "my_worklayer", "Work Layer")
        box.label(text="Set the max movement")
        box.prop(scene, "my_tolerance", "")
        box.label(text="Enter edit mode &amp; select target faces")
        box.label(text="Press the align objects button")
        box.operator("object.align_operator", text="align objects")
        box.label(text="Collision info available in console")


def register():
    bpy.utils.register_class(VIEW3D_PT_Solve_Intersect)
    bpy.utils.register_class(AlignOperator)
    makeprop()


def unregister():
    bpy.utils.unregister_class(VIEW3D_PT_Solve_Intersect)
    bpy.utils.unregister_class(AlignOperator)
    delprop()


if __name__ == "__main__":
    register()



Edited:
removed use of bmesh
corrected target face vertex positions to reference world position
additional check for collision withing the target face min max x,y
coordinates