再说不会用python计算地球表面多边形面积,可不能了!(记录五种可行方法)

这篇具有很好参考价值的文章主要介绍了再说不会用python计算地球表面多边形面积,可不能了!(记录五种可行方法)。希望对大家有所帮助。如果存在错误或未考虑完全的地方,请大家不吝赐教,您也可以点击"举报违法"按钮提交疑问。

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

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

分类方法

目前经过测试可行的方法共有五种,五种方法均假设地球为球体。我给他们分了三个大类,包括: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处理气象数据(六)数据的面积加权

四、总结

前四种方法适用于地球表面任意多边形的情况,其中使用积分的方法三精度最高,推荐所有人使用,方法二最简便,推荐给追求整洁的小伙伴使用。第五种方法仅适用于四边均与经纬线平行的区域,推荐给需要计算模式画好的网格面积的小伙伴。文章来源地址https://www.toymoban.com/news/detail-427399.html

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

到了这里,关于再说不会用python计算地球表面多边形面积,可不能了!(记录五种可行方法)的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处: 如若内容造成侵权/违法违规/事实不符,请点击违法举报进行投诉反馈,一经查实,立即删除!

领支付宝红包 赞助服务器费用

相关文章

  • 【计算机图形学 】扫描线多边形填充算法 | OpenGL+鼠标交互

    传送门 实现多边形扫描线填充算法,并和鼠标进行交互。 具体原理略过,会贴上完整代码,可直接运行。 环境: vs2019,OpenGL的库(可以搜索如何用vs使用OpenGL的库,可以使用vs自带的插件或者其他方法,很方便) 要点: 1.NET和AET的创建,改动 2.改变鼠标点击和鼠标拖拽的响应

    2023年04月08日
    浏览(73)
  • Elasticsearch计算距离,根据距离排序,地理点和地理多边形范围查找

    Elasticsearch 计算并返回距离一共有两种方法: sort 和 script_fields CentOS 7.6 Elasticsearch 7.10 响应结果如下, hits 下的 sort 字段就是距离,单位:km。 5.x 以前支持:distanceInKm(lat, lon) 函数,后来被废弃。现在只支持 arcDistance(lat, lon) 函数:计算两点距离,单位为:m。响应结果如下,

    2024年02月09日
    浏览(44)
  • 【Weiler-Atherton算法】 计算机图形学多边形裁剪算法

    源代码: https://github.com/ricar0/Weiler-Atherton-Alogrithm/tree/master 通常来说就是利用多边形来裁剪多边形的一种方法,一般情况下是利用矩形来裁剪凹凸多边形 凸多边形 凹多边形 上面红色划线部分就是裁剪出的部分 OPENGL基础语法 基本上就是一些画线和画多边形的操作,难度较低

    2023年04月09日
    浏览(63)
  • Python:opencv画点、圆、线、多边形、矩形

    简介 :机器学习视觉方向一般都需要在图像中添加标注框,标注框有着很大的用处,特别是对图像中某些需要关注的特征起到圈定的效果,方便对特征选择进行处理。 相关攻略: 机器学习:基本流程 Python:调用摄像头使用cv2库录制视频 Python:视频拆分成一帧一帧的图片

    2024年02月04日
    浏览(86)
  • python opencv 绘制矩形、圆、线、多边形

    👨‍💻 个人简介: 深度学习图像领域工作者 🎉 总结链接:              链接中主要是个人工作的总结,每个链接都是一些常用demo,代码直接复制运行即可。包括:                     📌 1.工作中常用深度学习脚本                     📌 2.to

    2024年02月03日
    浏览(102)
  • Python OpenCV实现鼠标绘制矩形框和多边形

    目录 Python OpenCV实现鼠标绘制矩形框和多边形 1. OpenCV鼠标事件操作说明 (1)setMouseCallback函数说明 (2)回调函数onMouse说明 (3)event 具体说明: (4)flags 具体说明 2. OpenCV实现鼠标绘制矩形框和多边形框 (1)绘制矩形框 (2)绘制多边形 (3)键盘控制 3. 完整的代码 本篇将

    2024年02月06日
    浏览(72)
  • Python的海龟 turtle 库使用详细介绍(画任意多边形,全网最详细)

    学Turtle库,其实就是学数学,而且还能提高对数学和学习的兴趣。Turtle库还能够帮助孩子更好地理解几何学和数学概念,比如角度、比例、几何图形的性质等等,是Python中一个很有趣的库。 Turtle库是Python中一个很有趣的库,可以用来绘制各种图形,比如直线、圆、正方形等等

    2024年04月13日
    浏览(59)
  • 计算机图形学实验——利用MFC对话框实现多边形绘制与填充(扫描线填充算法)附源码

    内容概括: 利用基于对话框的MFC项目 实现鼠标点击绘制多边形 实现扫描线算法填充多边形 源码见Yushan-Ji/ComputerGraphics: ECNU2023秋 计算机图形学课程实验代码 (github.com) 通过鼠标交互输入多边形 对各种多边形进行填充,包括边界自交的情况 利用 OnLButtonDown 和 OnRButtonDown 函数,

    2024年02月04日
    浏览(70)
  • [python][pcl]python-pcl案例之为平面模型构造凹凸外壳多边形

    测试环境: pcl==1.12.1 python-pcl==0.3.1 python==3.7 代码: 运行结果: PointCloud after filtering has: 139656 data points. PointCloud after segmentation has: built-in method count of list object at 0x0000028CA3429A48 inliers. PointCloud after projection has: 139656 data points. Concave hull has: 281 data points. table_scene_mug_stereo_textured.p

    2024年02月11日
    浏览(44)
  • python 如何判断点是否在多边形(三角形)内,或求点在3D面上的投影?

    方法1: 用shapely中的geometry包 1)polygon.covers(point) 如果point在多边形polygon上(包括边),返回True,否则False。 2)polygon.contains(point) 如果point在多边形polygon上(不包括边),返回True,否则False。 方法2: 用blender的内置python api。 将点投影到三角形平面上,并检查其是否在三角形

    2023年04月09日
    浏览(45)

觉得文章有用就打赏一下文章作者

支付宝扫一扫打赏

博客赞助

微信扫一扫打赏

请作者喝杯咖啡吧~博客赞助

支付宝扫一扫领取红包,优惠每天领

二维码1

领取红包

二维码2

领红包