Is there a way to get a list of coplanar faces using Python?
Using bmesh

As an alternative to bpy.ops.mesh.select_similar(type='COPLANAR')


  4 points are coplanar if the volume created by the points is 0. If any
  three points determine a plane then additional points can be checked
  for coplanarity by measuring the distance of the points from the
  plane, if the distance is 0 then the point is coplanar


Select all faces with normals within angle tolerance of face.normal, then check if the convex hull,  made from the verts of each of these faces with those of selected face, volume is zero (or within some tolerance).  If so it must be on same plane.

Quick script to test.  Runs thru all mesh objects in the scene each face in the objects mesh and the coplanar list of the indices of faces coplanar to that face.

"Volume" of face

import bpy
import bmesh
from mathutils import Vector

def coplanar_list(bm, index, TOL=0.001):
    bm.faces.ensure_lookup_table()
    face = bm.faces[index]
    coplanar = []
    # narrow down selection based on normals
    faces = [f for f in bm.faces if f.normal.angle(face.normal) &lt; TOL]
    # list of faces with same normal
    bm2 = bmesh.new()
    for f in faces:
        verts = [v.co for v in face.verts]
        if f.index != face.index:
            verts.extend(v.co for v in f.verts)
        # make a dummy bmesh with verts of each face
        o = sum(verts, Vector((0,0,0))) / len(verts)
        for v in verts:
            bm2.verts.new(v - o)
        #calculate a convex hull
        bmesh.ops.convex_hull(bm2, input=bm2.verts)
        # if the convex hull has no volume its on same plane.
        volume = bm2.calc_volume()
        if volume &lt; 0.001:
            coplanar.append(f.index)
        bm2.clear()
    bm2.free()
    return coplanar

context = bpy.context
scene = context.scene
meshobjs = [o for o in scene.objects if o.type == 'MESH']
print("-" * 72)
for o in meshobjs:
    print("Object", o.name)
    if o == context.edit_object:
        bm = bmesh.from_edit_mesh(o.data)
    else:
        bm = bmesh.new()
        bm.from_mesh(o.data)
    for f in bm.faces:
        print("\tface %3d " % (f.index), coplanar_list(bm, f.index))


Or using distance point to plane for 4th and higher verts of plane defined by first 3 verts.  Note this script merely lists faces whose verts are coplanar by way of example, not those that are coplanar to active face.

import bpy
import bmesh

from mathutils.geometry import distance_point_to_plane as dp2p

TOL = 1e-6

context = bpy.context
ob = context.edit_object #  run in edit mode
me = ob.data

bm = bmesh.from_edit_mesh(me)

coplanar = []
for f in bm.faces:
    if len(f.verts) == 3: #  tri
        coplanar.append(f)
        continue
    v1 = f.verts[1].co - f.verts[0].co
    v2 = f.verts[2].co - f.verts[0].co
    if all(abs(dp2p(v.co, f.verts[0].co, v1.cross(v2))) &lt; TOL for v in f.verts[3:]) :
        coplanar.append(f)

print([f.index for f in coplanar])