Assume we already have a shapely polygon of lat-lon coordinate system. We can project the coordinates to Albers Equal Area (aea) to calculate the area of the polygon.
from shapely.geometry import Polygon, shape
import pyproj#create a shapely polygon
poly = Polygon([(121,25), (121,24),(120,24)])# pa = pyproj.Proj("+proj=aea") # not work anymore.w,s,e,n = poly.bounds
pa = pyproj.Proj(proj='aea', lat_1 = s, lat_2 = n)#Extract the outer ring of the polygon
x, y = poly.exterior.coords.xy#projection
x, y = pa(list(x[:-1]), list(y[:-1]))#create a new geojson of the new projection polygon
cop = {"type": "Polygon", "coordinates": [zip(x, y)]}#geojson -> shapely polygon,
#then use the polygon's predefined function to calculate area.
area = shape(cop).area / 1000000 #km^2