random place on a face with python
Use the local space.
The default plane in blender has local coordinates where x and y range from -1 to 1 and z is zero.
&gt;&gt;&gt; me = C.object.data
&gt;&gt;&gt; for v in me.vertices:
...     v.co[:]
...     
(-1.0, -1.0, 0.0)
(1.0, -1.0, 0.0)
(-1.0, 1.0, 0.0)
(1.0, 1.0, 0.0)

With this arrangement can consider only the x, y coordinates effectively making it 2D within the plane z=0.
Hence any point with x and y randomly chosen from that range will be on the plane.
Simple test script, with an object that is the default plane &quot;Plane&quot; with any arbitrary global transform, and it will add an empty on its surface.
Multiplying this local point by the objects matrix world will ensure it is at a global location on the surface of planes current position.
import bpy
from mathutils import Vector, Matrix                    
from random import uniform
context = bpy.context
scene = context.scene

plane = scene.objects.get(&quot;Plane&quot;)
if plane:
    vec = Vector((uniform(-1, 1), uniform(-1, 1), 0))
    bpy.ops.object.empty_add(
            location = plane.matrix_world @ vec,
            ) 

Will re-iterate this relies on the plane having its original local or mesh coordinates in range explained above, making it a simple task to place an empty on its surface.
Via Parenting.
Another way to use local coordinates is via parenting. In this case need only set the local location of the empty within range above to have it on the plane. This will have the effect of aligning the axes of empty to plane.
if plane:
    vec = Vector((uniform(-1, 1), uniform(-1, 1), 0))
    bpy.ops.object.empty_add()
    mt = context.object
    mt.parent = plane
    mt.location = vec 

To make the influence of the parent planes scale to empty unity for any subsequent children of the empty,.   Does a child object inherit the matrix from the parent?
if plane:
    mw = plane.matrix_world
    _, _, scale = mw.decompose()
    S = Matrix.Diagonal(scale.to_4d())
    vec = Vector((uniform(-1, 1), uniform(-1, 1), 0))
    bpy.ops.object.empty_add()
    mt = context.object
    mt.parent = plane
    mt.matrix_parent_inverse = S.inverted()
    mt.location = S @ vec 

For an arbitrarily un axis aligned plane it is a bit more complex. https://stackoverflow.com/questions/42352622/finding-points-within-a-bounding-box-with-numpy