Area of bottom polygons using obj.to_mesh post modifier but without applying it First of all, I'm quite new to python so please bear with me if I'm asking dumb questions. I've seen similar questions answered but still can't get where I want.

Say we have an active object in the scene, which is a polyline, which is extruded in z direction and altered with solidify modifier to give it some thickness (like a wall for example). I'm trying to get the summary area of all the bottom faces of the object as visible in the scene without actually applying the modifier. From what I've read it should be possible via obj.to_mesh, but while I can get the area of all faces with sum(f.calc_area(), I'm missing something (possibly obvious) when it comes to getting just the sum of bottom ones. 

Here is my attempt:

import bpy
import bmesh

context = bpy.context
scene = context.scene
obj = bpy.context.active_object

# --- get a mesh from the object ---
apply_modifiers = True
settings = 'PREVIEW'
mesh = obj.to_mesh(scene, apply_modifiers, settings)

bm = bmesh.new()
bm.from_mesh(mesh)

for face in bm.faces:

    face_location = face.calc_center_median()

    z = face_location[2]

    bpy.ops.object.mode_set(mode = 'EDIT')

    if z ==0: 

        face.select==True

    bpy.ops.object.mode_set(mode = 'OBJECT')

    area = sum(f.calc_area() for f in bm.faces if f.select)

    print(area)

bm.free()

bpy.data.meshes.remove(mesh)


What I get is value "0" printed out 16 times (for each face in the object I guess as it has 16 faces).
I realise it's probably something trivial I'm missing but maybe someone could please explain it to me?
I guess, you are getting zero because the face.select is not assigned value True, but is being checked for condition True (face.select == True). Also if you are interested only in bottom faces, i.e faces with all vertices having z coordinate value zero, then checking for median z coordinate may not work in all cases. There could be faces with median z at zero but vertices lying on either side of x-y plane.
I can think of the following alternative (I am open to critique and corrections as I am also comparatively new to blender :)

import bpy
import bmesh


bm = bmesh.new()
#bm = bmesh.from_edit_mesh(bpy.context.active_object.data) (without modifier)
bm.from_object(bpy.context.active_object, bpy.context.scene)

#Floats are not always absolute zero, so comparing with a tolerance
total_area = sum(f.calc_area() for f in bm.faces if len([v for v in f.verts if v.co[2] &lt; 0.0001 and v.co[2] &gt; -0.0001]) == len(f.verts)) 
print('total area:', total_area)

bm.free()