### Python GIS分析：核心工具与实战代码 GIS空间分析的核心操作包括**空间叠加分析**、**缓冲区分析**、**可视域分析**等。Python生态通过`geopandas`、`shapely`、`rasterio`、`arcpy`等库提供了完整的解决方案。下表对比了主流库的核心功能与适用场景[ref_3][ref_4]： | 分析类型 | 核心库/工具 | 主要功能 | 适用场景 | | :--- | :--- | :--- | :--- | | **矢量叠加分析** | GeoPandas / Shapely | 交集、并集、差集、标识、联合 | 土地利用变化、行政区划合并 | | **栅格可视分析** | Rasterio / NumPy | 可视率、视线分析、地形计算 | 通信基站选址、景观可视性评估 | | **专业GIS工具** | ArcPy (ArcGIS) | 完整的GP工具链、高级空间统计 | 企业级GIS工作流、复杂空间建模 | | **空间关系与度量** | Shapely / PySAL | 缓冲区、距离计算、空间连接 | 设施服务范围、邻近度分析 | #### 1. 矢量数据叠加分析（以GeoPandas为例） 叠加分析是GIS的核心，用于计算两个或多个空间数据层的几何交集、并集等关系。以下代码演示了如何使用`geopandas`进行常见的叠加操作[ref_4]： ```python import geopandas as gpd import matplotlib.pyplot as plt # 1. 加载数据（示例为两个相交的面图层） gdf_parcels = gpd.read_file('parcels.shp') # 地块数据 gdf_floodzone = gpd.read_file('flood_zone.shp') # 洪泛区数据 # 确保使用相同的坐标系 gdf_floodzone = gdf_floodzone.to_crs(gdf_parcels.crs) # 2. 叠加分析操作 # a) 交集 (Intersection): 提取位于洪泛区内的地块部分 intersection = gpd.overlay(gdf_parcels, gdf_floodzone, how='intersection') # b) 并集 (Union): 合并两个图层的所有区域 union = gpd.overlay(gdf_parcels, gdf_floodzone, how='union') # c) 差集 (Difference): 提取不在洪泛区内的地块部分 difference = gpd.overlay(gdf_parcels, gdf_floodzone, how='difference') # d) 标识 (Identity): 保留第一个图层的所有要素，并用第二个图层进行切割和属性附加 identity = gpd.overlay(gdf_parcels, gdf_floodzone, how='identity') # 注：GeoPandas的`overlay`函数默认实现了这些操作[ref_4]。 # 3. 空间连接 (Spatial Join): 基于空间关系关联属性 # 为每个地块附加其所在洪泛区的属性（如洪水风险等级） joined = gpd.sjoin(gdf_parcels, gdf_floodzone, how='left', predicate='intersects') # 4. 可视化结果 fig, axes = plt.subplots(2, 3, figsize=(15, 10)) gdf_parcels.plot(ax=axes[0,0], color='lightgray', edgecolor='black') axes[0,0].set_title('原始地块') gdf_floodzone.plot(ax=axes[0,1], color='blue', alpha=0.5) axes[0,1].set_title('洪泛区') intersection.plot(ax=axes[0,2], color='red') axes[0,2].set_title('交集（受淹地块）') difference.plot(ax=axes[1,0], color='green') axes[1,0].set_title('差集（安全地块）') union.plot(ax=axes[1,1], color='purple', alpha=0.7) axes[1,1].set_title('并集') identity.plot(ax=axes[1,2], column='FLOOD_ZONE', legend=True, cmap='Reds') # 假设洪泛区有FLOOD_ZONE字段 axes[1,2].set_title('标识分析') plt.tight_layout() plt.show() ``` #### 2. 栅格数据分析（以可视率分析为例） 可视率分析用于计算从某观察点出发，地形上各点的可见性。其核心算法是**连线比较法**，即判断观察点到目标点的连线是否被中间的高程点所阻挡[ref_1]。 ```python import rasterio import numpy as np import matplotlib.pyplot as plt from math import sqrt def calculate_viewshed(dem_path, observer_point, observer_height=1.7, max_dist=5000): """ 计算数字高程模型（DEM）上某点的可视域。 :param dem_path: DEM文件路径 :param observer_point: (x, y) 观察点坐标（地图投影坐标） :param observer_height: 观察者离地高度（米） :param max_dist: 最大分析距离（米） :return: 二维布尔数组，True表示该像元可见。 """ with rasterio.open(dem_path) as src: dem = src.read(1) # 读取第一波段高程数据 profile = src.profile transform = src.transform # 将地理坐标转换为行列号 obs_row, obs_col = src.index(observer_point[0], observer_point[1]) height_obs = dem[obs_row, obs_col] + observer_height rows, cols = dem.shape viewshed = np.zeros((rows, cols), dtype=bool) # 遍历DEM中的每个像元（此处为简化示例，实际应用需优化性能） for r in range(rows): for c in range(cols): if (r == obs_row and c == obs_col): viewshed[r, c] = True # 观察点自身可见 continue # 计算两点间的距离和方向向量 dist = sqrt((r - obs_row)**2 + (c - obs_col)**2) * transform.a # 假设像元为正方形 if dist > max_dist: continue # 生成观察点到目标点连线经过的像元坐标（ Bresenham 直线算法简化） line_cells = bresenham_line(obs_row, obs_col, r, c) visible = True max_tan_angle = -np.inf for lr, lc in line_cells[1:-1]: # 跳过起点和终点 if 0 <= lr < rows and 0 <= lc < cols: # 计算中间点相对于观察点的高度角 cell_dist = sqrt((lr - obs_row)**2 + (lc - obs_col)**2) * transform.a height_cell = dem[lr, lc] tan_angle = (height_cell - height_obs) / cell_dist if cell_dist > 0 else np.inf # 如果中间点的高度角大于当前最大高度角，则更新并判断 if tan_angle > max_tan_angle: max_tan_angle = tan_angle # 计算目标点的最终高度角 tan_target = (dem[r, c] - height_obs) / dist if dist > 0 else np.inf if max_tan_angle > tan_target: visible = False break viewshed[r, c] = visible return viewshed, profile def bresenham_line(x0, y0, x1, y1): """简化版的Bresenham直线算法，返回连线经过的像素坐标列表。""" points = [] dx = abs(x1 - x0) dy = abs(y1 - y0) sx = 1 if x0 < x1 else -1 sy = 1 if y0 < y1 else -1 err = dx - dy while True: points.append((x0, y0)) if x0 == x1 and y0 == y1: break e2 = 2 * err if e2 > -dy: err -= dy x0 += sx if e2 < dx: err += dx y0 += sy return points # 使用示例 if __name__ == '__main__': dem_file = 'path/to/your/dem.tif' observer = (456200.0, 4428000.0) # 示例UTM坐标 visible_cells, dem_profile = calculate_viewshed(dem_file, observer, max_dist=3000) # 可视化 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,5)) with rasterio.open(dem_file) as src: dem_img = src.read(1) ax1.imshow(dem_img, cmap='terrain') ax1.scatter(observer[1], observer[0], c='red', s=50, marker='^') # 注意行列号与坐标的转换 ax1.set_title('DEM地形') ax2.imshow(visible_cells, cmap='binary') ax2.scatter(observer[1], observer[0], c='red', s=50, marker='^') ax2.set_title('可视域分析结果（白色为可见）') plt.show() ``` **关键点**：实际应用中，上述双循环算法效率较低，需通过**向量化计算**（使用NumPy）或**并行计算**进行优化[ref_1]。此外，专业分析还需考虑**地球曲率**和**大气折射**的修正。 #### 3. 使用ArcPy进行专业叠加分析 在ArcGIS环境中，`arcpy`提供了与桌面工具完全对应的功能，适合自动化处理企业级GIS工作流[ref_2][ref_5][ref_6]。 ```python import arcpy import os # 设置工作空间和许可 arcpy.env.workspace = r"C:\Project\Data.gdb" arcpy.CheckOutExtension("Spatial") # 1. 叠加分析（以“联合”和“交集取反”为例） input_features = ["Parcels", "FloodZones"] output_union = "Parcels_Flood_Union" output_erase = "Parcels_Outside_Flood" # a) 联合 (Union): 计算所有输入要素的几何并集 arcpy.analysis.Union(input_features, output_union) # 此操作将两个图层的所有区域合并，并保留所有属性[ref_6]。 # b) 交集取反 (Erase): 从输入要素中擦除与擦除要素重叠的部分 # 即找出不在洪泛区内的地块 arcpy.analysis.Erase("Parcels", "FloodZones", output_erase) # 此操作是叠加分析中“差集”的一种实现[ref_5]。 # 2. 标识 (Identity): 保留输入要素的所有部分，并用标识要素进行切割和属性附加 output_identity = "Parcels_Identity_Flood" arcpy.analysis.Identity("Parcels", "FloodZones", output_identity) # 结果中，与洪泛区相交的地块被分割并附加洪泛区属性，不相交的保持原状[ref_2]。 # 3. 空间连接 (Spatial Join) output_join = "Parcels_Joined" arcpy.analysis.SpatialJoin("Parcels", "FloodZones", output_join, join_operation="JOIN_ONE_TO_ONE", match_option="INTERSECT") print("叠加分析完成。") ``` #### 4. 缓冲区分析与空间度量 缓冲区分析用于创建要素周围一定距离内的区域，是邻近度分析的基础[ref_4]。 ```python import geopandas as gpd from shapely.geometry import Point, LineString # 创建示例点、线数据 points = gpd.GeoSeries([Point(0, 0), Point(2, 3)]) lines = gpd.GeoSeries([LineString([(1,1), (4,5)])]) # 创建缓冲区 point_buffers = points.buffer(distance=1) # 1单位距离的缓冲区 line_buffers = lines.buffer(distance=0.5) # 计算面积和距离 area = point_buffers.area distance_between = points[0].distance(points[1]) print(f"缓冲区面积: {area.values}") print(f"点间距离: {distance_between}") ``` **总结**：Python GIS分析的选择取决于具体任务。对于快速原型开发、开源环境或轻量级分析，`geopandas`和`rasterio`组合是首选[ref_3][ref_4]；对于需要执行复杂视线分析（如可视率）的场景，需结合`NumPy`进行底层数组运算[ref_1]；而在已部署ArcGIS的企业环境中，使用`arcpy`调用成熟的GP工具则是最高效、最稳定的方式[ref_2][ref_5][ref_6]。所有操作都需特别注意**坐标参考系统（CRS）的一致性**，这是保证空间分析结果正确性的前提[ref_4]。