给定基站，求出每个基站的覆盖范围边界坐标序列

任务目标：

• 利用提供的NumPy数组中的经纬度数据，确定通信基站的信号覆盖范围。
• 确保每个基站的覆盖区域是不重叠的。
• 输出每个基站的覆盖范围边界的经纬序列。
• 可视化城市边界、基站点位置以及每个基站的覆盖范围边界。

数据描述：

• bound_array：一个NumPy数组，包含城市边界数据，城市边界由一系列经纬度坐标点组成，形成一个不规则多边形区域。
• point_array：一个NumPy数组，包含通信基站位置数据，每个基站由一个经纬度坐标点表示。

处理步骤：

1. 从bound_array中提取城市边界的经纬度坐标点序列。
2. 从point_array中提取通信基站的经纬度坐标点。
3. 根据基站坐标点，计算城市中每个点到最近基站的距离。
4. 根据距离确定每个点所属的基站覆盖范围。
5. 为每个基站定义其覆盖范围的边界，即不重叠的多边形区域。
6. 输出每个基站覆盖范围的边界经纬序列。
7. 使用Python代码进行可视化，绘制城市边界、基站点位置和每个基站的覆盖范围边界。

结果输出：

• 输出文件包含每个基站覆盖范围的边界经纬序列。
• 结果为dict，对应为每个基站所覆盖的范围的边界
• 可视化结果：一张图，展示城市边界、基站点位置和每个基站的覆盖范围边界

  import numpy as np
  from scipy.spatial import cKDTree
  import matplotlib.pyplot as plt
  from matplotlib.path import Path
  
  # 定义不规则多边形边界的顶点
  irregular_polygon = np.array([[1, 3], [4, 2], [5, 5], [2, 4], [1, 3]])
  # irregular_polygon = np.load('bound.npy')
  
  # 定义栅格的分辨率（网格大小）
  grid_resolution = 0.1
  
  # 计算栅格的坐标点范围
  min_x, max_x = np.min(irregular_polygon[:, 0]), np.max(irregular_polygon[:, 0])
  min_y, max_y = np.min(irregular_polygon[:, 1]), np.max(irregular_polygon[:, 1])
  
  # 创建栅格的坐标点
  x_values = np.arange(min_x, max_x + grid_resolution, grid_resolution)
  y_values = np.arange(min_y, max_y + grid_resolution, grid_resolution)
  grid_x, grid_y = np.meshgrid(x_values, y_values)
  
  # 初始化栅格分配
  grid_assignment = np.zeros((grid_x.shape[0], grid_y.shape[1]), dtype=int)
  
  # 定义中心点
  center_points = np.array([[3, 3], [4, 4]])
  # center_points = np.load('point.npy')
  
  # 计算每个栅格点到所有中心点的距离
  tree = cKDTree(center_points)
  distances, indices = tree.query(np.column_stack((grid_x.ravel(), grid_y.ravel())))
  
  # 分配每个栅格点最近的中心点
  grid_assignment.flat = np.arange(len(center_points))[indices]
  
  # 确定哪些栅格点在多边形内部
  polygon_path = Path(irregular_polygon)
  inside_mask = polygon_path.contains_points(np.column_stack((grid_x.ravel(), grid_y.ravel())))
  
  
  # 可视化结果，只显示在多边形内的栅格点
  plt.figure(figsize=(10, 10))
  sc = plt.scatter(grid_x.ravel()[inside_mask], grid_y.ravel()[inside_mask], c=grid_assignment.ravel()[inside_mask], cmap='viridis', alpha=0.5)
  plt.colorbar(sc)
  plt.plot(irregular_polygon[:, 0], irregular_polygon[:, 1], 'k-', linewidth=2, label='Polygon Boundary')
  plt.scatter(center_points[:, 0], center_points[:, 1], c='red', marker='o', s=100, label='Center Points')
  plt.title('Grid Assignment to Center Points (Inside Polygon Only)')
  plt.xlabel('X Coordinate')
  plt.ylabel('Y Coordinate')
  plt.legend()
  
  bound_list = []
  inside_mask = inside_mask.reshape(grid_assignment.shape)
  directions = [(-1, 0), (1, 0), (0, -1), (0, 1), (-1, -1), (-1, 1), (1, -1), (1, 1)]
  def is_adjacent(x, y):
      # 检查坐标是否在数组范围内
      if x < 0 or x >= grid_assignment.shape[0]-1 or y < 0 or y >= grid_assignment.shape[1]-1:
          return False
      # 检查相邻像素
      return True
  id = 0
  for x in range(grid_assignment.shape[0]):
      for y in range(grid_assignment.shape[1]):
          id = id+1
          for dx, dy in directions:
  #             print(x,y)
              if is_adjacent(x,y):
                  if ((inside_mask[x][y]==True and grid_assignment[x + dx][y + dy]!=grid_assignment[x][y]) or (inside_mask[x][y]==True and inside_mask[x + dx][y + dy]==False)):
                      bound_list.append((x + dx,y + dy))
  #                     print(grid_assignment[x + dx][y + dy],grid_assignment[x][y])
  #                     print(inside_mask[x + dx][y + dy])
  #                     print(x,y)
                      break
                  
  bound_list = list(set(bound_list))
  # plt.figure(figsize=(10, 10))
  # # plt.scatter(bound_list[:, 0], bound_list[:, 1], c='green', marker='o', s=100)
  # plt.scatter([x for x, y in bound_list], [y for x, y in bound_list],c='green')
  
  # plt.scatter(grid_x.ravel()[inside_mask], grid_y.ravel()[inside_mask], c=grid_assignment.ravel()[inside_mask], cmap='viridis', alpha=0.5)
  
  ans_bound = []
  
  for x,y in bound_list:
  #     print(x,y)
  #     print(grid_x[0][y],grid_y[x][0])
      ans_bound.append((grid_x[0][y],grid_y[x][0]))
  plt.scatter([x for x, y in ans_bound], [y for x, y in ans_bound],c='green',s=5)
  
  
  
  plt.show()
      
  import numpy as np
  from scipy.spatial import cKDTree
  
  # 假设我们有以下点和中心点的坐标列表
  points = ans_bound
  
  # 使用cKDTree为中心点创建一个KD树，以便快速查询最近邻
  tree = cKDTree(center_points)
  
  # 计算每个点到所有中心点的距离，并找到最近的中心点
  distances, indices = tree.query(points)
  
  # 创建一个字典来存储分配结果
  assignment = {}
  
  # 遍历每个点和它的最近中心点的索引
  for point, index in zip(points, indices):
      # 使用元组表示点的坐标，以便可以作为字典的键
      point_key = tuple(point)
      # 使用元组表示中心点的坐标
      center_point = tuple(center_points[index])
      
      # 如果中心点还没有在字典中，初始化一个空列表
      if center_point not in assignment:
          assignment[center_point] = []
      # 将点添加到对应中心点的列表中
      assignment[center_point].append(point_key)
  
  # 输出结果
  print(assignment)
  
  import pickle
  with open('ans_dict.pickle', 'wb') as file:
      pickle.dump(assignment, file)