泰森多边形计算流域平均降雨量
环境：python3.9
arcgis的操作方法可以参照：https://blog.youkuaiyun.com/yinchoushi8780/article/details/148742106
输入：
流域边界的geojson文件（整个流域边界"id"为"Watershed0"，其他流域id可以随意写），气象站点的坐标。
# sub_catchment.geojson实例
{
"type": "FeatureCollection",
"features": [
{
"type": "Feature",
"geometry": {
"type": "Polygon",
"coordinates": [
[
[
118.61069444444433,
28.41736111111112
],
[
118.61069444444433,
28.41708333333334
]
]
]
},
"properties": {
"id": "Watershed11.1",
}
},
{
"type": "Feature",
"geometry": {
"type": "Polygon",
"coordinates": [
[<类同>]
}}}
# 定义点坐标示例
points_dict = {
"气象站1": {"1": [118.6366, 28.8926]},
"气象站2": {"2": [118.6322, 28.7912]},
"气象站3": {"3": [118.6695, 28.6349]},
"气象站4": {"4": [118.5166, 28.5115]},
"气象站5": {"5": [118.6602, 28.4312]},
"气象站6": {"6": [118.4709, 28.9467]}
}用到的函数：
# 提取泰森多边形的所有边线（包括有限边和无限边）
def extract_all_voronoi_edges(vor):
"""提取泰森多边形的所有边线，包括有限边和无限边"""
edges = []
# 提取有限边
for simplex in vor.ridge_vertices:
if simplex[0] != -1 and simplex[1] != -1:
edge = [vor.vertices[simplex[0]], vor.vertices[simplex[1]]]
edges.append(edge)
# 提取无限边
point_center = vor.points.mean(axis=0)
for pointidx, simplex in zip(vor.ridge_points, vor.ridge_vertices):
if simplex[0] == -1 or simplex[1] == -1:
# 找到有限顶点
i = simplex[0] if simplex[0] != -1 else simplex[1]
if i == -1: # 如果两个顶点都是无限远点，跳过
continue
# 计算无限边的方向向量
t = vor.points[pointidx[1]] - vor.points[pointidx[0]]
t = t / np.linalg.norm(t)
n = np.array([-t[1], t[0]]) # 法向量
# 确定方向
midpoint = vor.points[pointidx].mean(axis=0)
direction = np.sign(np.dot(midpoint - point_center, n)) * n
far_point = vor.vertices[i] + direction * 100 # 延长线
edge = [vor.vertices[i], far_point]
edges.append(edge)
return edges
# 站点编号关联泰森多边形
def assign_point_ids_to_polygons(poly_gdf, points, point_ids):
"""
为每个泰森多边形分配对应的点编号
通过检查点是否在多边形内来确定对应关系
"""
polygon_point_ids = []
for i, polygon in enumerate(poly_gdf.geometry):
assigned = False
# 检查每个点是否在当前多边形内
for j, point in enumerate(points):
point_geom = Point(point[0], point[1])
if polygon.contains(point_geom):
polygon_point_ids.append(point_ids[j])
assigned = True
break
# 如果没有找到包含点的多边形，使用最近点的方法
if not assigned:
print(f"多边形 {i} 内无point点")
return polygon_point_ids
# 生成泰森多边形
def generate_thiessen_polygons(geojson_file, points_dict):
"""
生成泰森多边形并返回流域GeoDataFrame和泰森多边形GeoDataFrame
参数:
geojson_file (str): 流域边界GeoJSON文件路径
points_dict (dict): 点坐标字典，格式如 {"气象站1": {"1": [118.6366, 28.8926]}, ...}
返回:
tuple: (catchment_gdf, poly_gdf_new) 流域GeoDataFrame和泰森多边形GeoDataFrame
"""
# 读取流域边界文件
catchment_gdf = gpd.read_file(geojson_file)
# 提取点坐标和编号
point_coords = []
point_ids = []
point_names = []
for name, info in points_dict.items():
for point_id, coord in info.items():
point_coords.append(coord)
point_ids.append(point_id)
point_names.append(name)
points = np.array(point_coords)
# 生成泰森多边形
vor = Voronoi(points)
# 提取所有边线
all_edges = extract_all_voronoi_edges(vor)
# 将边线转换为LineString对象
voronoi_lines = [LineString(edge) for edge in all_edges]
# 创建线段的GeoDataFrame
line_gdf = gpd.GeoDataFrame({'id': range(len(voronoi_lines))},
geometry=voronoi_lines, crs=catchment_gdf.crs)
# 使用line_gdf中的线段构造成面
# 首先提取所有线段几何并合并
all_lines_from_gdf = unary_union(line_gdf.geometry)
# 1. 生成包含所有点和流域的矩形边界
# 获取所有泰森点的边界
points_bounds = np.array([points[:, 0].min(), points[:, 1].min(),
points[:, 0].max(), points[:, 1].max()])
# 获取流域边界
watershed_bounds = catchment_gdf.total_bounds
# 合并边界，创建包含所有要素的矩形
min_x = min(points_bounds[0], watershed_bounds[0])
min_y = min(points_bounds[1], watershed_bounds[1])
max_x = max(points_bounds[2], watershed_bounds[2])
max_y = max(points_bounds[3], watershed_bounds[3])
# 添加一些缓冲区确保完全包含
buffer = 0.01 # 根据坐标系调整缓冲区大小
min_x -= buffer
min_y -= buffer
max_x += buffer
max_y += buffer
# 2. 创建矩形边界线
rectangle_coords = [(min_x, min_y), (max_x, min_y), (max_x, max_y), (min_x, max_y), (min_x, min_y)]
rectangle_line = LineString(rectangle_coords)
# 3. 合并矩形线和泰森多边形线段
combined_lines_new = unary_union([all_lines_from_gdf, rectangle_line])
# 4. 使用polygonize构造面
polygons_new = list(polygonize(combined_lines_new))
# 5. 使用流域边界裁剪生成的多边形
catchment_polygon = catchment_gdf[catchment_gdf['id'] == 'Watershed0'].geometry.iloc[0]
filtered_polygons_new = []
for polygon in polygons_new:
# 使用流域多边形裁剪
clipped = catchment_polygon.intersection(polygon)
if not clipped.is_empty:
filtered_polygons_new.append(clipped)
poly_gdf_new = gpd.GeoDataFrame(geometry=filtered_polygons_new, crs=catchment_gdf.crs)
# 为多边形分配点编号
assigned_point_ids = assign_point_ids_to_polygons(poly_gdf_new, points, point_ids)
# 添加点编号和点名称列
poly_gdf_new['point_id'] = assigned_point_ids
# 根据point_id查找对应的point_name
assigned_point_names = []
for point_id in assigned_point_ids:
# 在point_ids中查找对应的索引
idx = point_ids.index(point_id)
assigned_point_names.append(point_names[idx])
poly_gdf_new['point_name'] = assigned_point_names
print("多边形与点编号对应关系:")
for i, (point_id, point_name) in enumerate(zip(assigned_point_ids, assigned_point_names)):
print(f"多边形 {i} -> 点编号 {point_id} -> 点名称 {point_name}")
return catchment_gdf, poly_gdf_new
# 计算子流域与泰森多边形的相交分析，并将结果保存为CSV文件
def calculate_intersection_analysis(sub_catchments, poly_gdf_new, output_path='thiessen_results.csv'):
"""
计算子流域与泰森多边形的相交分析，并将结果保存为CSV文件
参数:
sub_catchments (GeoDataFrame): 子流域数据
poly_gdf_new (GeoDataFrame): 泰森多边形数据
output_path (str): 输出CSV文件的路径，默认为'intersection_results.csv'
返回:
str: 保存的CSV文件路径
"""
# 创建结果属性表
intersection_results = []
# 为每个子流域计算与泰森多边形的相交关系
for idx, catchment_row in sub_catchments.iterrows():
catchment_geom = catchment_row.geometry
# 获取真实的流域ID字段
actual_catchment_id = catchment_row['id'] if 'id' in sub_catchments.columns else idx
# 获取原始流域面积（使用.area方法计算，与S1和S2单位一致）
s = catchment_geom.area
# 遍历所有泰森多边形
for poly_idx, poly_row in poly_gdf_new.iterrows():
poly_geom = poly_row.geometry
# 确保正确获取 point_id 和 point_name
point_id = poly_row['point_id'] if 'point_id' in poly_row else None
point_name = poly_row['point_name'] if 'point_name' in poly_row else None
# 计算与当前流域的交集
intersection = catchment_geom.intersection(poly_geom)
if not intersection.is_empty and intersection.area > 0:
# 计算交集面积(S1)
s1 = intersection.area
# 计算完整泰森多边形面积(S2)
s2 = poly_geom.area
# 计算权重
weight = s1 / s if s > 0 else 0
# 添加到结果
intersection_results.append({
'流域id': actual_catchment_id,
'泰森多边形的point_id': point_id,
'站点名称': point_name,
'相交后的面积S1': s1,
'原始流域面积': s, # 使用.area计算，单位与S1和S2一致
'对应泰森多边形面积S2': s2,
'权重': weight
})
# 转换为GeoDataFrame并保存结果
if intersection_results:
result_gdf = gpd.GeoDataFrame(intersection_results, columns=[
'流域id', '泰森多边形的point_id', '站点名称', '相交后的面积S1',
'原始流域面积', '对应泰森多边形面积S2', '权重'
])
print("\n相交分析结果表:")
print(result_gdf.to_string(index=False))
# 保存结果到CSV文件
# result_gdf.to_csv(output_path, index=False, encoding='utf-8-sig')
# print(f"相交分析结果已保存至: {output_path}")
rainfall_weight = {
'basin_id': result_gdf['流域id'].tolist(),
'station': result_gdf['站点名称'].tolist(),
'rainfall_weight': result_gdf['权重'].tolist()
}
return rainfall_weight
else:
print("没有找到相交区域")
return None主程序：
我将整个流域的外边界记录了'id'为'Watershed0'的Polygon了，把它和子流域边界一起保存进了sub_catchment.geojson文件，所以在做相交分析的时候需要先把Watershed0去除了。
import numpy as np
from scipy.spatial import Voronoi
from shapely.geometry import Polygon, LineString, Point
import geopandas as gpd
from shapely.ops import unary_union, polygonize
import csv
# 定义点坐标
points_dict = {
"气象站1": {"1": [118.6366, 28.8926]},
"气象站2": {"2": [118.6322, 28.7912]},
"气象站3": {"3": [118.6695, 28.6349]},
"气象站4": {"4": [118.5166, 28.5115]},
"气象站5": {"5": [118.6602, 28.4312]},
"气象站6": {"6": [118.4709, 28.9467]}
}
thiessen_output = 'intersection_results.csv'
# 调用函数
catchment_gdf, poly_gdf_new = generate_thiessen_polygons('sub_catchment.geojson', points_dict)
# 获取除去Watershed0以外的子流域
sub_catchments = catchment_gdf[catchment_gdf['id'] != 'Watershed0'].copy()
# 子流域与泰森多边形的相交分析，计算权重，并将结果保存为CSV文件
calculate_intersection_analysis(sub_catchments, poly_gdf_new, thiessen_output)结果：
此部分面积单位为平方度，我没有进行投影计算面积，因为我只需要得到权重，精度已经够用了。
后面将对应的站点的降水乘以这个权重就可以得到流域平均降雨量了
完整代码：
import numpy as np
from scipy.spatial import Voronoi
from shapely.geometry import Polygon, LineString, Point
import geopandas as gpd
from shapely.ops import unary_union, polygonize
import matplotlib.pyplot as plt
# 提取泰森多边形的所有边线（包括有限边和无限边）
def extract_all_voronoi_edges(vor):
"""提取泰森多边形的所有边线，包括有限边和无限边"""
edges = []
# 提取有限边
for simplex in vor.ridge_vertices:
if simplex[0] != -1 and simplex[1] != -1:
edge = [vor.vertices[simplex[0]], vor.vertices[simplex[1]]]
edges.append(edge)
# 提取无限边
point_center = vor.points.mean(axis=0)
for pointidx, simplex in zip(vor.ridge_points, vor.ridge_vertices):
if simplex[0] == -1 or simplex[1] == -1:
# 找到有限顶点
i = simplex[0] if simplex[0] != -1 else simplex[1]
if i == -1: # 如果两个顶点都是无限远点，跳过
continue
# 计算无限边的方向向量
t = vor.points[pointidx[1]] - vor.points[pointidx[0]]
t = t / np.linalg.norm(t)
n = np.array([-t[1], t[0]]) # 法向量
# 确定方向
midpoint = vor.points[pointidx].mean(axis=0)
direction = np.sign(np.dot(midpoint - point_center, n)) * n
far_point = vor.vertices[i] + direction * 100 # 延长线
edge = [vor.vertices[i], far_point]
edges.append(edge)
return edges
# 站点编号关联泰森多边形
def assign_point_ids_to_polygons(poly_gdf, points, point_ids):
"""
为每个泰森多边形分配对应的点编号
通过检查点是否在多边形内来确定对应关系
"""
polygon_point_ids = []
for i, polygon in enumerate(poly_gdf.geometry):
assigned = False
# 检查每个点是否在当前多边形内
for j, point in enumerate(points):
point_geom = Point(point[0], point[1])
if polygon.contains(point_geom):
polygon_point_ids.append(point_ids[j])
assigned = True
break
# 如果没有找到包含点的多边形，使用最近点的方法
if not assigned:
print(f"多边形 {i} 内无point点")
return polygon_point_ids
# 生成泰森多边形
def generate_thiessen_polygons(geojson_file, points_dict):
"""
生成泰森多边形并返回流域GeoDataFrame和泰森多边形GeoDataFrame
参数:
geojson_file (str): 流域边界GeoJSON文件路径
points_dict (dict): 点坐标字典，格式如 {"气象站1": {"1": [118.6366, 28.8926]}, ...}
返回:
tuple: (catchment_gdf, poly_gdf_new) 流域GeoDataFrame和泰森多边形GeoDataFrame
"""
# 读取流域边界文件
catchment_gdf = gpd.read_file(geojson_file)
# 提取点坐标和编号
point_coords = []
point_ids = []
point_names = []
for name, info in points_dict.items():
for point_id, coord in info.items():
point_coords.append(coord)
point_ids.append(point_id)
point_names.append(name)
points = np.array(point_coords)
# 生成泰森多边形
vor = Voronoi(points)
# 提取所有边线
all_edges = extract_all_voronoi_edges(vor)
# 将边线转换为LineString对象
voronoi_lines = [LineString(edge) for edge in all_edges]
# 创建线段的GeoDataFrame
line_gdf = gpd.GeoDataFrame({'id': range(len(voronoi_lines))},
geometry=voronoi_lines, crs=catchment_gdf.crs)
# 使用line_gdf中的线段构造成面
# 首先提取所有线段几何并合并
all_lines_from_gdf = unary_union(line_gdf.geometry)
# 1. 生成包含所有点和流域的矩形边界
# 获取所有泰森点的边界
points_bounds = np.array([points[:, 0].min(), points[:, 1].min(),
points[:, 0].max(), points[:, 1].max()])
# 获取流域边界
watershed_bounds = catchment_gdf.total_bounds
# 合并边界，创建包含所有要素的矩形
min_x = min(points_bounds[0], watershed_bounds[0])
min_y = min(points_bounds[1], watershed_bounds[1])
max_x = max(points_bounds[2], watershed_bounds[2])
max_y = max(points_bounds[3], watershed_bounds[3])
# 添加一些缓冲区确保完全包含
buffer = 0.01 # 根据坐标系调整缓冲区大小
min_x -= buffer
min_y -= buffer
max_x += buffer
max_y += buffer
# 2. 创建矩形边界线
rectangle_coords = [(min_x, min_y), (max_x, min_y), (max_x, max_y), (min_x, max_y), (min_x, min_y)]
rectangle_line = LineString(rectangle_coords)
# 3. 合并矩形线和泰森多边形线段
combined_lines_new = unary_union([all_lines_from_gdf, rectangle_line])
# 4. 使用polygonize构造面
polygons_new = list(polygonize(combined_lines_new))
# 5. 使用流域边界裁剪生成的多边形
catchment_polygon = catchment_gdf[catchment_gdf['id'] == 'Watershed0'].geometry.iloc[0]
filtered_polygons_new = []
for polygon in polygons_new:
# 使用流域多边形裁剪
clipped = catchment_polygon.intersection(polygon)
if not clipped.is_empty:
filtered_polygons_new.append(clipped)
poly_gdf_new = gpd.GeoDataFrame(geometry=filtered_polygons_new, crs=catchment_gdf.crs)
# 为多边形分配点编号
assigned_point_ids = assign_point_ids_to_polygons(poly_gdf_new, points, point_ids)
# 添加点编号和点名称列
poly_gdf_new['point_id'] = assigned_point_ids
# 根据point_id查找对应的point_name
assigned_point_names = []
for point_id in assigned_point_ids:
# 在point_ids中查找对应的索引
idx = point_ids.index(point_id)
assigned_point_names.append(point_names[idx])
poly_gdf_new['point_name'] = assigned_point_names
print("多边形与点编号对应关系:")
for i, (point_id, point_name) in enumerate(zip(assigned_point_ids, assigned_point_names)):
print(f"多边形 {i} -> 点编号 {point_id} -> 点名称 {point_name}")
return catchment_gdf, poly_gdf_new
# 计算子流域与泰森多边形的相交分析，并将结果保存为CSV文件
def calculate_intersection_analysis(sub_catchments, poly_gdf_new, output_path='thiessen_results.csv'):
"""
计算子流域与泰森多边形的相交分析，并将结果保存为CSV文件
参数:
sub_catchments (GeoDataFrame): 子流域数据
poly_gdf_new (GeoDataFrame): 泰森多边形数据
output_path (str): 输出CSV文件的路径，默认为'intersection_results.csv'
返回:
str: 保存的CSV文件路径
"""
# 创建结果属性表
intersection_results = []
# 为每个子流域计算与泰森多边形的相交关系
for idx, catchment_row in sub_catchments.iterrows():
catchment_geom = catchment_row.geometry
# 获取真实的流域ID字段
actual_catchment_id = catchment_row['id'] if 'id' in sub_catchments.columns else idx
# 获取原始流域面积（使用.area方法计算，与S1和S2单位一致）
s = catchment_geom.area
# 遍历所有泰森多边形
for poly_idx, poly_row in poly_gdf_new.iterrows():
poly_geom = poly_row.geometry
# 确保正确获取 point_id 和 point_name
point_id = poly_row['point_id'] if 'point_id' in poly_row else None
point_name = poly_row['point_name'] if 'point_name' in poly_row else None
# 计算与当前流域的交集
intersection = catchment_geom.intersection(poly_geom)
if not intersection.is_empty and intersection.area > 0:
# 计算交集面积(S1)
s1 = intersection.area
# 计算完整泰森多边形面积(S2)
s2 = poly_geom.area
# 计算权重
weight = s1 / s if s > 0 else 0
# 添加到结果
intersection_results.append({
'流域id': actual_catchment_id,
'泰森多边形的point_id': point_id,
'站点名称': point_name,
'相交后的面积S1': s1,
'原始流域面积': s, # 使用.area计算，单位与S1和S2一致
'对应泰森多边形面积S2': s2,
'权重': weight
})
# 转换为GeoDataFrame并保存结果
if intersection_results:
result_gdf = gpd.GeoDataFrame(intersection_results, columns=[
'流域id', '泰森多边形的point_id', '站点名称', '相交后的面积S1',
'原始流域面积', '对应泰森多边形面积S2', '权重'
])
print("\n相交分析结果表:")
print(result_gdf.to_string(index=False))
# 保存结果到CSV文件
# result_gdf.to_csv(output_path, index=False, encoding='utf-8-sig')
# print(f"相交分析结果已保存至: {output_path}")
rainfall_weight = {
'basin_id': result_gdf['流域id'].tolist(),
'station': result_gdf['站点名称'].tolist(),
'rainfall_weight': result_gdf['权重'].tolist()
}
return rainfall_weight
else:
print("没有找到相交区域")
return None
#--------------------------------------------------------------------------------------------
# 定义点坐标
points_dict = {
"气象站1": {"1": [118.6366, 28.8926]},
"气象站2": {"2": [118.6322, 28.7912]},
"气象站3": {"3": [118.6695, 28.6349]},
"气象站4": {"4": [118.5166, 28.5115]},
"气象站5": {"5": [118.6602, 28.4312]},
"气象站6": {"6": [118.4709, 28.9467]}
}
# thiessen_output = 'intersection_results.csv'
# 调用函数
catchment_gdf, poly_gdf_new = generate_thiessen_polygons('sub_catchment.geojson', points_dict)
# 获取除去Watershed0以外的子流域
sub_catchments = catchment_gdf[catchment_gdf['id'] != 'Watershed0'].copy()
# 结果
rainfall_weight = calculate_intersection_analysis(sub_catchments, poly_gdf_new)
print(rainfall_weight)补充：
一个shp转geojson文件的方法
import geopandas as gpd
def shp_to_geojson_advanced(shp_file_path, geojson_file_path, from_encoding='gbk', to_encoding='utf-8'):
"""
将 shapefile 转换为 GeoJSON 格式（支持编码处理）
参数:
shp_file_path (str): 输入的 shapefile 文件路径
geojson_file_path (str): 输出的 GeoJSON 文件路径
from_encoding (str): 输入 shapefile 的编码格式，默认为 'gbk'
to_encoding (str): 输出 GeoJSON 的编码格式，默认为 'utf-8'
返回:
str: 输出的 GeoJSON 文件路径
"""
try:
# 读取 shapefile 文件
gdf = gpd.read_file(shp_file_path, encoding=from_encoding)
# 保存为 GeoJSON 格式
gdf.to_file(geojson_file_path, driver='GeoJSON', encoding=to_encoding)
print(f"已成功将 {shp_file_path} 转换为 {geojson_file_path}")
return geojson_file_path
except Exception as e:
print(f"转换过程中出现错误: {e}")
return None
# 使用示例
# shp_to_geojson_advanced('input.shp', 'output.geojson')
扫一扫
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
