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

mijujun

已于 2022-06-15 21:30:56 修改

阅读量3.3k

点赞数 5

分类专栏：学习记录文章标签： python 几何学

于 2022-06-15 21:17:46 首次发布

本文链接：https://blog.csdn.net/mijujun/article/details/125301553

版权

由于地理投影导致导致每个像元实际地面面积不同，越靠近北极实际面积越小，越靠近赤道实际面积越大，如果不进行面积加权就简单平均，会导致温度较实际温度偏低。
直接使用卫星地图的计算面积功能就会遇到这样的问题，多数卫星地图的计算面积功能是将地图假设为平面来计算，经纬度变化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[: