博士来稿！如何计算自定义的风暴面积

用户11172986

发布于 2024-06-20 19:08:14

17000

代码可运行

文章被收录于专栏：气python风雨气python风雨

运行总次数：0

代码可运行

前言

前段时间正苦于没有素材，正好有位博士老兄投稿，让我们看看有什么新花样吧

台风是经久不衰的研究课题，我们经常会找一些指标去描述台风的强度大小。

而这位老哥的导师想要以台风的某条闭合等值线为准，计算其包围的面积。

这怎么搞呢？他找到一个使用polygon计算matplotlib绘图对象面积的方法

下面让我们开始吧

⏰ 温馨

由于可视化代码过长隐藏，可点击运行Fork查看 🔜🔜若没有成功加载可视化图，点击运行可以查看 ps:隐藏代码在【代码已被隐藏】所在行，点击所在行，可以看到该行的最右角，会出现个三角形，点击查看即可

当然，这位读者的数据已丢失，我们用台风模拟的wrfout数据替代一下画一个图。

相当位温计算与绘图代码可参考如何计算WRF台风模拟的假相当位温

假定我们要计算的区域是370k以内的区域，那么我们假定它就是核心的风暴区域。下面开始计算

获取路径

In [8]:

fig = plt.figure(figsize=(15, 12))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
# Plot Titles
plt.title(r'850-hPa Heights (m) 相当位温', loc='left')
ax.set_extent([116., 130., 20.,32.], ccrs.PlateCarree())
# Plot Background
add_china_map_2basemap(ax, name="river", edgecolor='k', lw=0.5, encoding='gbk', grid0=None)  # "河流"
add_china_map_2basemap(ax, name="nation",edgecolor='k', lw=0.5, encoding='gbk', grid0=None)  # "国界"
add_china_map_2basemap(ax, name="province", edgecolor='k', lw=0.5, encoding='gbk', grid0=None)  # "省界"
cs1 = ax.contour(lons, lats, se2,levels=[372],colors='k',transform=ccrs.PlateCarree())

path0 = cs1.collections[0].get_paths()[0]
vertices = path0.vertices
plt.plot(vertices[:,0],vertices[:,1])

Out[8]:

[<matplotlib.lines.Line2D at 0x7f52f9e68400>]

# 获取第一条曲线的路径数据    
path0 = cs1.collections[0].get_paths()[0]    

# 提取路径的顶点坐标    
vertices = path0.vertices    

# 使用matplotlib绘制顶点的横纵坐标    
plt.plot(vertices[:,0], vertices[:,1])

这段代码的主要目的是提取一条曲线的路径顶点坐标，并使用matplotlib将其绘制出来。首先从cs1对象（可能是由matplotlib绘制的图形或图像）的首个collections元素中获取第一条路径数据。接着，从路径数据中提取出所有顶点的坐标，保存在二维数组vertices中。最后，利用matplotlib的plot函数，以顶点的横坐标和纵坐标为输入，绘制出表示这些顶点连线的折线图。

投影转换与面积计算

In [19]:

import cartopy.crs as ccrs

# 定义投影方式：兰伯特等角投影
# 参数设置：
#   central_longitude=122：中央经度为122度
#   central_latitude=28：中央纬度为28度
#   standard_parallels=(22, 23)：标准纬线为26度和28度
projection = ccrs.LambertConformal(
    central_longitude=122,
    central_latitude=28,
    standard_parallels=(26, 28)
)

# 将地理坐标（经度122度，纬度28度）从PlateCarree投影（经纬度投影）转换到当前定义的LambertConformal投影
projection.transform_point(122, 28, src_crs=ccrs.PlateCarree())

# 将之前提取的顶点坐标（vertices）从PlateCarree投影转换到LambertConformal投影
xyz = projection.transform_points(ccrs.PlateCarree(), vertices[:, 0], vertices[:, 1])

# 引入shapely库中的Polygon类，用于构建多边形几何对象
from shapely.geometry import Polygon

# 使用转换后的坐标xyz创建一个shapely多边形
polygon = Polygon(xyz[:, 0:2])

# 计算多边形的面积
area = polygon.area

# 打印多边形面积，保留三位小数，科学记数法表示
print('{:4.3e}'.format(area))