再说不会用python计算地球表面多边形面积，可不能了！（记录五种可行方法）

由于地理投影导致导致每个像元实际地面面积不同，越靠近北极实际面积越小，越靠近赤道实际面积越大，如果不进行面积加权就简单平均，会导致温度较实际温度偏低。
直接使用卫星地图的计算面积功能就会遇到这样的问题，多数卫星地图的计算面积功能是将地图假设为平面来计算，经纬度变化1度时默认距离变化为10km。带来极大误差。使用谷歌卫星地图截取的（110，39），（115，40），（110，41）三个点之间的区域面积，后文中的坐标也是该坐标下进行的计算。

学习地球科学的小伙伴经常会因为无法计算网格的面积大小或者其他更加不规则的区域的面积而让头发惨死在手中，因此今天对计算地球表面多边形面积的方法进行了整理。

分类方法

目前经过测试可行的方法共有五种，五种方法均假设地球为球体。我给他们分了三个大类，包括：1.使用首尾闭环坐标（GeoJSON格式）；2.使用非闭环坐标；3.使用最大最小经纬度（仅能计算经纬度与地球网格平行的情况）。

一、使用首尾闭环坐标（GeoJSON格式坐标）

使用GeoJSON格式坐标的方法有两种，一种使用basemap包进行计算，一种使用area包。适合使用卫星云图直接导出坐标。

方法一：

import numpy
from mpl_toolkits.basemap import Basemap

#############   数据
coordinates=numpy.array([[110, 39],
                         [115, 40],
                         [110, 41],
                         [110, 39]])######### 坐标需要为闭环

############   计算
lats=coordinates[:,1]
lons=coordinates[:,0]

lat1=numpy.min(lats)
lat2=numpy.max(lats)
lon1=numpy.min(lons)
lon2=numpy.max(lons)

bmap=Basemap(projection='cea',llcrnrlat=lat1,llcrnrlon=lon1,urcrnrlat=lat2,urcrnrlon=lon2)
xs,ys=bmap(lons,lats)

area=numpy.abs(0.5*numpy.sum(ys[:-1]*numpy.diff(xs)-xs[:-1]*numpy.diff(ys)))
S=area/1e6

############   输出结果
print("面积为"+str(S)+"平方公里")

输出结果为

方法一计算面积为47355.612488090286平方公里

方法二：

使用函数area，未安装的筒子使用

pip install area

安装即可

from area import area

#############   数据
obj = {'type':'Polygon','coordinates':[[[110, 39],
                                        [115, 40],
                                        [110, 41],
                                        [110, 39]]]}########   闭环
############   计算
S3 = area(obj)/1e6

############   输出结果
print("方法二计算面积为"+str(S3)+"平方公里")

输出结果为

方法二计算面积为47461.8151872328平方公里

二、使用非闭环坐标（自定义函数计算）

使用两种思路定义函数，第一种函数使用积分方法，第二种函数将纬度乘以一个纬度的长度，然后将经度乘以纬度的长度和纬度的余弦。

方法三：

def polygon_area(lons,lats, radius = 6378137):
    """
    Computes area of spherical polygon, assuming spherical Earth.
    Returns result in ratio of the sphere's area if the radius is specified.
    Otherwise, in the units of provided radius.
    lats and lons are in degrees.
    """
    import numpy as np
    from numpy import arctan2, cos, sin, sqrt, pi, power, append, diff, deg2rad
    lons = np.deg2rad(lons)
    lats = np.deg2rad(lats)

    # Line integral based on Green's Theorem, assumes spherical Earth

    #close polygon
    if lats[0]!=lats[-1]:
        lats = append(lats, lats[0])
        lons = append(lons, lons[0])

    #colatitudes relative to (0,0)
    a = sin(lats/2)**2 + cos(lats)* sin(lons/2)**2
    colat = 2*arctan2( sqrt(a), sqrt(1-a) )

    #azimuths relative to (0,0)
    az = arctan2(cos(lats) * sin(lons), sin(lats)) % (2*pi)

    # Calculate diffs
    # daz = diff(az) % (2*pi)
    daz = diff(az)
    daz = (daz + pi) % (2 * pi) - pi

    deltas=diff(colat)/2
    colat=colat[0:-1]+deltas

    # Perform integral
    integrands = (1-cos(colat)) * daz

    # Integrate
    area = abs(sum(integrands))/(4*pi)

    area = min(area,1-area)
    if radius is not None: #return in units of radius
        return area * 4*pi*radius**2
    else: #return in ratio of sphere total area
        return area

lon = [110,115,110]    ######  经度
lat = [39,40,41]       ######  纬度

############   调用函数
S = polygon_area(lon,lat)
S1 = S/(1e6)

############   输出结果
print("方法三计算面积为"+str(S1)+"平方公里")

输出结果为

方法三计算面积为47262.22164026087平方公里

方法四：

def area_of_polygon(longitude,latitude):

    from math import pi, cos, radians
    earth_radius = 6371009 # in meters
    lat_dist = pi * earth_radius / 180.0

    y = [lat * lat_dist for lat in latitude]
    x = [long * lat_dist * cos(radians(lat))
                for lat, long in zip(latitude, longitude)]
    area = 0.0
    for i in range(-1, len(x) - 1):
        area += x[i] * (y[i + 1] - y[i - 1])
    return abs(area) / 2.0


#############   数据
lon = [110,115,110]
lat = [39,40,41]

############   调用函数
S4 = area_of_polygon(lon,lat)/1e6

############   输出结果
print("方法四计算面积为"+str(S4)+"平方公里")

输出结果为

方法四计算面积为47516.878614105466平方公里

以上四种方法参考了

码农家园——关于几何：如何使用python计算地球表面上多边形的面积？

三、使用最大最小经纬度（仅能计算经纬度与地球网格平行的情况）

方法五：

def calc_area_5(lat1,lat2,lon1,lon2):
    import numpy as np
    pi = 3.1415926
    R = 6371007.181

    area = (pi/180.0)*R*R*abs(np.sin(lat1/180.0*pi) - np.sin(lat2/180.0*pi)) * abs(lon1 - lon2)
    return area

### 设置网格坐标
lat1 = 39
lat2 = 41
lon1 = 110
lon2 = 115

############   调用函数
S5 = calc_area_5(lat1,lat2,lon1,lon2)/1e6

############   输出结果
print("方法五计算面积为"+str(S5)+"平方公里")

输出结果为

方法五计算面积为94711.52539324109平方公里

该方法被应用于许多模式中计算网格面积，如计算水循环的vic模型。优点在于计算简便，缺点也很明显，仅能计算四边均与经纬线平行的区域。

该方法引用于
Matlab处理气象数据（六）数据的面积加权

四、总结

前四种方法适用于地球表面任意多边形的情况，其中使用积分的方法三精度最高，推荐所有人使用，方法二最简便，推荐给追求整洁的小伙伴使用。第五种方法仅适用于四边均与经纬线平行的区域，推荐给需要计算模式画好的网格面积的小伙伴。

PS：所有以上方法均在假设地球为球体的条件下适用，若需要计算椭球体情况下的面积，五种方法都有误差。