如何用Python为MapInfo创建基站及扇形饼图？代码问题求助

核心问题点 #

• 角度单位不匹配：方位角（AZIMUTH）和波束宽度（BW）一般是角度制，但代码中使用cos/sin计算坐标时，仅转换了步长角度，未将起始/结束角度转为弧度，导致坐标计算完全错误，这是图形显示异常的主要原因。
• 数值列被强制转字符串：代码把EASTING、NORTHING、AZIMUTH等需要参与数学运算的数值列转成了字符串，直接引发运算报错，是部分场景运行失败的关键原因。
• 无效冗余代码：最后修改radius *=3但未实际使用新半径，重复设置angle_step，属于无效操作。
• 多边形构建精度问题：用np.arange生成角度序列可能漏掉最后一个端点，导致多边形不闭合；buffer和simplify的参数设置过小，对米制坐标系（EPSG:27700）几乎无意义，还可能破坏图形。

具体修复步骤 #

1. 统一角度单位：将起始/结束角度全部转为弧度后再参与坐标计算。
2. 调整数据类型转换逻辑：仅将非数值型列转为字符串，保留需要计算的列的数值类型。
3. 移除无效代码：删掉未用到的半径修改、重复的角度步长设置。
4. 优化多边形构建：用np.linspace替代np.arange生成角度序列，确保包含首尾端点；调整buffer/simplify参数或按需移除，提前修复无效多边形。

修正后的完整代码 #

import numpy as np
from shapely.geometry import Polygon
import geopandas as gpd
import os

def create_sector_pie(row):
    try:
        # 强制转换为浮点型，避免字符串运算报错
        azimuth = float(row['AZIMUTH'])
        beamwidth = float(row['BW'])
        radius = float(row['Radius'])
        easting = float(row['EASTING'])
        northing = float(row['NORTHING'])

        # 计算扇形起始/结束角度，统一转为弧度制
        start_angle_rad = np.radians(azimuth - beamwidth / 2)
        end_angle_rad = np.radians(azimuth + beamwidth / 2)

        # 用linspace生成角度序列，确保包含首尾端点，避免多边形不闭合
        angle_step = np.radians(0.5)
        num_steps = int(np.ceil((end_angle_rad - start_angle_rad) / angle_step)) + 1
        angles = np.linspace(start_angle_rad, end_angle_rad, num_steps)

        # 生成扇形坐标点
        center = (easting, northing)
        polygon_coords = [center]
        for angle in angles:
            x = center[0] + radius * np.cos(angle)
            y = center[1] + radius * np.sin(angle)
            polygon_coords.append((x, y))
        polygon_coords.append(center)  # 闭合多边形

        # 创建多边形并优先修复有效性
        polygon = Polygon(polygon_coords)
        if not polygon.is_valid:
            polygon = polygon.buffer(0)

        # 按需保留buffer和simplify（EPSG:27700为米制，参数建议用米级数值）
        # polygon = polygon.buffer(0.1)
        # polygon = polygon.simplify(tolerance=0.1)

        return polygon
    except Exception as e:
        print(f"错误信息: {e}")
        print(f"出错的行数据:\n{row}")
        return None

# 生成几何列
data['geometry'] = data.apply(create_sector_pie, axis=1)

# 过滤无效几何行
data = data.dropna(subset=['geometry'])

# 仅转换非数值型列为字符串
columns_to_convert = ['CSR', 'CELL_ID', 'SECTOR', 'MI_MYCOM_ID', 'GENERATION', 'CELL_TYPE', 
                      'FDD', 'CELL_3G_CARRIER', 'AREA', 'JV_ID', 'PCI', 'TAC', 'Band']
data[columns_to_convert] = data[columns_to_convert].astype(str)

# 创建GeoDataFrame并设置坐标系
gdf = gpd.GeoDataFrame(data, geometry='geometry')
gdf.crs = 'EPSG:27700'

# 保存为MapInfo TAB文件
tab_file = os.path.join(output_folder, 'N5_sectors.tab')
gdf.to_file(tab_file, driver='MapInfo File')

额外说明 #

• 如果你的方位角是顺时针从北开始定义（GIS常见规则），而数学计算默认是逆时针从X轴开始，需要调整角度转换逻辑，比如start_angle_rad = np.radians(90 - azimuth - beamwidth / 2)，具体需匹配你的数据定义。
• EPSG:27700为米制坐标系，radius单位是米，buffer/simplify的参数建议设置为0.1米以上才有实际效果。

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