Skip to main content

shp文件处理

shp文件基础​

Shape files 是 ESRI 提供的一种矢量数据格式，它没有拓扑信息，一个 Shape files 由一组文件组成，其中必要的基本文件包括坐标文件（ .shp ）、索引文件（ .shx ）和属性文件（ .dbf ）三个文件。

坐标文件 (.shp) 用于记录空间坐标信息。它由头文件和实体信息两部分构成，坐标文件的文件头是一个长度固定 (100 bytes) 的记录段，一共有 9 个 int 型和 7 个 double 型数据，

文件头​

起始位置名称数值类型位序
0File Code9994Integerbig
4Unused0Integerbig
8Unused0Integerbig
12Unused0Integerbig
16Unused0Integerbig
20Unused0Integerbig
24文件长度文件的实际长度Integerbig
28版本号1000IntegerLittle
32几何类型表示这个 Shapefile 文件所记录的空间数据的几何类型IntegerLittle
36Xmin空间数据所占空间范围的 X 方向最小值DoubleLittle
44Ymin空间数据所占空间范围的 Y 方向最小值DoubleLittle
52Xmax空间数据所占空间范围的 X 方向最大值DoubleLittle
60Ymax空间数据所占空间范围的 Y 方向最大值DoubleLittle
68*Zmin空间数据所占空间范围的 Z 方向最小值DoubleLittle
76*Zmax空间数据所占空间范围的 Z 方向最大值DoubleLittle
84*Mmin最小 Measure 值DoubleLittle
92*Mmax最大 Measure 值DoubleLittle

几何类型表​

编号几何类型
0Null Shape （表示这个 Shapefile 文件不含坐标）
1Point （表示 Shapefile 文件记录的是点状目标，但不是多点）
3PolyLine （表示 Shapefile 文件记录的是线状目标）
5Polygon （表示 Shapefile 文件记录的是面状目标）
8MultiPoint （表示 Shapefile 文件记录的是多点，即点集合）
11PointZ （表示 Shapefile 文件记录的是三维点状目标）
13PolyLineZ （表示 Shapefile 文件记录的是三维线状目标）
15PolygonZ （表示 Shapefile 文件记录的是三维面状目标）
18MultiPointZ （表示 Shapefile 文件记录的是三维点集合目标）
21PointM （表示含有 Measure 值的点状目标）
23PolyLineM （表示含有 Measure 值的线状目标）
25PolygonM （表示含有 Measure 值的面状目标）
28MultiPointM （表示含有 Measure 值的多点目标）
31MultiPatch （表示复合目标）

常用的坐标系​

{
  "global": [
    {
      "name": "WGS 84",
      "epsg": 4326,
      "description": "最常用的坐标系，用经度和纬度表示地球上的点。"
    },
    {
      "name": "Web Mercator",
      "epsg": 3857,
      "description": "用于 Web 地图投影的坐标系，主要用于 Google Maps、OpenStreetMap 等在线地图。"
    }
  ],
  "china": [
    {
      "name": "GCJ-02",
      "description": "中国国内常用的坐标系，用于在线地图和导航应用，但不适用于地理信息系统（GIS）数据。"
    },
    {
      "name": "BD-09",
      "description": "用于百度地图的坐标系，与标准 WGS 84 和 GCJ-02 坐标系有一定的差异。"
    },
    {
      "name": "Xian 1980",
      "epsg": 4610,
      "description": "中国西安80坐标系，用于国内测绘和GIS应用。"
    },
    {
      "name": "Beijing 1954",
      "epsg": 4214,
      "description": "中国北京54坐标系，用于早期的测绘和GIS应用。"
    }
  ],
  "cgcs2000": [
    {
      "name": "CGCS2000",
      "epsg": 4490,
      "description": "中国大地坐标系统2000，用于中国国内的测绘和GIS应用。"
    }
  ]
}

转换坐标系​

gdf = gdf.to_crs("EPSG:4326") # WGS 84

gdf = gdf.to_crs("EPSG:4490") # CGCS2000

常用模块​

import geopandas as gpd
import shapely

from shapely.geometry import Polygon, Point

读取一个shp文件​


# 检查是否存在 .cpg 后缀文件
cpg_file_path = shp_file_path.replace(".shp", ".cpg")
encoding = "utf-8"  # 默认编码为 utf-8

try:
    with open(cpg_file_path, "r") as cpg_file:
        encoding = cpg_file.read().strip()
except FileNotFoundError:
    print("打开文件错误")
    return []
    
# 读取面SHP文件
gdf = gpd.read_file(shp_file_path, encoding=encoding)

创建矩形​

依赖​

import geopandas as gpd
from shapely.geometry import Polygon
from shapely.affinity import rotate
from typing import Tuple, Optional

根据【中心点】​

def create_rectangle_by_point(center: Tuple[float, float], width: float, height: float, rotation: float = 0.0,output_shp_path: Optional[str] = None, crs: Optional[str] = None) -> Polygon:
    """
    创建一个矩形的 shp 文件，根据给定的中心点、宽度和高度以及可选的旋转角度和坐标系。
    
    Args:
        center (Tuple[float, float]): 矩形中心坐标 (x, y)。
        width (float): 矩形宽度。
        height (float): 矩形高度。
        rotation (float, optional): 矩形的旋转角度（默认为 0）。
        output_shp_path (str, optional): 输出的 shp 文件路径。如果为 None，则不输出文件（默认为 None）。
        crs (str, optional): 输出的 shp 文件的坐标系（默认为 None，使用默认坐标系）。

    Returns:
        Polygon: 创建的矩形多边形。
    """
    half_width = width / 2.0
    half_height = height / 2.0

    # 计算矩形的四个角点
    corners = [
        (center[0] - half_width, center[1] - half_height),
        (center[0] + half_width, center[1] - half_height),
        (center[0] + half_width, center[1] + half_height),
        (center[0] - half_width, center[1] + half_height)
    ]
    
    # 将入参导出为字段
    data = {
        "page_id":[1],
        "x": [center[0]],
        "y": [center[1]],
        "width": [width],
        "height": [height],
        "rotation": [rotation]
    }

    # 创建一个 GeoDataFrame 包含矩形多边形
    rectangle_polygon = Polygon(corners)
    gdf = gpd.GeoDataFrame(data, geometry=[rectangle_polygon])
    
    # 设置坐标系
    if crs:
        gdf.crs = crs
            
    # 如果有旋转角度，则进行矩形旋转
    if rotation != 0.0:
        rectangle_polygon = rotate(rectangle_polygon, rotation, origin=center)

    if output_shp_path:
        # 导出矩形 shp 文件
        gdf.to_file(output_shp_path, driver="ESRI Shapefile")

    return gdf


if __name__ == "__main__":
    center_coords = (10.0, 20.0)  # 矩形中心坐标 (x, y)
    rectangle_width = 8.0  # 矩形宽度
    rectangle_height = 6.0  # 矩形高度
    rectangle_rotation = 30.0  # 旋转角度（单位：度）
    output_rectangle_shp = "./test/rectangle.shp"  # 输出的矩形 shp 文件名

    # 设置输出 shp 文件的坐标系（可选）
    output_crs = "EPSG:4326"

    # 创建矩形并导出 shp 文件
    created_rectangle = create_rectangle_shp(center_coords, rectangle_width, rectangle_height, rectangle_rotation,
                                             output_rectangle_shp, output_crs)
    print("创建的矩形多边形：", created_rectangle)

根据shp外轮廓​

def create_rectangle(input_target: str | list, output_shp_path: Optional[str] = None):
    """
    @Description 根据一个shp轮廓或者指定的，生成一个项目合理的显示范围，生成点集最大轮廓的矩形并返回矩形的 SHP 文件。

    - param input_target         :{str}           shp文件绝对路径或者[startX,stattY,endX,endY]点
    - param output_shp_path=None :{Optional[str]} 输出的矩形 SHP 文件路径。如果为 None，则不输出文件（默认为 None）。

    @returns `{}` {description}

    """
    corners = []
    crs = None

    print("input_target: ", input_target)

    if isinstance(input_target, str):
        gdf = gpd.read_file(input_target)

        # 将点集合并成一个 MultiPoint
        multi_point = unary_union(gdf["geometry"])

        # 获取最大轮廓的外接矩形
        envelope = multi_point.envelope

        # 获取矩形的四个角点
        min_x, min_y, max_x, max_y = envelope.bounds

        corners = [(min_x, min_y), (max_x, min_y), (max_x, max_y), (min_x, max_y)]

        crs = gdf.crs

    elif isinstance(input_target, list) and len(input_target) == 4:
        min_x, min_y, max_x, max_y = input_target
        corners = [(min_x, min_y), (max_x, min_y), (max_x, max_y), (min_x, max_y)]

    else:
        raise "不支持的 input_target"

    rectangle = Polygon(corners)

    if output_shp_path:
        new_gdf = gpd.GeoDataFrame({"show_name": ["rect"]}, geometry=[rectangle])

        # 设置坐标系为输入点集的坐标系
        if crs:
            new_gdf.crs = crs

        # 导出矩形 SHP 文件
        new_gdf.to_file(output_shp_path, driver="ESRI Shapefile")

    return rectangle

生成圆形​

def create_circle_shp(
    center_point: Tuple[float, float],
    radius: float,
    output_shp_path: str,
    crs: str = None,
) -> None:
    """
    @Description {description}

    - param radius          :{float} 圆的半径。
    - param output_shp_path :{str}   文件路径。
    - param crs             :{str}   坐标系： "EPSG:4326"

    @returns `{ None }` {description}

    """

    if not isinstance(center_point, Point):
        center_point = Point(center_point)

    circle = center_point.buffer(radius)  # 创建缓冲区，即圆形
    circle_polygon = Polygon(circle.exterior.coords)

    # 创建一个 GeoDataFrame 包含圆形多边形
    gdf = gpd.GeoDataFrame(geometry=[circle_polygon])

    # 设置坐标系
    if crs:
        gdf.crs = crs

    # 导出圆形 shp 文件
    if output_shp_path:
        gdf.to_file(output_shp_path, driver="ESRI Shapefile")

    return circle_polygon

检查shp文件的dtype​

def shapefile_fields(shp_path, change_dtype=None, output_path=None):
    """
    读取指定的 SHP 文件并打印出字段信息，可选修改指定字段的数据类型，并可选保存修改后的 GeoDataFrame。

    Args:
        shp_path (str): SHP 文件的路径。
        change_dtype (dict, optional): 要修改的字段数据类型的字典，键是字段名称，值是目标数据类型字符串（例如："int64"）。
        output_path (str, optional): 保存修改后 GeoDataFrame 的路径。如果为 None，则不保存（默认为 None）。
    """
    gdf = gpd.read_file(shp_path)
    res = True

    if change_dtype:
        for field_name, new_dtype in change_dtype.items():
            if field_name in gdf.columns:
                if gdf[field_name].dtype != new_dtype:
                    print(f"注意，当前字段{field_name}已存在，但类型不是{new_dtype}: ")
                    res = False
            else:
                # 如果字段不存在，则根据提供的字段类型新建字段
                gdf[field_name] = gdf.apply(lambda row: None, axis=1)  # 创建空字段
                gdf[field_name] = gdf[field_name].astype(new_dtype)

    # 打印字段信息
    # 获取字段名称和数据类型信息
    # field_info = gdf.dtypes
    # print("字段信息：")
    # for field_name, data_type in field_info.iteritems():
    #     print(f"字段名：{field_name}，数据类型：{data_type}")

    if output_path:
        gdf.to_file(output_path, driver="ESRI Shapefile")

    return res

判断是点、线、面​

import geopandas as gpd

def check_geometry_type(shp_file_path: str) -> str:
    """
    判断给定SHP文件的几何类型。

    参数：
    - shp_file_path: SHP文件路径

    返回：
    - 几何类型的字符串表示，可能的值为 ：
    - 'Point'、 点
    - 'Polygon'、 面
    - 'LineString'、 线
    - 'MultiPoint'、 多个点
    - 'MultiPolygon'、多个面
    - 'MultiLineString' 多个线
    - 'Unknown'
    """
    # 使用 GeoPandas 直接读取文件
    gdf = gpd.read_file(shp_file_path)

    # 获取第一个几何对象
    first_geometry = gdf.geometry.iloc[0] if not gdf.empty else None

    # 判断几何类型
    if first_geometry is not None:
        geom_type = first_geometry.geom_type
        if geom_type == 'Polygon':
            return 'Polygon'
        elif geom_type == 'LineString':
            return 'LineString'
        elif geom_type == 'Point':
            return 'Point'
        elif geom_type == 'MultiPolygon':
            return 'MultiPolygon'
        elif geom_type == 'MultiLineString':
            return 'MultiLineString'
        elif geom_type == 'MultiPoint':
            return 'MultiPoint'
    return 'Unknown'

# 示例用法
shp_file_path = "/path/to/your/file.shp"
geometry_type = check_geometry_type(shp_file_path)
print("几何类型:", geometry_type)


挖空​

两个相交的部分进行挖空，要挖空的shp支持多个多边形

def polygon_clip_polygon(polygon_shp: str, clip_rang_shp: str, output_shp: str) -> str:
    """
    @Description 两个相交的部分进行挖空，要挖空的shp支持多个多边形

    - param polygon_shp   :{str} 要挖空的shp
    - param clip_rang_shp :{str} 要挖空的范围，支持多个多边形的shp
    - param output_shp    :{str} 要输出的shp路径
    """
    try:
        # 读取shp文件
        input_polygon = gpd.read_file(polygon_shp)
        subtract_polygon = gpd.read_file(clip_rang_shp)

        res_polygon = input_polygon
        for i, row in subtract_polygon.iterrows():
            each_polygon = row.geometry

            # 进行挖空操作
            res_polygon = res_polygon.difference(each_polygon)

        res_polygon.to_file(output_shp, encoding="utf-8")

        return output_shp
    except Exception as e:
        print("polygon_clip_polygon fail")

        return ""

获取shp文件的全部文件​

shp文件一般由多种不同格式的同名文件组成，返回一个shp文件的所有相关文件列表，所有文件经过path.exists()判断过必然合法存在才返回结果，否则返回空列表

def get_shp_file_list(shp_name:str) -> list[str]:
    """
    @Description shp文件一般由多种不同格式的同名文件组成，返回一个shp文件的所有相关文件列表

    - param shp_name :{param} shp文件的绝对路径，回判断真实存在才返回

    @returns `{ list[str] }` {description}

    """
    exts = [".cpg", ".dbf", ".sbn", ".sbx", ".shp", ".shx", ".shp.xml"]
    name, _ = path.splitext(shp_name)
    res = []
    for ext in exts:
        tar = "{}{}".format(name, ext).decode("utf-8")
        if path.exists(tar):
            res.append(tar)

    return res

获取shp文件的范围坐标​

import geopandas as gpd
from pydantic import BaseModel
from rasterio.transform import Affine

class Coordinates(BaseModel):
    left_top: tuple[float, float]
    right_top: tuple[float, float]
    left_bottom: tuple[float, float]
    right_bottom: tuple[float, float]
    
def get_shapefile_corners(filename: str) -> Coordinates:
    try:
        df = gpd.read_file(filename)

        # 获取第一个几何对象
        geometry = df.geometry[0]

        print(geometry.bounds)

        # 获取四个坐标点
        bounds = geometry.bounds

        return Coordinates(
            left_top=(bounds[0], bounds[3]),
            right_top=(bounds[2], bounds[3]),
            left_bottom=(bounds[0], bounds[1]),
            right_bottom=(bounds[2], bounds[1]),
        )
    except Exception as e:
        print(f"Error occurred: {str(e)}")
        return None

平移 shp​

方案1：​

支持shp文件或者直接点坐标，并带坐标系修正

def shp_translate_by_shp(
    move_shp: str,
    refs: str | tuple[float, float],
    output_shp: str,
    fix_crs: str | bool = None,
) -> str:
    """
    @Description 平移一个shp文件，或者将两个shp文件移动到同一视图中

    - param move_shp   :{str}                   要移动的shp
    - param refs       :{str}                   shp文件路径
    - param refs       :{tuple[float, float]}   支持tuple(x,y)要移动到的的坐标位置，以左下角为准
    - param output_shp :{str}                   导出的shp绝对路径
    - param fix_crs    :{str|bool}              坐标系，是否修复坐标系

    @returns `{ str}` 成功返回输出路径，失败返回空
    @example
    ```python
    input_shp1 = r"F:\要移动的.shp"
    input_shp2 = r"Z:\参考标注.shp"
    output_shp = path.join(path.dirname(input_shp2), "output.shp")

    shp_translate_by_shp(input_shp1, input_shp2, output_shp)
    # or
    shp_translate_by_shp(input_shp1, (x,y), output_shp)
    ```
    """
    """
    try:

        # 读取Shapefile文件
        move_gdf = gpd.read_file(move_shp)
        dx = dy = None

        if isinstance(refs, str) and refs.endswith(".shp"):
            refs_gdf = gpd.read_file(refs)
            dx = refs_gdf.total_bounds[0] - move_gdf.total_bounds[0]
            dy = refs_gdf.total_bounds[1] - move_gdf.total_bounds[1]

            # 检查坐标系
            if refs_gdf.crs != move_gdf.crs:
                if fix_crs is None:
                    print(f"【警告】，当前shp文件坐标不一致: {refs_gdf.crs}=>{move_gdf.crs}")
                elif isinstance(fix_crs, str):
                    move_gdf.crs = fix_crs
                elif isinstance(fix_crs, bool):
                    move_gdf.crs = refs_gdf.crs
                else:
                    print(f"【警告】，当前shp文件坐标不一致: {refs_gdf.crs}=>{move_gdf.crs}")

        elif isinstance(refs, tuple):
            dx = refs[0] - move_gdf.total_bounds[0]
            dy = refs[1] - move_gdf.total_bounds[1]
            if fix_crs:
                move_gdf.crs = fix_crs

        else:
            print("暂不支持的输入类型")
            return ""

        move_gdf.geometry = move_gdf.geometry.translate(dx, dy)
        move_gdf.to_file(output_shp)

        return output_shp

    except Exception as e:
        print(f"shp_translate_by_shp fail: {e}")
        return ""
        
input_shp1 = r"F:\要移动的.shp"
input_shp2 = r"Z:\参考标注.shp"
output_shp = path.join(path.dirname(input_shp2), "output.shp")

shp_translate_by_shp(input_shp1, input_shp2, output_shp)
# or
shp_translate_by_shp(input_shp1, (x,y), output_shp)

方案2：​

先获取点再平移

import geopandas as gpd
from pydantic import BaseModel
from rasterio.transform import Affine

def shift_and_align_shapefile(
    shp_file: str, output_shp: str, bottom_left: tuple[float, float]
) -> str:
    """
    用来将两个不在同一视图的shp文件，平移进同一视图，并输保存成新的shp

    - param shp_file    :{param} 要移动的shp
    - param output_shp :{param} 新shp的绝对路径
    - param bottom_left :{param} 要偏移的坐标位置，以左下角为准
    """
    try:
        # 读取Shapefile文件
        gdf = gpd.read_file(shp_file)

        # 计算Shapefile当前的边界框
        bounds = gdf.total_bounds

        # 计算Shapefile当前的宽度和高度
        width = bounds[2] - bounds[0]
        height = bounds[3] - bounds[1]

        # 计算目标边界框的左下角坐标
        new_bottom_left = bottom_left

        # 计算目标边界框的平移量
        dx = new_bottom_left[0] - bounds[0]
        dy = new_bottom_left[1] - bounds[1]

        # 创建仿射变换对象
        transform = Affine.translation(dx, dy)

        # 进行平移操作
        gdf.geometry = gdf.geometry.translate(dx, dy)

        # 将平移后的数据保存为新的Shapefile文件
        gdf.to_file(output_shp)

        print(f"Shapefile文件已保存为: {output_shp}")
        return output_shp

    except Exception as e:
        print(f"Error: {e}")
        return ""
# 将test.shp的范围移动到河道范围附近
test_file = r"./test.shp"

shp_coords = get_shapefile_corners(r"./河道范围.shp")

output_shp = "./output.shp"
shift_and_align_shapefile(test_file, output_shp, shp_coords.left_bottom)

shp_coords = get_shapefile_corners(output_shp)
print("output_shp: ", shp_coords)

mesh转换成shp​

from mikeio import Mesh

def mesh_to_shp(file_path: str, output_path: str) -> str:
    msh = Mesh(file_path)

    # 获取shp对象
    shp = msh.to_shapely()
    buffer = shp.buffer(0)

    # 导出shp文件
    gdf = gpd.GeoSeries([buffer])
    gdf.to_file(output_path)

裁剪shp​

优化版​

def clip_shp(clip_geometry: str, input_data: str, output_path: str) -> str:
    """
    @Description 裁剪矢量数据。使用指定的裁剪几何体（如多边形）来裁剪输入的矢量数据（如点、线、多边形）。

    - param clip_geometry :{str} 裁剪几何体的文件路径（如多边形.shp）
    - param input_data    :{str} 需要被裁剪的矢量数据文件路径
    - param output_path   :{str}  裁剪后结果的输出文件路径

    @returns `{ str}` 裁剪后保存的文件路径

    """
    try:
        # 读取裁剪几何体和输入数据
        clip_gdf = gpd.read_file(clip_geometry)
        input_gdf = gpd.read_file(input_data)

        # 使用裁剪几何体对输入数据进行裁剪
        # 注意：这里使用overlay函数，'intersection'表示求交集
        clipped_gdf = gpd.overlay(input_gdf, clip_gdf, how="intersection")

        # 如果裁剪结果为空，则输出提示
        if clipped_gdf.empty:
            print("裁剪结果为空，请检查裁剪几何体和输入数据是否有重叠部分。")
            return ""

        # 保存裁剪后的结果
        clipped_gdf.to_file(output_path, encoding="utf-8")

        return output_path

    except Exception as e:
        print(f"clip_vector_data 失败: {e}")
        return ""

初版​

def polygon_clip_line(polygon_shp: str, line_shp: str, output_shp: str) -> str:
    """
    @Description 复刻arcMap的裁剪功能，这里使用面shp对线段shp进行裁剪
    """
    try:
        clip_polygon = gpd.read_file(polygon_shp)
        line = gpd.read_file(line_shp)

        # 根据项目范围裁剪线段
        clipped_line = gpd.clip(line, clip_polygon)

        # 保存为新的shp文件
        clipped_line.to_file(output_shp, encoding="utf-8")

        return output_shp

    except Exception as e:
        print("polygon_clip_line fail")

        return ""

挖空shp​

使用绿色的方块对紫色的shp进行挖空挖空完成后
def polygon_clip_polygon(polygon_shp: str, clip_rang_shp: str, output_shp: str) -> str:
    """
    @Description 两个相交的部分进行挖空，要挖空的shp支持多个多边形

    - param polygon_shp   :{str} 要挖空的shp
    - param clip_rang_shp :{str} 要挖空的范围，支持多个多边形的shp
    - param output_shp    :{str} 要输出的shp路径
    """
    try:
        # 读取shp文件
        input_polygon = gpd.read_file(polygon_shp)
        subtract_polygon = gpd.read_file(clip_rang_shp)

        res_polygon = input_polygon
        for i, row in subtract_polygon.iterrows():
            each_polygon = row.geometry

            # 进行挖空操作
            res_polygon = res_polygon.difference(each_polygon)

        res_polygon.to_file(output_shp, encoding="utf-8")

        return output_shp
    except Exception as e:
        print("polygon_clip_polygon fail")

        return ""

裁剪​

复刻arcMap的裁剪功能，这里使用面shp对线段shp进行裁剪

def polygon_clip_polygon(polygon_shp: str, clip_rang_shp: str, output_shp: str) -> str:
    """
    @Description 两个相交的部分进行挖空，要挖空的shp支持多个多边形

    - param polygon_shp   :{str} 要挖空的shp
    - param clip_rang_shp :{str} 要挖空的范围，支持多个多边形的shp
    - param output_shp    :{str} 要输出的shp路径
    """
    try:
        # 读取shp文件
        input_polygon = gpd.read_file(polygon_shp)
        subtract_polygon = gpd.read_file(clip_rang_shp)

        res_polygon = input_polygon
        for i, row in subtract_polygon.iterrows():
            each_polygon = row.geometry

            # 进行挖空操作
            res_polygon = res_polygon.difference(each_polygon)

        res_polygon.to_file(output_shp, encoding="utf-8")

        return output_shp
    except Exception as e:
        print("polygon_clip_polygon fail")

        return ""

判断多边形是否交集​

方法1：调用shapely工具箱--[只能处理凸多边形]​

import numpy as np
from shapely.geometry import Polygon # 多边形
 
def Cal_area_2poly(data1,data2):
    """
    任意两个图形的相交面积的计算
    :param data1: 当前物体
    :param data2: 待比较的物体
    :return: 当前物体与待比较的物体的面积交集
    """
 
    poly1 = Polygon(data1).convex_hull # Polygon：多边形对象
    poly2 = Polygon(data2).convex_hull
 
    if not poly1.intersects(poly2):
        inter_area = 0 # 如果两多边形不相交
    else:
        inter_area = poly1.intersection(poly2).area # 相交面积
    return inter_area

data1 = [] # 带比较的第一个物体的顶点坐标
data2 = [] #待比较的第二个物体的顶点坐标
area = Cal_area_2poly(data1,data2)

方法2：将2个图形分别映射到单通道的图片上，进行操作​

import numpy as np
import cv2

def get_iou_2poly(data1, data2, h, w): 
    '''
    get the iou of 2 poly
    '''
    data1 = np.array(data1, np.int32)  ##the area of poly1 > poly2
    data2 = np.array(data2, np.int32)

    mask1 = np.zeros((h, w), np.uint8) 
    mask2 = np.zeros((h, w), np.uint8) 
    # cv2.imwrite("mask.jpg", mask) 

    cv2.fillPoly(mask1, [data1], 1)
    data1_count = (mask1==1).sum()

    cv2.fillPoly(mask2, [data2], 1)
    mask = mask1 + mask2
    inter = (mask==2).sum()  

    if inter > 0:
        return inter
    else:
        return 0 

data1 = []  # 带比较的第一个物体的顶点坐标
data2 = []   #待比较的第二个物体的顶点坐标
h, w = img.shape[:2]
area = Cal_area_2poly(data1,data2，h,w)

获取shp面文件的中心点​

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

from typing import List, Tuple


def calculate_polygon_center(polygon: Polygon) -> Tuple[float, float]:
    """计算多边形的中心点坐标"""
    x, y = polygon.centroid.x, polygon.centroid.y
    return x, y


def process_multipolygon(multipolygon: MultiPolygon) -> List[Tuple[float, float]]:
    """处理 MultiPolygon，并返回其中所有多边形的中心点坐标列表"""
    center_points = [calculate_polygon_center(polygon) for polygon in multipolygon]
    return center_points


def process_shapefile(
    input_shapefile: str, output_shapefile: str = None
) -> List[Tuple[float, float]]:
    """
    处理输入的 shapefile，将各个多边形的中心点导出为一个新的 shapefile（可选）。

    Args:
        input_shapefile (str): 输入的 shapefile 文件路径。
        output_shapefile (str, optional): 输出的 shapefile 文件路径。默认为 None，不导出输出文件。

    Returns:
        List[Tuple[float, float]]: 所有多边形的中心点坐标列表。
    """
    gdf = gpd.read_file(input_shapefile)
    center_points = []

    for idx, row in gdf.iterrows():
        geom = row["geometry"]
        if isinstance(geom, Polygon):
            center_points.append(calculate_polygon_center(geom))
        elif isinstance(geom, MultiPolygon):
            center_points.extend(process_multipolygon(geom))

    center_gdf = gpd.GeoDataFrame(geometry=gpd.points_from_xy(*zip(*center_points)))
    if output_shapefile:
        center_gdf.to_file(output_shapefile, driver="ESRI Shapefile")

    return center_gdf


def calculate_and_export_center(input_target: str, output_shp_path: str) -> None:
    """计算输入 shapefile 中所有多边形的中心点，并导出为一个新的 shapefile"""
    if isinstance(input_target, str):
        gdf = gpd.read_file(input_target)
    else:
        gdf = input_target

    # 计算所有多边形的联合中心点
    center_point = gdf.unary_union.centroid

    # 打印中心点
    # print("中心点坐标：", center_point.x, center_point.y)

    # 创建一个包含中心点的 GeoDataFrame
    center_gdf = gpd.GeoDataFrame(geometry=[center_point], crs=gdf.crs)

    # 导出包含中心点的 shapefile
    center_gdf.to_file(output_shp_path)


# 主函数
if __name__ == "__main__":
    input_shp = "./test/项目范围_面.shp"  # 多个面的shp文件
    output_shp = "./test/output.shp"  # 面文件的所有中心点
    output_point_shp = "./test/center.shp"  # 唯一中心点

    # 调用函数处理多边形 shapefile
    # center_points = process_shapefile(input_shp, output_shp)

    # 调用函数计算中心点并导出为 shapefile
    # calculate_and_export_center(center_points, output_point_shp)

    input_shp_gdf = gpd.read_file(input_shp)
    output_shp_gdf = gpd.read_file(output_shp)

    print(isinstance(output_shp_gdf, MultiPolygon))
    print(isinstance(input_shp_gdf, MultiPolygon))

格式转换​

excel转换成shp​


坐标系获取​


面shp和点shp包含情况判断​

# -*- coding: utf-8 -*-
#
# @Author: CPS
# @email: 373704015@qq.com
# @Date: 2023-10-20 17:55:27.881455
# @Last Modified by: CPS
# @Last Modified time: 2023-10-20 17:55:27.881455
# @file_path "W:\CPS\MyProject\gis-api"
# @Filename "test2.py"
# @Description: 判断面shp和点shp文件的包含情况，根据指定的字段返回字典
#
import os, sys

sys.path.append("..")

import geopandas as gpd
from shapely.geometry import Point


def collect_fields_within_polygons(
    points_shp, polygons_shp, point_fields, polygon_fields
):
    # 读取点集和面集的shapefile文件
    points = gpd.read_file(points_shp)
    polygons = gpd.read_file(polygons_shp)

    # 检查坐标系是否一致
    if points.crs != polygons.crs:
        print(
            "Warning: Coordinate systems are different. It is recommended to use the same coordinate system for accurate spatial operations."
        )

    # 将点集和面集的数据进行合并
    merged_data = gpd.sjoin(points, polygons, op="within")

    # 创建一个字典用于存储结果
    result_dict = {}

    # 遍历每个面，将包含的点添加到字典中
    for index, row in polygons.iterrows():
        page_number = row["PageNumber"]

        # 获取每个点的字段值列表
        point_data = merged_data.loc[
            merged_data["index_right"] == index, point_fields
        ].values.tolist()

        # 将字段值列表添加到结果字典
        result_dict[page_number] = point_data

    return result_dict


# 使用示例
point_fields = ["编号", "X坐标", "Y坐标"]  # 你需要提取的点的字段列表
polygon_fields = ["PageNumber"]  # 你需要提取的面的字段列表


# 使用示例
target = r"谷圩河"


points_file = r"Z:\work\2023\项目\封开县重点中小河流水域岸线保护与利用规划（二期）\出图\shp\岸线规划图层\谷圩河\04谷圩河控制点\一期控制点\右岸控制点.shp"
polygons_file = r"Z:/work/2023/项目/封开县重点中小河流水域岸线保护与利用规划（二期）/出图/shp/显示范围1比2000_fix/显示范围_{0}.shp".format(
    target
)
result = collect_fields_within_polygons(
    points_file, polygons_file, point_fields, polygon_fields
)

for key in result:
    result[key].sort(key=(lambda x: x[0]))
    result[key] = [point_fields] + result[key]

print(result)

SHP文件添加字段​


分割多线段图层​

import geopandas as gpd
from shapely.geometry import LineString, MultiLineString

def split_multilinestring_to_segments(input_shp_path: str, output_shp_path: str = None):
    """
    将一个MultiLineString几何对象拆分为多个LineString，并将每个LineString保存为一个独立的SHP文件。

    参数：
    - input_shp_path: 输入的线SHP文件路径（包含一个MultiLineString几何对象）
    - output_folder: 输出的SHP文件夹路径
    """
    # 读取输入的线SHP文件
    gdf = gpd.read_file(input_shp_path)

    # 确保输入文件包含一个MultiLineString几何对象
    if len(gdf) != 1 or not isinstance(gdf["geometry"].iloc[0], MultiLineString):
        raise ValueError("输入的SHP文件必须包含一个MultiLineString几何对象。")

    # 获取MultiLineString几何对象
    multilinestring = gdf["geometry"].iloc[0]

    # 拆分MultiLineString为多个LineString
    line_segments = [LineString(line.coords) for line in multilinestring]

    # 创建新的GeoDataFrame保存拆分后的线段
    segments_gdf = gpd.GeoDataFrame(geometry=line_segments)

    for col in gdf.columns:
        if col != "geometry":
            segments_gdf[col] = gdf[col].iloc[0]

    # 将所有 LineString 保存为同一个 SHP 文件
    if not output_shp_path is None:
        segments_gdf.to_file(output_shp_path)

    return segments_gdf
# 示例用法
input_shp_path = r"Z:\work\2024\项目\黄埔至南沙东部快速通道沙湾水道特大桥段工程对三沙口水文站水文监测影响分析报告\模型出图\mxd\test.shp"
output_folder = path.join(path.dirname(input_shp_path), "output.shp")

split_multilinestring_to_segments(input_shp_path, output_folder)

获取点集的X和Y 坐标​

使用geopandas​

import geopandas as gpd
import numpy as np

input_shp = r"./clip_range.shp"
gdf = gpd.read_file(raw_shp)
Z = gdf["speed"].values
X = gdf["geometry"].apply(lambda geom: geom.x).values
Y = gdf["geometry"].apply(lambda geom: geom.x).values