如何使用Python将Shapefile文件拆分为两个独立文件

问题描述 #

我有一个Shapefile文件，希望将其拆分为两个不同的Shapefile文件（拆分出指定经纬度范围的区域和剩余区域）。我尝试的代码输出的是由多个点组成的连线，无法得到预期的面/线要素，代码如下：

import geopandas as gpd ; import shapefile ; import numpy as np ; import pandas as pd
from shapely.geometry import Polygon, Point
import matplotlib.pyplot as plt ;  

lat_1=6 ; lat_2=20 ; lon_1=70 ; lon_2=90
A=gpd.read_file(shp_file)
points=[]
for shape_rec in A.shapeRecords():
    #print(shape_rec)
     pts = pd.DataFrame(np.array(shape_rec.shape.points))
     pts.columns=['lon','lat']
     points.append(pts) 
points_1=pd.concat(points)

points_2=points_1[(points_1.lon >=lon_1)&(points_1.lon <=lon_2) & (points_1.lat >=lat_1)&(points_1.lat <=lat_2)] 


gdf = gpd.GeoDataFrame( points_2, geometry=gpd.points_from_xy(points_2.lon , points_2.lat ), crs='EPSG:4326')
gdf.drop(['lat', 'lon'], axis=1, inplace=True)  # optional
gdf.to_file(main+'shp_files/test.shp') 

问题原因 #

你这段代码错误地把Shapefile中的几何要素拆解成了单个点，筛选后保存的是点要素，自然只会得到点连线，而不是保留原有的多边形/线几何类型。

解决方案 #

下面提供两种符合需求的实现方式，根据你的实际场景选择：

方式1：筛选完全落在指定范围的整要素

如果目标是保留完全处于指定经纬度范围内的完整要素，直接基于GeoDataFrame的几何属性筛选即可：

import geopandas as gpd

# 配置参数
lat_min, lat_max = 6, 20
lon_min, lon_max = 70, 90
shp_path = "你的原始Shapefile路径"
output_path = "输出文件夹路径/"

# 读取原始数据
gdf = gpd.read_file(shp_path)

# 筛选质心在指定范围内的要素（也可改用geometry.bounds判断完全包含）
filtered = gdf[
    (gdf.geometry.centroid.x >= lon_min) &
    (gdf.geometry.centroid.x <= lon_max) &
    (gdf.geometry.centroid.y >= lat_min) &
    (gdf.geometry.centroid.y <= lat_max)
]

# 筛选剩余要素
remaining = gdf[~gdf.index.isin(filtered.index)]

# 保存结果
filtered.to_file(f"{output_path}filtered_area.shp")
remaining.to_file(f"{output_path}remaining_area.shp")

方式2：裁剪出与范围相交的要素部分

如果目标是提取要素与指定范围重叠的部分（比如一个多边形只有局部在范围内，需要切割出这部分），用geopandas.clip工具：

import geopandas as gpd
from shapely.geometry import Polygon

# 配置参数
lat_min, lat_max = 6, 20
lon_min, lon_max = 70, 90
shp_path = "你的原始Shapefile路径"
output_path = "输出文件夹路径/"

# 创建裁剪范围的多边形
clip_poly = Polygon([
    (lon_min, lat_min), (lon_max, lat_min),
    (lon_max, lat_max), (lon_min, lat_max), (lon_min, lat_min)
])
clip_gdf = gpd.GeoDataFrame({"geometry": [clip_poly]}, crs="EPSG:4326")

# 读取原始数据，确保CRS匹配
gdf = gpd.read_file(shp_path).to_crs("EPSG:4326")

# 裁剪相交区域
clipped = gpd.clip(gdf, clip_gdf)

# 提取剩余区域
remaining = gdf.difference(clip_poly)
remaining = gpd.GeoDataFrame(geometry=remaining, crs=gdf.crs)
remaining = remaining[~remaining.geometry.is_empty]  # 移除空几何

# 保存结果
clipped.to_file(f"{output_path}clipped_area.shp")
remaining.to_file(f"{output_path}remaining_area.shp")

注意事项 #

• 确保原始Shapefile的坐标参考系（CRS）是EPSG:4326（经纬度），如果不是，先通过gdf.to_crs("EPSG:4326")转换。
• 方式1适合筛选完整要素，方式2适合切割要素局部区域，按需选择。

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