%pip install cartopy geopy -i https://mirrors.aliyun.com/pypi/simple/

画圆形

import matplotlib.pyplot as plt
# 设置中文字体（Windows系统）
plt.rcParams['font.sans-serif'] = ['SimHei']  # 黑体
# 解决负号显示问题
plt.rcParams['axes.unicode_minus'] = False  

fig, ax = plt.subplots(figsize=(6, 6))

# 绘制圆形
circle = plt.Circle((0, 0), radius=1000, color='blue', fill=False, linewidth=1)
ax.add_patch(circle)

# 设置图形属性
ax.set_xlim(-1100, 1100)
ax.set_ylim(-1100, 1100)
ax.set_aspect('equal')  # 确保圆不变成椭圆
ax.set_title('圆形')
plt.grid(True)
plt.show()

png

偏移点计算

import math

def calculate_offset(lat_lon, distance_m, bearing_deg):
    """
    根据方位角和距离计算偏移点
    参数:
        lat, lon: 起点坐标(度)
        distance_m: 距离(米)
        bearing_deg: 方位角(度，正北为0，顺时针增加)
    返回:
        (new_lat, new_lon): 新坐标
    """
    # 地球半径(米)
    earth_radius = 6371000
    
    # 转换为弧度
    lat_rad = math.radians(lat_lon[0])
    lon_rad = math.radians(lat_lon[1])
    bearing_rad = math.radians(bearing_deg)
    
    # 计算新纬度
    new_lat_rad = math.asin(math.sin(lat_rad) * math.cos(distance_m/earth_radius) + 
                  math.cos(lat_rad) * math.sin(distance_m/earth_radius) * math.cos(bearing_rad))
    
    # 计算新经度
    new_lon_rad = lon_rad + math.atan2(math.sin(bearing_rad) * math.sin(distance_m/earth_radius) * math.cos(lat_rad),
                              math.cos(distance_m/earth_radius) - math.sin(lat_rad) * math.sin(new_lat_rad))
    
    # 转换回度数
    new_lat = math.degrees(new_lat_rad)
    new_lon = math.degrees(new_lon_rad)
    
    return (new_lat, new_lon)

# 示例：从北京向东北方向(45度)移动1000米
new_lat, new_lon = calculate_offset((39.9042, 116.4074), 1000, 45)
print(f"新坐标: {new_lat:.6f}, {new_lon:.6f}")

搜索点计算

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from geopy.distance import geodesic

def find_search_points(center, maxR, r):
    colors = ['blue', 'green', 'black', 'pink', 'orange', 'gray', 'purple', 'brown', 'cyan', 'magenta', 'white', 'yellow', 'red']

    center_lon, center_lat = center
    iterTimes = math.ceil(maxR/r) # 迭代次数

    print(f'{center=},{maxR=},search r={r},{iterTimes=}')

    search_points = []

    # 创建地图
    fig = plt.figure(figsize=(10, 10))
    ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

    # 绘制经纬度线
    ax.gridlines(crs=ccrs.PlateCarree(), linewidth=0.5, linestyle='-', color='gray', alpha=0.5, draw_labels=True)

    def plot_circle(lon, lat, r, color='blue', alpha=0.5,linestyle='--'):
        width_km = 2*r/1000/(111*math.cos(math.radians(lat)))
        height_km = 2*r/1000/111
        rotation_angle = math.degrees(math.atan(width_km/height_km))
        # print(f"{rotation_angle=}")
        ax.add_patch(patches.Ellipse((lon,lat),width_km,height_km
                                     ,color=color,linestyle=linestyle,fill=False,angle=rotation_angle,alpha=alpha,transform=ccrs.PlateCarree()))

    # 绘制搜索范围
    plot_circle(center_lon,center_lat,maxR)

    # 绘制迭代搜索范围
    for t in range(iterTimes):
        if t==0:
            # 绘制中心点
            ax.scatter(center_lon,center_lat, color='red', s=20, transform=ccrs.PlateCarree())
            # 绘制迭代圆标签
            ax.text(center_lon, center_lat, "(0,0)", fontsize=12, transform=ccrs.PlateCarree())
            plot_circle(center_lon,center_lat,r)

            search_points.append((t, 1, "(0,0)", center_lon, center_lat, maxR, r, 0))
            continue

        points_num = 4 * t
        for n in range(points_num):
            p_lat,p_lon = calculate_offset((center_lat,center_lon), r*t, 360/points_num*n)
            p_name = f"({t},{360/points_num*n})"

            distance = geodesic((center_lat,center_lon), (p_lat, p_lon)).meters

            # 绘制迭代圆心
            ax.scatter(p_lon, p_lat, color=colors[t], s=20, transform=ccrs.PlateCarree())
            # 绘制迭代圆标签
            ax.text(p_lon, p_lat, p_name, fontsize=12, transform=ccrs.PlateCarree())
            # 绘制迭代搜索范围
            plot_circle(p_lon, p_lat, r, color=colors[t])

            search_points.append((t, n, p_name, p_lon, p_lat, maxR, r, distance))
    print(f"search points = {len(search_points)}")

    plt.title('GCJ-02')
    plt.grid(True)
    plt.show()

    return [(p[3],p[4]) for p in search_points], search_points

maxR = 1000 # 搜索范围
r=350 # 搜索步长
center = (116.4074, 39.9042) # 中心点
geos, search_results = find_search_points(center, maxR, r)
for i in geos:
    print(i)

png

(116.4074, 39.9042)
(116.40740000000001, 39.90734762562073)
(116.41150318241657, 39.90419992769774)
(116.40740000000001, 39.90105237437928)
(116.40329681758345, 39.90419992769774)
(116.40740000000001, 39.910495251241436)
(116.4132031532617, 39.908651270223686)
(116.4156063648158, 39.90419971079095)
(116.4132023992473, 39.89974844056726)
(116.40740000000001, 39.89790474875857)
(116.40159760075272, 39.89974844056726)
(116.39919363518422, 39.90419971079095)
(116.40159684673833, 39.908651270223686)