如何用经纬度计算多边形的面积？

import math

# 经纬度转墨卡托投影
def lonlat_to_mercator(lon, lat):
    r_major = 6378137.0  # 地球的半长轴
    x = r_major * math.radians(lon)
    scale = x / lon
    y = 180.0 / math.pi * math.log(math.tan(math.pi / 4.0 + lat * (math.pi / 180.0) / 2.0)) * scale
    return x, y

# 计算多边形的面积
def calculate_polygon_area(points):
    area = 0.0
    if len(points) < 3:
        return area

    for i in range(len(points)):
        p1 = points[i]
        p2 = points[(i + 1) % len(points)]
        area += (p2[0] - p1[0]) * (p2[1] + p1[1])

    return abs(area / 2.0)

# 示例数据，多边形的顶点坐标（经纬度）
polygon = [(116.397, 39.917), (116.397, 39.947), (116.437, 39.947), (116.437, 39.917)]

# 将经纬度转换为墨卡托投影坐标
polygon_mercator = [lonlat_to_mercator(lon, lat) for lon, lat in polygon]

# 计算多边形的面积
area = calculate_polygon_area(polygon_mercator)

print("多边形的面积为：", area, "平方米")