Python实现多边形内最大轴对齐矩形以掩膜图像的方法求助

1. 依赖库安装 #

先安装所需的处理库：

pip install geopandas shapely rasterio

2. 读取Shapefile并提取目标多边形 #

用geopandas读取矢量文件，提取需要处理的多边形：

import geopandas as gpd
from shapely.geometry import Polygon, box

# 读取shapefile
gdf = gpd.read_file("你的矢量文件.shp")
# 获取第一个多边形（可根据需求调整索引或遍历所有多边形）
target_polygon = gdf.geometry.iloc[0]

3. 寻找多边形内的最大轴对齐矩形 #

以下是直观的实现方式，通过枚举多边形顶点的坐标组合，筛选出完全位于多边形内部的面积最大轴对齐矩形：

def find_max_axis_aligned_rect(polygon):
    # 提取多边形所有顶点的x、y坐标
    coords = list(polygon.exterior.coords)
    xs = [x for x, y in coords]
    ys = [y for x, y in coords]
    
    # 排序去重，减少候选坐标数量
    xs_sorted = sorted(list(set(xs)))
    ys_sorted = sorted(list(set(ys)))
    
    max_area = 0
    max_rect = None
    
    # 枚举所有可能的x、y区间组合
    for i in range(len(xs_sorted)):
        x1 = xs_sorted[i]
        for j in range(i+1, len(xs_sorted)):
            x2 = xs_sorted[j]
            for k in range(len(ys_sorted)):
                y1 = ys_sorted[k]
                for l in range(k+1, len(ys_sorted)):
                    y2 = ys_sorted[l]
                    # 创建轴对齐矩形
                    rect = box(x1, y1, x2, y2)
                    # 判断矩形是否完全包含在多边形内
                    if polygon.contains(rect):
                        area = rect.area
                        if area > max_area:
                            max_area = area
                            max_rect = rect
    return max_rect

# 获取最大轴对齐矩形
max_rect = find_max_axis_aligned_rect(target_polygon)

注：如果多边形顶点数量极多，枚举法效率较低，可参考计算几何中的单调链分解、扫描线等优化算法提升性能。

4. 用矩形掩膜图像 #

处理带地理坐标系的栅格图像（如TIF） #

用rasterio直接基于地理坐标掩膜：

import rasterio
from rasterio.mask import mask

# 读取栅格图像
with rasterio.open("你的栅格图像.tif") as src:
    # 将矩形转换为rasterio要求的GeoJSON格式
    rect_geojson = [max_rect.__geo_interface__]
    # 执行掩膜并裁剪图像
    out_image, out_transform = mask(src, rect_geojson, crop=True)
    
    # 保存掩膜后的图像
    out_meta = src.meta.copy()
    out_meta.update({
        "height": out_image.shape[1],
        "width": out_image.shape[2],
        "transform": out_transform
    })
    
    with rasterio.open("掩膜结果.tif", "w", **out_meta) as dest:
        dest.write(out_image)

处理普通非地理图像（如JPG/PNG） #

需先将地理坐标转换为图像像素坐标，再用opencv生成掩膜：

import cv2
import numpy as np

# 已知图像的地理参数：左上角地理坐标(x0,y0)、分辨率(res_x, res_y)
x0, y0 = 100.0, 200.0  # 示例值，替换为实际参数
res_x, res_y = 0.5, -0.5  # res_y通常为负，因图像y轴向下

# 获取矩形的地理边界
minx, miny, maxx, maxy = max_rect.bounds

# 转换为像素坐标
x1 = int((minx - x0) / res_x)
y1 = int((y0 - maxy) / abs(res_y))
x2 = int((maxx - x0) / res_x)
y2 = int((y0 - miny) / abs(res_y))

# 读取图像并生成掩膜
img = cv2.imread("你的普通图像.jpg")
mask = np.zeros(img.shape[:2], dtype=np.uint8)
mask[y1:y2, x1:x2] = 255
masked_img = cv2.bitwise_and(img, img, mask=mask)

# 保存结果
cv2.imwrite("掩膜结果.jpg", masked_img)

注意事项 #

• 确保Shapefile与图像的坐标系一致，若不一致需先用gdf.to_crs()转换坐标
• 若多边形包含孔洞，polygon.contains(rect)会自动判断矩形是否位于外部多边形内且不在孔洞中

内容的提问来源于stack exchange，提问作者Vinícius D'Lucas