如何用Python创建可在Google Earth中打开的KML格式风玫瑰图

目前没有直接生成Google Earth兼容KML风玫瑰的工具包，但你可以结合已掌握的windrose和simplekml模块，通过提取风玫瑰统计数据→转换为KML多边形→匹配样式三步实现，具体如下：

1. 从windrose提取风玫瑰核心数据 #

windrose生成的可视化图背后包含风向区间、各风速段频率（对应扇形半径）等关键数值，你需要从WindroseAxes对象中提取这些数据，而非仅生成matplotlib图：

from windrose import WindroseAxes
import numpy as np

# 替换为你的实际风速(ws)、风向(wd)数据
ws = np.random.random(1000) * 10  # 风速示例
wd = np.random.random(1000) * 360  # 风向示例（0-360度，北为0）

# 初始化风玫瑰轴并生成统计数据
ax = WindroseAxes.from_ax()
# 设置和示例图一致的风速区间、扇区开口等参数
ax.bar(wd, ws, bins=[0,2,4,6,8,10], normed=True, opening=0.8)

# 提取关键数据
directions = ax.directions  # 每个扇区的中心风向（如[0,22.5,45,...]）
sector_angle = ax._sector  # 每个扇区的角度宽度（默认22.5度，16个扇区）
freq = ax._info['freq']    # 各风速段在每个扇区的频率（shape: 风速段数 × 扇区数）
bins = ax.bins             # 风速区间（如[0,2,4,...10]）

2. 用simplekml将统计数据转换为KML扇形多边形 #

将每个风速段的扇区转换为KML多边形，需要把相对频率（半径比例）转换为实际地理坐标（经纬度），这里用UTM投影保证米到经纬度转换的准确性：

import simplekml
import pyproj

# 初始化KML对象
kml = simplekml.Kml()
# 风玫瑰中心的经纬度（替换为你的目标地点）
center_lon, center_lat = -74.0060, 40.7128

# 定义投影：WGS84经纬度 ↔ UTM平面坐标（根据中心纬度选UTM带，示例为18带）
proj_wgs84 = pyproj.Proj('EPSG:4326')
proj_utm = pyproj.Proj(proj='utm', zone=18, ellps='WGS84')

# 设置风玫瑰的最大实际半径（单位：米，可根据需求调整）
max_radius_m = 500
# 计算频率到实际米数的转换系数
max_freq = freq.max()
radius_per_freq = max_radius_m / max_freq

# 遍历每个风向扇区
for dir_idx, dir_center in enumerate(directions):
    # 计算扇区的起始/结束角度（从北顺时针）
    start_angle = dir_center - sector_angle / 2
    end_angle = dir_center + sector_angle / 2

    # 从外到内绘制每个风速段的扇形（避免遮挡）
    for bin_idx in reversed(range(len(bins)-1)):
        # 获取当前风速段的内外半径（米）
        r_inner = freq[bin_idx, dir_idx] * radius_per_freq if bin_idx > 0 else 0
        r_outer = freq[bin_idx+1, dir_idx] * radius_per_freq

        # 生成扇形的坐标点序列
        coords = [(center_lon, center_lat)]  # 先加入中心坐标
        # 绘制外圈点（从起始到结束角度）
        for angle in np.arange(start_angle, end_angle + 1, 2):  # 每2度取一个点，平衡精度和文件大小
            x_utm, y_utm = pyproj.transform(proj_wgs84, proj_utm, center_lon, center_lat)
            # 计算UTM坐标下的偏移（sin对应经度方向，cos对应纬度方向）
            x_utm += r_outer * np.sin(np.radians(angle))
            y_utm += r_outer * np.cos(np.radians(angle))
            lon, lat = pyproj.transform(proj_utm, proj_wgs84, x_utm, y_utm)
            coords.append((lon, lat))
        # 绘制内圈点（反向，闭合多边形）
        for angle in np.arange(end_angle, start_angle - 1, -2):
            x_utm, y_utm = pyproj.transform(proj_wgs84, proj_utm, center_lon, center_lat)
            x_utm += r_inner * np.sin(np.radians(angle))
            y_utm += r_inner * np.cos(np.radians(angle))
            lon, lat = pyproj.transform(proj_utm, proj_wgs84, x_utm, y_utm)
            coords.append((lon, lat))

        # 创建KML多边形
        pol = kml.newpolygon(name=f"Wind {bins[bin_idx]}-{bins[bin_idx+1]}km/h | Dir {dir_center}°")
        pol.outerboundaryis = coords

        # 提取windrose的颜色，转换为KML格式（AABBGGRR）
        face_color = ax.patches[dir_idx * len(bins) + bin_idx].get_facecolor()
        kml_color = simplekml.Color.rgb(
            int(face_color[0]*255), int(face_color[1]*255),
            int(face_color[2]*255), int(face_color[3]*255)
        )
        pol.style.polystyle.color = kml_color
        pol.style.polystyle.outline = 1
        pol.style.linestyle.color = simplekml.Color.black
        pol.style.linestyle.width = 1

# 保存KML文件
kml.save("windrose_google_earth.kml")

3. 关键注意事项 #

• 风向一致性：确保你的风向数据是从北顺时针定义（Google Earth和windrose默认规则），如果数据是其他格式（如从东逆时针），需提前转换角度。
• 坐标转换精度：UTM投影比直接计算经纬度偏移更准确，尤其在高纬度地区，需根据中心坐标选择正确的UTM带。
• 样式匹配：直接复用windrose生成的matplotlib颜色，能保证KML风玫瑰和示例图的视觉风格一致。
• 文件大小优化：减少扇形的角度采样点（如每5度一个点），避免KML文件过大影响加载速度。

内容的提问来源于stack exchange，提问作者Lucho