雷达系列 | 如何绘制极坐标下的雷达数据

用户11172986

发布于 2024-11-04 13:45:38

36000

代码可运行

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

运行总次数：0

代码可运行

雷达系列 | 如何绘制极坐标下的雷达数据

温馨提示

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

前言

一位读者朋友私信说不知道怎么处理极坐标下的雷达数据，那么我们今天来了解一下

项目目的

本项目旨在解决在气象作图过程中将雷达数据的极坐标转为经纬度的问题

需要注意的是，你必须知道雷达的坐标、方位角与库长

项目方法

azimuth_range_to_lat_lon 是 MetPy 库中的一个函数，用于将极坐标系统中的方位角和距离位置转换为经纬度坐标。该函数的参数和返回值如下：

参数

azimuths (array-like)：定义网格的一系列方位角。如果这不是一个 pint.Quantity 对象，则假定单位是度。
ranges (array-like)：从极点（即坐标系统的原点）到各点的距离数组。通常以米为单位。
center_lat (float)：极点的纬度，以十进制度数表示。
center_lon (float)：极点的经度，以十进制度数表示。
geod (pyproj.Geod or None, 可选)：用于前向方位角和距离计算的 PyProj Geod 对象。如果为 None，则使用默认的球形椭球体。

返回值

lon, lat (2D arrays)：与原始位置相对应的经度和纬度二维数组。

注意

这个函数对于处理雷达数据或任何其他以极坐标形式提供的地理空间数据非常有用，因为它允许用户将这些数据转换成更常见的经纬度格式，以便进行进一步的分析或可视化。在实际应用中，你可能需要安装 MetPy 和 pyproj 这两个库来使用这个功能。如果你的数据不是以 pint.Quantity 的形式提供，确保它们是以正确的单位（例如，方位角为度，距离为米）给出的。

导入库

import cartopy.crs as ccrs
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import numpy as np

from metpy.calc import azimuth_range_to_lat_lon
from metpy.cbook import get_test_data
from metpy.io import Level3File
from metpy.plots import add_timestamp, colortables, USCOUNTIES
from metpy.units import units

数据处理

# 打开雷达数据文件
name = '/home/mw/project/KOUN_SDUS54_N0QTLX_201305202016'
f = Level3File(name)

# 从文件对象中提取数据字典
datadict = f.sym_block[0][0]

# 根据文件指定的比例尺将数据转换为数组
data = f.map_data(datadict['data'])
print("雷达原始数据：",data)
# 获取方位角和距离，并为其分配单位
az = units.Quantity(np.array(datadict['start_az'] + [datadict['end_az'][-1]]), 'degrees')
rng = units.Quantity(np.linspace(0, f.max_range, data.shape[-1] + 1), 'kilometers')
print("方位角",az)
print("距离",rng)
# 从文件中提取中心经度和纬度
cent_lon = f.lon
cent_lat = f.lat

# 将方位角和距离转换为地理坐标
xlocs, ylocs = azimuth_range_to_lat_lon(az, rng, cent_lon, cent_lat)
print("经纬度数据：",xlocs)

雷达原始数据：[[nan nan 5.5 ... nan nan nan]
 [nan nan 4.  ... nan nan nan]
 [nan nan 5.  ... nan nan nan]
 ...
 [nan nan 3.5 ... nan nan nan]
 [nan nan 3.  ... nan nan nan]
 [nan nan 3.5 ... nan nan nan]]
方位角 [123.0 124.0 125.0 126.0 127.0 128.0 129.0 130.0 131.0 132.0 133.0 134.0 135.0 136.0 137.0 138.0 139.0 140.0 141.0 142.0 143.0 144.0 145.0 146.0 146.9 148.0 149.0 150.0 151.0 152.0 152.9 154.0 155.0 156.0 157.0 158.0 159.0 160.0 161.0 162.0 163.0 164.0 165.0 166.0 167.0 168.0 169.0 170.0 171.0 172.0 173.0 174.0 175.0 176.0 177.0 178.0 179.0 180.0 181.0 182.0 183.0 184.0 185.0 186.0 187.0 188.0 189.0 190.0 191.0 192.0 193.0 194.0 195.0 196.0 197.0 198.0 199.0 200.0 201.0 202.0 203.0 204.0 205.0 206.0 207.0 208.0 209.0 210.0 211.0 212.0 213.0 214.0 215.0 216.0 217.0 218.0 219.0 220.0 221.0 222.0 223.0 223.9 225.0 226.0 227.0 228.0 229.0 229.9 231.0 232.0 233.0 234.0 235.0 236.0 237.0 238.0 239.0 240.0 241.0 242.0 243.0 244.0 245.0 246.0 247.0 248.0 249.0 250.0 251.0 252.0 253.0 254.0 255.0 256.0 257.0 258.0 259.0 260.0 261.0 262.0 263.0 264.0 265.0 266.0 267.0 268.0 269.0 270.0 271.0 272.0 273.0 274.0 275.0 276.0 277.0 278.0 279.0 280.0 281.0 282.0 283.0 284.0 285.0 286.0 287.0 288.0 289.0 290.0 290.90000000000003 292.0 293.0 294.0 295.0 296.0 297.0 298.0 299.0 300.0 301.0 302.0 303.0 303.90000000000003 305.0 306.0 307.0 308.0 309.0 310.0 311.0 312.0 313.0 314.0 315.0 316.0 317.0 318.0 319.0 320.0 321.0 322.0 323.0 323.90000000000003 325.0 326.0 327.0 328.0 329.0 330.0 331.0 332.0 333.0 334.0 335.0 336.0 337.0 338.0 339.0 340.0 341.0 342.0 343.0 344.0 345.0 346.0 347.0 348.0 349.0 350.0 351.0 352.0 353.0 354.0 355.0 356.0 357.0 358.0 359.0 0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 12.0 13.0 14.0 15.0 16.0 17.0 18.0 19.0 20.0 21.0 22.0 23.0 24.0 25.0 26.0 27.0 28.0 29.0 30.0 31.0 32.0 33.0 34.0 35.0 36.0 37.0 38.0 39.0 40.0 40.900000000000006 42.0 43.0 44.0 45.0 46.0 47.0 48.0 49.0 50.0 51.0 52.0 53.0 54.0 55.0 56.0 57.0 58.0 59.0 60.0 61.0 62.0 63.0 64.0 65.0 66.0 67.0 68.0 69.0 70.0 71.0 72.0 73.0 74.0 75.0 76.0 77.0 78.0 79.0 80.0 81.0 82.0 83.0 84.0 85.0 86.0 87.0 88.0 89.0 90.0 91.0 92.0 93.0 94.0 95.0 96.0 97.0 98.0 99.0 100.0 101.0 102.0 102.9 104.0 105.0 106.0 107.0 108.0 109.0 110.0 111.0 112.0 113.0 114.0 115.0 116.0 117.0 118.0 119.0 120.0 121.0 122.0 123.0] degree
距离 [0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 12.0 13.0 14.0 15.0 16.0 17.0 18.0 19.0 20.0 21.0 22.0 23.0 24.0 25.0 26.0 27.0 28.0 29.0 30.0 31.0 32.0 33.0 34.0 35.0 36.0 37.0 38.0 39.0 40.0 41.0 42.0 43.0 44.0 45.0 46.0 47.0 48.0 49.0 50.0 51.0 52.0 53.0 54.0 55.0 56.0 57.0 58.0 59.0 60.0 61.0 62.0 63.0 64.0 65.0 66.0 67.0 68.0 69.0 70.0 71.0 72.0 73.0 74.0 75.0 76.0 77.0 78.0 79.0 80.0 81.0 82.0 83.0 84.0 85.0 86.0 87.0 88.0 89.0 90.0 91.0 92.0 93.0 94.0 95.0 96.0 97.0 98.0 99.0 100.0 101.0 102.0 103.0 104.0 105.0 106.0 107.0 108.0 109.0 110.0 111.0 112.0 113.0 114.0 115.0 116.0 117.0 118.0 119.0 120.0 121.0 122.0 123.0 124.0 125.0 126.0 127.0 128.0 129.0 130.0 131.0 132.0 133.0 134.0 135.0 136.0 137.0 138.0 139.0 140.0 141.0 142.0 143.0 144.0 145.0 146.0 147.0 148.0 149.0 150.0 151.0 152.0 153.0 154.0 155.0 156.0 157.0 158.0 159.0 160.0 161.0 162.0 163.0 164.0 165.0 166.0 167.0 168.0 169.0 170.0 171.0 172.0 173.0 174.0 175.0 176.0 177.0 178.0 179.0 180.0 181.0 182.0 183.0 184.0 185.0 186.0 187.0 188.0 189.0 190.0 191.0 192.0 193.0 194.0 195.0 196.0 197.0 198.0 199.0 200.0 201.0 202.0 203.0 204.0 205.0 206.0 207.0 208.0 209.0 210.0 211.0 212.0 213.0 214.0 215.0 216.0 217.0 218.0 219.0 220.0 221.0 222.0 223.0 224.0 225.0 226.0 227.0 228.0 229.0 230.0 231.0 232.0 233.0 234.0 235.0 236.0 237.0 238.0 239.0 240.0 241.0 242.0 243.0 244.0 245.0 246.0 247.0 248.0 249.0 250.0 251.0 252.0 253.0 254.0 255.0 256.0 257.0 258.0 259.0 260.0 261.0 262.0 263.0 264.0 265.0 266.0 267.0 268.0 269.0 270.0 271.0 272.0 273.0 274.0 275.0 276.0 277.0 278.0 279.0 280.0 281.0 282.0 283.0 284.0 285.0 286.0 287.0 288.0 289.0 290.0 291.0 292.0 293.0 294.0 295.0 296.0 297.0 298.0 299.0 300.0 301.0 302.0 303.0 304.0 305.0 306.0 307.0 308.0 309.0 310.0 311.0 312.0 313.0 314.0 315.0 316.0 317.0 318.0 319.0 320.0 321.0 322.0 323.0 324.0 325.0 326.0 327.0 328.0 329.0 330.0 331.0 332.0 333.0 334.0 335.0 336.0 337.0 338.0 339.0 340.0 341.0 342.0 343.0 344.0 345.0 346.0 347.0 348.0 349.0 350.0 351.0 352.0 353.0 354.0 355.0 356.0 357.0 358.0 359.0 360.0 361.0 362.0 363.0 364.0 365.0 366.0 367.0 368.0 369.0 370.0 371.0 372.0 373.0 374.0 375.0 376.0 377.0 378.0 379.0 380.0 381.0 382.0 383.0 384.0 385.0 386.0 387.0 388.0 389.0 390.0 391.0 392.0 393.0 394.0 395.0 396.0 397.0 398.0 399.0 400.0 401.0 402.0 403.0 404.0 405.0 406.0 407.0 408.0 409.0 410.0 411.0 412.0 413.0 414.0 415.0 416.0 417.0 418.0 419.0 420.0 421.0 422.0 423.0 424.0 425.0 426.0 427.0 428.0 429.0 430.0 431.0 432.0 433.0 434.0 435.0 436.0 437.0 438.0 439.0 440.0 441.0 442.0 443.0 444.0 445.0 446.0 447.0 448.0 449.0 450.0 451.0 452.0 453.0 454.0 455.0 456.0 457.0 458.0 459.0 460.0] kilometer
经纬度数据：[[-97.278      -97.26875527 -97.25951166 ... -93.15820683 -93.14945551
  -93.14070523]
 [-97.278      -97.26886147 -97.25972408 ... -93.20829874 -93.19965936
  -93.19102103]
 [-97.278      -97.26897045 -97.25994206 ... -93.25952548 -93.25100023
  -93.24247603]
 ...
 [-97.278      -97.26855135 -97.25910378 ... -93.06149175 -93.05252325
  -93.04355576]
 [-97.278      -97.26865189 -97.25930488 ... -93.10926587 -93.10040483
  -93.09154482]
 [-97.278      -97.26875527 -97.25951166 ... -93.15820683 -93.14945551
  -93.14070523]]

绘图

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy
# 设置图形布局
fig = plt.figure(figsize=(7.5, 8))  

# 定义反射率的色表和数据范围
ctable = ('NWSStormClearReflectivity', -20, 0.5)  # dBZ


# 创建子图
crs = ccrs.LambertConformal()
ax = fig.add_subplot(1, 1, 1, projection=crs)

# 获取色表和归一化器
norm, cmap = colortables.get_with_steps(*ctable)

# 绘制反射率数据
mesh = ax.pcolormesh(xlocs, ylocs, data, norm=norm, cmap=cmap, transform=ccrs.PlateCarree())

# 设置地图范围
ax.set_extent([cent_lon - 1.5, cent_lon + 1.5, cent_lat - 1.5, cent_lat + 1.5])

# 设置纵横比
ax.set_aspect('equal', 'datalim')

# 添加经纬度刻度
gl = ax.gridlines(
    crs=ccrs.PlateCarree(),  # 指定网格线的坐标系
    draw_labels=True,  # 显示标签文字
    linewidth=1,  # 网格线宽度
    color="gray",  # 网格线颜色
    alpha=0.5,  # 图形内网格的透明度
    linestyle="--",  # 网格线的线型
    x_inline=False,  # 禁止x标签显示在图框内部
    y_inline=False,  # 禁止y标签显示在图框内部
    rotate_labels=False,  # 禁止标签文字旋转
)

# 关闭顶部和右侧的标签
gl.top_labels = False  # 关闭顶部标签
gl.right_labels = False  # 关闭右侧标签

# 设置经度和纬度的格式
gl.xformatter = cartopy.mpl.gridliner.LONGITUDE_FORMATTER
gl.yformatter = cartopy.mpl.gridliner.LATITUDE_FORMATTER

# 设置标签字体大小
gl.xlabel_style = {'size': 12}  # 经度标签大小
gl.ylabel_style = {'size': 12}  # 纬度标签大小
# 添加时间戳
add_timestamp(ax, f.metadata['prod_time'], y=0.02, high_contrast=True)


# 显示图形
plt.show()