I am using Python to solve a problem with many polygons that could be really stressful with other software related to my field (mining), however i havent found a way to extract the Z coordinates from the Shapefile that I am working with.
You can use exterior and interior methods to extract the coordinates as tuples (x,y,z). I belive this should work for both polygon and multipolygons with any number of interior rings:
import geopandas as gpd df = gpd.read_file('/home/bera/GIS/test/zpoly.shp') def give_z(x): if x.type == 'Polygon': x = [x] zlist =  for polygon in x: zlist.extend([c[-1] for c in polygon.exterior.coords[:-1]]) for inner_ring in polygon.interiors: zlist.extend([c[-1] for c in inner_ring.coords[:-1]]) return zlist #return sum(zlist)/len(zlist) #In your case to get mean. Or just return zlist if they are all the same df['z'] = df.geometry.apply(give_z)