python 地图上面标记站点外加画圆圈

第一个excel是一堆站点经纬度数据

第二个excel是两个雷达站点经纬度数据

第三个是地图shp数据

想画出雷达探测范围230km，和经纬度站点，放在地图里

#!/usr/bin/env python
# -*- coding:utf-8 -*-
"""
@author: Suyue
@file: lianxi.py
@time: 2025/05/22
@desc:研究区概况
"""
import geopandas as gpd
import pandas as pd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
import matplotlib.patches as mpatches

# 设置中文字体，解决乱码问题
plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用黑体
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题

# 设置投影 - 使用WGS1984地理坐标系
crs_wgs84 = ccrs.PlateCarree()

# 1. 读取盟市界SHP数据
shp_path = "G:/标准地图/内蒙古shp干旱业务机/盟市界.shp"  # 替换为您的实际路径
gdf = gpd.read_file(shp_path)

# 确保SHP数据在WGS84下
if gdf.crs != 'EPSG:4326':
    gdf = gdf.to_crs('EPSG:4326')

# 2. 读取站点数据
stations_path = "G:/装备科/论文/03-雨滴谱/雨滴谱数据/站号.xlsx"  # 替换为您的实际路径
stations_df = pd.read_excel(stations_path)
# 假设列名为'经度'和'纬度'
stations_gdf = gpd.GeoDataFrame(
    stations_df,
    geometry=gpd.points_from_xy(stations_df['经度'], stations_df['纬度']),
    crs="EPSG:4326"  # WGS84
)

# 3. 读取两个特殊站点的数据
special_stations_path = "G:/装备科/论文/03-雨滴谱/雨滴谱数据/新一代天气雷达.xlsx"  # 替换为您的实际路径
special_df = pd.read_excel(special_stations_path)
special_gdf = gpd.GeoDataFrame(
    special_df,
    geometry=gpd.points_from_xy(special_df['经度'], special_df['纬度']),
    crs="EPSG:4326"  # WGS84
)

# 4. 创建缓冲区 (230km半径)
# 先将数据转换为等距投影(如UTM)以创建精确的圆形缓冲区
utm_crs = special_gdf.estimate_utm_crs()
special_gdf_utm = special_gdf.to_crs(utm_crs)

buffers = []
for geom in special_gdf_utm.geometry:
    buffer = geom.buffer(230000)  # 230km = 230000米
    buffers.append(buffer)

# 将缓冲区转换回WGS84
buffers_gdf = gpd.GeoDataFrame(geometry=buffers, crs=utm_crs).to_crs("EPSG:4326")

# 5. 绘制地图
fig, ax = plt.subplots(figsize=(10, 10), subplot_kw={'projection': crs_wgs84})

# 添加盟市界
ax.add_geometries(gdf.geometry, crs=crs_wgs84, edgecolor='black', facecolor='none', linewidth=1)

# 计算合适的显示范围
all_points = pd.concat([
    stations_gdf.geometry.bounds,
    special_gdf.geometry.bounds,
    buffers_gdf.geometry.bounds
])
minx, miny = all_points[['minx', 'miny']].min()
maxx, maxy = all_points[['maxx', 'maxy']].max()

# 添加更大的边距（例如15%或20%）
margin = 0.2 * max(maxy - miny, maxx - minx)  # 增加到20%
ax.set_extent([
    minx - margin,
    maxx + margin,
    miny - margin,
    maxy + margin
], crs=crs_wgs84)

# 设置地图为正方形
ax.set_aspect('equal')

# 绘制普通站点
stations_gdf.plot(ax=ax, color='black', markersize=50, label='Station', transform=crs_wgs84)

# 绘制特殊站点
special_gdf.plot(ax=ax, color='blue', markersize=80, marker='*', label='Radar', transform=crs_wgs84)

# 绘制缓冲区圆
buffers_gdf.boundary.plot(ax=ax, color='blue', linewidth=2, alpha=0.5, transform=crs_wgs84)

# 添加指北针
def add_north(ax, labelsize=18, loc_x=0.95, loc_y=0.95, width=0.06, height=0.09, pad=0.14):
    """
    画一个比例尺带'N'文字注释
    主要参数如下
    :param ax: 要画的坐标区域 Axes实例 plt.gca()获取即可
    :param labelsize: 显示'N'文字的大小
    :param loc_x: 以文字下部为中心的占整个ax横向比例
    :param loc_y: 以文字下部为中心的占整个ax纵向比例
    :param width: 指南针占ax比例宽度
    :param height: 指南针占ax比例高度
    :param pad: 文字符号占ax比例间隙
    :return: None
    """
    minx, maxx = ax.get_xlim()
    miny, maxy = ax.get_ylim()
    ylen = maxy - miny
    xlen = maxx - minx
    left = [minx + xlen*(loc_x - width*.5), miny + ylen*(loc_y - pad)]
    right = [minx + xlen*(loc_x + width*.5), miny + ylen*(loc_y - pad)]
    top = [minx + xlen*loc_x, miny + ylen*(loc_y - pad + height)]
    center = [minx + xlen*loc_x, left[1] + (top[1] - left[1])*.4]
    triangle = mpatches.Polygon([left, top, right, center], color='k')
    ax.add_patch(triangle)
    ax.text(s='N',
            x=minx + xlen*loc_x,
            y=miny + ylen*(loc_y - pad + height),
            fontsize=labelsize,
            horizontalalignment='center',
            verticalalignment='bottom')

# 在右上角添加指北针
add_north(ax, labelsize=16, loc_x=0.95, loc_y=0.95, width=0.05, height=0.08, pad=0.1)

# 添加图例并放置在右下角
ax.legend(loc='lower right')

# 添加标题
plt.title("Station Distribution and Radar Detection Range Map")

# 添加经纬度网格
gl = ax.gridlines(draw_labels=True, linestyle='--', alpha=0.5)
gl.top_labels = False
gl.right_labels = False

# 手动添加一条在最右边的经线
# 获取当前的最大经度值，并向上取整到最近的整数
max_longitude = np.ceil(maxx)
gl.xlocator = plt.FixedLocator(np.arange(round(minx-1), round(maxx+2), 1))  # 调整xlocator以包含最右边的经线

plt.tight_layout()
plt.show()