请教：如何用Python根据已知的经纬度坐标在地图上给定半径画圆

import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import Point
from shapely.ops import transform
import pyproj

def draw_circle_on_map(center_lat, center_lon, radius_km, num_points=100):
    """
    根据经纬度中心点和半径(公里)在地图上画圆
    
    参数:
    center_lat: 中心点纬度
    center_lon: 中心点经度
    radius_km: 圆半径(公里)
    num_points: 圆的点数(默认100)
    """
    # 创建中心点
    center_point = Point(center_lon, center_lat)
    
    # 定义投影转换：WGS84经纬度 -> UTM投影(用于计算距离)
    wgs84 = pyproj.CRS('EPSG:4326')
    # 自动选择UTM区域
    utm_zone = int((center_lon + 180) / 6) + 1
    utm_crs = pyproj.CRS(f'EPSG:326{utm_zone:02d}' if center_lat >= 0 else f'EPSG:327{utm_zone:02d}')
    
    # 创建投影转换器
    transformer = pyproj.Transformer.from_crs(wgs84, utm_crs, always_xy=True)
    
    # 将中心点转换到UTM坐标
    center_utm = transform(transformer.transform, center_point)
    
    # 在UTM坐标系中创建圆
    angles = np.linspace(0, 2 * np.pi, num_points)
    circle_points_utm = []
    
    for angle in angles:
        # 计算圆上的点(UTM坐标)
        x = center_utm.x + radius_km * 1000 * np.cos(angle)  # 公里转米
        y = center_utm.y + radius_km * 1000 * np.sin(angle)
        circle_points_utm.append(Point(x, y))
    
    # 将圆点转换回经纬度
    transformer_back = pyproj.Transformer.from_crs(utm_crs, wgs84, always_xy=True)
    circle_points_wgs84 = [transform(transformer_back.transform, pt) for pt in circle_points_utm]
    
    # 提取经纬度坐标
    lons = [pt.x for pt in circle_points_wgs84]
    lats = [pt.y for pt in circle_points_wgs84]
    
    return lons, lats

# 使用示例
if __name__ == "__main__":
    # 示例：以北京天安门为中心，画一个50公里半径的圆
    center_lat = 39.9042  # 纬度
    center_lon = 116.4074  # 经度
    radius_km = 50
    
    # 获取圆的经纬度坐标
    circle_lons, circle_lats = draw_circle_on_map(center_lat, center_lon, radius_km)
    
    # 可视化
    plt.figure(figsize=(10, 8))
    plt.plot(circle_lons, circle_lats, 'b-', linewidth=2, label=f'{radius_km}km半径')
    plt.plot(center_lon, center_lat, 'ro', markersize=10, label='中心点')
    plt.xlabel('经度')
    plt.ylabel('纬度')
    plt.title(f'以({center_lat}, {center_lon})为中心，{radius_km}公里半径的圆')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.axis('equal')
    plt.show()
    
    # 打印前5个点的坐标
    print("圆的前5个点坐标：")
    for i in range(5):
        print(f"点{i+1}: 经度={circle_lons[i]:.6f}, 纬度={circle_lats[i]:.6f}")

核心思路：

1. 使用pyproj进行坐标投影转换，因为经纬度是球面坐标，不能直接使用平面几何
2. 先将中心点转换到UTM投影坐标系（平面直角坐标系）
3. 在UTM坐标系中用三角函数计算圆上的点
4. 将圆点坐标转换回WGS84经纬度坐标系

需要的库：

pip install numpy matplotlib shapely pyproj

关键点：

• 必须进行投影转换，直接在经纬度上画圆会变形
• UTM投影适合小范围区域的距离计算
• 使用shapely处理几何对象更方便

一句话建议：用投影转换解决球面坐标画圆问题。

圆上的点到圆心距离是相同的，你每隔一个角度计算出另一端的坐标就可以了，一般用等角线算法，距离很长的话可以用大圆算法