python填充多边形，获取所有内部点

shepely

简单示例

from shapely.geometry import Polygon

# 创建一个多边形
polygon = Polygon([(0, 0), (0, 5), (5, 5), (5, 0)])

# 提取多边形中的点
points = list(polygon.exterior.coords)

# 打印结果
print(points)

实际示例

import numpy as np
from shapely.geometry import Polygon, Point, MultiPolygon
import matplotlib.pyplot as plt
from matplotlib.path import Path
plt.rcParams['figure.facecolor'] = 'white'
plt.rcParams['axes.facecolor'] = 'silver'
# 网格线
plt.rcParams['xtick.direction'] = 'in'
print(plt.rcParams["font.sans-serif"])
# 设置全局字体为黑体
# plt.rcParams['font.sans-serif'] = ['simhei'] # 替换为其他中文字体如

def get_integer_points_in_polygon(polygon, include_boundary=True):
    """
    获取多边形内部的所有整数坐标点
    
    参数:
        polygon: Shapely Polygon对象（可以包含孔洞）
        include_boundary: 是否包含边界上的点
        
    返回:
        list: 多边形内部的整数坐标点列表 [(x1, y1), (x2, y2), ...]
    """
    # 获取多边形的边界框
    min_x, min_y, max_x, max_y = polygon.bounds
    
    # 将边界框坐标转换为整数
    min_x, min_y, max_x, max_y = map(int, [min_x, min_y, max_x, max_y])
    
    # 生成边界框内的所有整数点
    x_coords = np.arange(min_x, max_x + 1)
    y_coords = np.arange(min_y, max_y + 1)
    
    # 创建网格点
    xx, yy = np.meshgrid(x_coords, y_coords)
    grid_points = np.vstack((xx.flatten(), yy.flatten())).T
    
    # 检查每个点是否在多边形内部
    interior_points = []
    for point in grid_points:
        shapely_point = Point(point)
        if polygon.contains(shapely_point) or (include_boundary and polygon.touches(shapely_point)):
            interior_points.append(tuple(point))
    
    return interior_points


def visualize_polygon_and_points(polygon, points):
    """
    可视化多边形和内部点
    
    参数:
        polygon: Shapely Polygon对象
        points: 点列表
    """
    fig, ax = plt.subplots(figsize=(10, 8), facecolor='white')
    
    # 绘制多边形
    if hasattr(polygon, 'exterior'):
        # 单个多边形
        x, y = polygon.exterior.xy
        ax.plot(x, y, 'b-', linewidth=2)
        
        # 绘制孔洞
        for interior in polygon.interiors:
            x, y = interior.xy
            ax.plot(x, y, 'r-', linewidth=2)
    else:
        # 多多边形
        for poly in polygon.geoms:
            x, y = poly.exterior.xy
            ax.plot(x, y, 'b-', linewidth=2)
            
            for interior in poly.interiors:
                x, y = interior.xy
                ax.plot(x, y, 'r-', linewidth=2)
    
    # 绘制点
    if points:
        x_pts, y_pts = zip(*points)
        ax.plot(x_pts, y_pts, 'go', markersize=2)
    
    ax.set_aspect('equal')
    plt.title('多边形及其内部整数点')
    plt.xlabel('X坐标')
    plt.ylabel('Y坐标')
    plt.grid(True)
    plt.show()

# 使用示例
if __name__ == "__main__":
    # 创建一个带孔洞的多边形示例
    exterior = [(0, 0), (10, 0), (10, 10), (0, 10), (0, 0)]
    interior = [(2, 2), (8, 2), (8, 8), (2, 8), (2, 2)]
    
    polygon = Polygon(exterior, [interior])
    
    # 获取多边形内部的所有整数点
    points = get_integer_points_in_polygon(polygon)
    
    print(f"找到 {len(points)} 个内部整数点")
    print("前10个点:", points[:10])
    
    # 可视化
    visualize_polygon_and_points(polygon, points)
    

geopandas

处理非标准geojson

import json
import geopandas as gpd
import pandas as pd
from shapely.geometry import Polygon, GeometryCollection
import matplotlib.pyplot as plt

# 方法1：直接解析GeometryCollection
def parse_geometry_collection(geojson_data):
    """解析GeometryCollection格式的GeoJSON"""
    features = []
    
    if geojson_data.get('type') != 'GeometryCollection':
        raise ValueError("不是GeometryCollection格式")
    
    for geom in geojson_data.get('geometries', []):
        # 提取几何类型和坐标
        geom_type = geom.get('type')
        coords = geom.get('coordinates', [])
        
        # 提取属性
        properties = geom.get('properties', {})
        
        # 创建对应的shapely几何对象
        if geom_type == 'Polygon':
            geometry = Polygon(coords[0])  # 多边形外环
        elif geom_type == 'Point':
            from shapely.geometry import Point
            geometry = Point(coords)
        elif geom_type == 'LineString':
            from shapely.geometry import LineString
            geometry = LineString(coords)
        else:
            continue  # 跳过不支持的几何类型
        
        features.append({
            'type': 'Feature',
            'geometry': geometry,
            'properties': properties
        })
    print(features)
    return features

# 方法2：转换为标准GeoJSON格式
def convert_to_standard_geojson(geojson_data):
    """将GeometryCollection转换为标准FeatureCollection格式"""
    features = parse_geometry_collection(geojson_data)
    
    return {
        'type': 'FeatureCollection',
        'features': features
    }

# 方法3：创建GeoDataFrame
def create_geodataframe(geojson_data):
    """从GeometryCollection创建GeoDataFrame"""
    features = parse_geometry_collection(geojson_data)
    
    # 提取几何和属性
    geometries = [f['geometry'] for f in features]
    properties = [f['properties'] for f in features]
    
    # 创建GeoDataFrame
    gdf = gpd.GeoDataFrame(properties, geometry=geometries)
    
    return gdf

# 使用示例
if __name__ == "__main__":
    # 示例数据（根据您提供的格式）
    geojson_data = {
        "type": "GeometryCollection",
        "geometries": [
            {
                "type": "Polygon", 
                "properties": {"uID": 0, "label": "1"}, 
                "coordinates": [[[10436, 6198], [10481, 6198], [10504, 6221], [10436, 6198]]]
            },
            {
                "type": "Polygon", 
                "properties": {"uID": 1, "label": "2"}, 
                "coordinates": [[[10550, 6250], [10600, 6250], [10600, 6300], [10550, 6300], [10550, 6250]]]
            }
        ]
    }

    # 解析GeometryCollection
    features = parse_geometry_collection(geojson_data)
    print("解析出的要素数量:", len(features))
    
    # 转换为标准GeoJSON
    standard_geojson = convert_to_standard_geojson(geojson_data)
    print("标准GeoJSON格式已创建")
    
    # 创建GeoDataFrame
    gdf = create_geodataframe(geojson_data)
    print("GeoDataFrame信息:")
    print(gdf.head())
    
    # 可视化
    fig, ax = plt.subplots(figsize=(10, 8), facecolor='white')
    gdf.plot(ax=ax, color='lightblue', edgecolor='black')
    
    # 添加标签
    for idx, row in gdf.iterrows():
        centroid = row.geometry.centroid
        ax.text(centroid.x, centroid.y, row['label'], 
                fontsize=12, ha='center', va='center')
    
    plt.title("GeometryCollection中的多边形")
    plt.xlabel("X坐标")
    plt.ylabel("Y坐标")
    plt.grid(True)
    plt.show()
    
    # 保存为标准GeoJSON文件
    gdf.to_file("converted.geojson", driver="GeoJSON")
    print("已保存为标准GeoJSON文件: converted.geojson")

处理实际文件

import json

# 从文件读取GeometryCollection格式的GeoJSON
def read_geometry_collection_file(file_path):
    """从文件读取GeometryCollection格式的GeoJSON"""
    with open(file_path, 'r') as f:
        geojson_data = json.load(f)
    
    return geojson_data

# 使用示例
file_path = "your_file.geojson"  # 替换为实际文件路径
geojson_data = read_geometry_collection_file(file_path)

# 处理数据
gdf = create_geodataframe(geojson_data)
print(gdf.info())

确保多边形坐标闭合

def ensure_polygon_closed(coordinates):
    """确保多边形坐标闭合（首尾点相同）"""
    if not coordinates:
        return coordinates
    
    # 对于多边形，coordinates是一个列表的列表
    # 第一个列表是外环，后续列表是内环（孔）
    closed_coords = []
    for ring in coordinates:
        if len(ring) < 3:
            continue  # 跳过无效环（至少需要3个点）
        
        # 检查是否闭合（第一个点和最后一个点相同）
        first_point = ring[0]
        last_point = ring[-1]
        
        if first_point != last_point:
            # 添加第一个点作为最后一个点以闭合环
            closed_ring = ring + [first_point]
        else:
            closed_ring = ring
            
        closed_coords.append(closed_ring)
    
    return closed_coords

提取坐标信息

def extract_coordinates_from_geometry_collection(geojson_data):
    """从GeometryCollection中提取所有坐标"""
    coordinates = []
    
    if geojson_data.get('type') != 'GeometryCollection':
        return coordinates
    
    for geom in geojson_data.get('geometries', []):
        geom_type = geom.get('type')
        coords = geom.get('coordinates', [])
        properties = geom.get('properties', {})
        
        if geom_type == 'Polygon':
            # 多边形坐标 - 第一个列表是外环，后续列表是内环（孔）
            for ring in coords:
                coordinates.extend(ring)
        elif geom_type == 'Point':
            coordinates.append(coords)
        elif geom_type == 'LineString':
            coordinates.extend(coords)
        elif geom_type == 'MultiPolygon':
            for polygon in coords:
                for ring in polygon:
                    coordinates.extend(ring)
    
    return coordinates

# 使用示例
all_coords = extract_coordinates_from_geometry_collection(geojson_data)
print(f"提取到 {len(all_coords)} 个坐标点")
print("前10个坐标:", all_coords[:10])

生成mask

import numpy as np
from shapely.geometry import Polygon, MultiPolygon, Point
import matplotlib.pyplot as plt
from matplotlib.path import Path
import matplotlib.patches as patches
import json

def get_integer_points_from_geometry_collection(geojson_data):
    """
    从GeometryCollection中提取所有多边形内部的整数坐标点
    
    参数:
        geojson_data: GeometryCollection格式的GeoJSON数据
        
    返回:
        set: 所有多边形内部的整数坐标点集合
    """
    interior_points = set()
    
    if geojson_data.get('type') != 'GeometryCollection':
        return interior_points
    
    for geom in geojson_data.get('geometries', []):
        geom_type = geom.get('type')
        coords = geom.get('coordinates', [])
        properties = geom.get('properties', {})
        
        if geom_type == 'Polygon':
            # 创建Shapely多边形对象
            polygon = Polygon(coords[0], coords[1:])
            # 获取多边形内部的整数点
            points = get_integer_points_in_polygon(polygon)
            interior_points.update(points)
            
        elif geom_type == 'MultiPolygon':
            # 处理多多边形
            for polygon_coords in coords:
                polygon = Polygon(polygon_coords[0], polygon_coords[1:])
                points = get_integer_points_in_polygon(polygon)
                interior_points.update(points)
    
    return interior_points

def get_integer_points_in_polygon(polygon):
    """
    获取多边形内部的所有整数坐标点
    
    参数:
        polygon: Shapely Polygon对象
        
    返回:
        set: 多边形内部的整数坐标点集合
    """
    # 获取多边形的边界框
    min_x, min_y, max_x, max_y = polygon.bounds
    
    # 将边界框坐标转换为整数
    min_x, min_y, max_x, max_y = map(int, [np.floor(min_x), np.floor(min_y), 
                                          np.ceil(max_x), np.ceil(max_y)])
    
    # 生成边界框内的所有整数点
    x_coords = np.arange(min_x, max_x + 1)
    y_coords = np.arange(min_y, max_y + 1)
    
    # 创建网格点
    xx, yy = np.meshgrid(x_coords, y_coords)
    grid_points = np.vstack((xx.flatten(), yy.flatten())).T

    
    # 检查每个点是否在多边形内部
    include_boundary = True
    interior_points = []
    for point in grid_points:
        shapely_point = Point(point)
        if polygon.contains(shapely_point) or (include_boundary and polygon.touches(shapely_point)):
            interior_points.append(tuple(point))
    
    return interior_points

    # # 修正后的路径创建方法 - 确保代码数组是一维且长度匹配
    # exterior = np.array(polygon.exterior.coords)
    
    # # 创建路径代码 - 确保是一维数组
    # n_points = len(exterior)
    # codes = np.ones(n_points, dtype=Path.code_type) * Path.LINETO
    # codes[0] = Path.MOVETO
    
    # # 创建路径
    # path = Path(exterior, codes)
    
    # # 检查每个点是否在多边形内部
    # contains = path.contains_points(grid_points)
    # interior_points = set(tuple(map(int, point)) for point, inside in zip(grid_points, contains) if inside)
    
    # return interior_points

def create_mask_from_points(points, shape=None):
    """
    从点集合创建二值mask
    
    参数:
        points: 点集合
        shape: mask的形状 (height, width)。如果为None，则根据点坐标自动确定
        
    返回:
        numpy.array: 二值mask
    """
    if not points:
        return np.array([]), (0, 0)
    
    # 如果没有指定形状，则根据点坐标确定
    if shape is None:
        all_x = [p[0] for p in points]
        all_y = [p[1] for p in points]
        min_x, max_x = min(all_x), max(all_x)
        min_y, max_y = min(all_y), max(all_y)
        width = max_x - min_x + 1
        height = max_y - min_y + 1
        origin = (min_x, min_y)
    else:
        height, width = shape
        origin = (0, 0)
    
    # 创建空的mask
    mask = np.zeros((height, width), dtype=np.uint8)
    
    # 填充mask
    for point in points:
        x, y = point
        mask_y = y - origin[1]
        mask_x = x - origin[0]
        
        if 0 <= mask_y < height and 0 <= mask_x < width:
            mask[mask_y, mask_x] = 1
    
    return mask, origin

def visualize_polygon_and_mask(polygon, mask, origin, points=None):
    """
    可视化多边形和mask
    
    参数:
        polygon: Shapely Polygon对象
        mask: 二值mask
        origin: mask的原点坐标
        points: 点集合（可选）
    """
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))
    
    # 绘制多边形
    if hasattr(polygon, 'exterior'):
        # 单个多边形
        x, y = polygon.exterior.xy
        ax1.plot(x, y, 'b-', linewidth=2)
        
        # 绘制孔洞
        for interior in polygon.interiors:
            x, y = interior.xy
            ax1.plot(x, y, 'r-', linewidth=2)
    else:
        # 多多边形
        for poly in polygon.geoms:
            x, y = poly.exterior.xy
            ax1.plot(x, y, 'b-', linewidth=2)
            
            for interior in poly.interiors:
                x, y = interior.xy
                ax1.plot(x, y, 'r-', linewidth=2)
    
    # 绘制点
    if points:
        x_pts, y_pts = zip(*points)
        ax1.plot(x_pts, y_pts, 'go', markersize=2)
    
    ax1.set_aspect('equal')
    ax1.set_title('多边形及其内部整数点')
    ax1.set_xlabel('X坐标')
    ax1.set_ylabel('Y坐标')
    ax1.grid(True)
    
    # 绘制mask
    ax2.imshow(mask, cmap='gray', origin='lower')
    ax2.set_title('生成的Mask')
    ax2.set_xlabel(f'X (原点: {origin[0]})')
    ax2.set_ylabel(f'Y (原点: {origin[1]})')
    
    plt.tight_layout()
    plt.show()

def process_geojson_data(geojson_data):
    """
    处理GeoJSON数据并生成mask
    
    参数:
        geojson_data: GeoJSON数据
        
    返回:
        mask: 二值mask
        origin: mask的原点坐标
        interior_points: 内部点集合
    """
    # 提取多边形内部的所有整数点
    interior_points = get_integer_points_from_geometry_collection(geojson_data)
    print(f"找到 {len(interior_points)} 个内部整数点")
    
    # 创建mask
    mask, origin = create_mask_from_points(interior_points)
    print(f"Mask形状: {mask.shape}, 原点: {origin}")
    
    # 可视化
    # 注意：这里简化处理，只可视化第一个多边形
    if geojson_data.get('type') == 'GeometryCollection' and geojson_data.get('geometries'):
        first_geom = geojson_data['geometries'][0]
        if first_geom.get('type') == 'Polygon':
            polygon = Polygon(first_geom['coordinates'][0], first_geom['coordinates'][1:])
            visualize_polygon_and_mask(polygon, mask, origin, interior_points)
    
    return mask, origin, interior_points

# 使用示例
if __name__ == "__main__":
    # 示例GeometryCollection数据
    geojson_data = {
        "type": "GeometryCollection",
        "geometries": [
            {
                "type": "Polygon", 
                "properties": {"uID": 0, "label": "1"}, 
                "coordinates": [
                    [[10, 10], [50, 10], [50, 50], [10, 50], [10, 10]],  # 外环
                    [[20, 20], [40, 20], [40, 40], [20, 40], [20, 20]]   # 内环（孔洞）
                ]
            },
            {
                "type": "Polygon", 
                "properties": {"uID": 1, "label": "2"}, 
                "coordinates": [
                    [[60, 60], [80, 60], [80, 80], [60, 80], [60, 60]]
                ]
            }
        ]
    }
    
    # 处理数据
    mask, origin, interior_points = process_geojson_data(geojson_data)
    
    # 保存mask为图像
    plt.imsave('mask.png', mask, cmap='gray')
    print("Mask已保存为 mask.png")
    
    # 保存点坐标到文件
    with open('interior_points.txt', 'w') as f:
        for point in sorted(interior_points):
            f.write(f"{point[0]}, {point[1]}\n")
    print("点坐标已保存为 interior_points.txt")