用PythonArcGIS手把手教你计算土地利用强度指数附完整代码与数据土地利用强度指数Land Use Intensity Index, LUI是量化人类活动对自然环境干扰程度的重要指标。对于地理信息科学、生态学和城市规划领域的研究者来说掌握LUI的计算方法不仅能评估区域开发强度还能为生态保护、国土空间规划提供数据支撑。本文将带你从零开始通过Python和ArcGIS的组合完成从原始栅格数据到可视化成果的全流程分析。1. 数据准备与预处理计算LUI的第一步是获取并预处理土地利用数据。目前常用的全球土地利用数据集包括Globeland30、FROM-GLC等分辨率从30米到300米不等。以Globeland30为例其将地表覆盖分为10大类类别代码土地利用类型典型强度赋值范围10耕地4-720林地2-330草地2-340灌木地1-250湿地1-260水域1-370人造地表5-1080裸地0-190冰雪0-1注意强度赋值需根据研究区域特点调整城市地区建设用地的权重通常高于农业区在ArcGIS中预处理数据的典型步骤投影转换确保所有数据使用统一坐标系统裁剪研究区使用掩膜提取目标区域重分类将原始分类转换为强度值# ArcPy代码示例重分类土地利用数据 import arcpy from arcpy.sa import * # 设置工作环境 arcpy.env.workspace D:/LUI_Project/Data arcpy.env.overwriteOutput True # 定义重分类规则旧值:新值 remap RemapValue([[10,4],[20,2],[30,2],[40,1], [50,1],[60,1],[70,8],[80,0],[90,0]]) # 执行重分类 out_reclass Reclassify(landuse.tif, Value, remap) out_reclass.save(landuse_intensity.tif)2. LUI指数计算原理与实现土地利用强度指数的核心公式为$$ LUI \frac{\sum_{i1}^{n}(A_i \times C_i)}{A} \times 100 $$其中$A_i$第i类土地利用类型的面积$C_i$对应类型的强度赋值$A$区域总面积$n$土地利用类型数量在Python中实现该计算的高效方法是结合rasterio和numpyimport rasterio import numpy as np def calculate_lui(intensity_raster): with rasterio.open(intensity_raster) as src: data src.read(1) # 计算各类面积 unique, counts np.unique(data, return_countsTrue) # 转换为字典便于查询 area_dict dict(zip(unique, counts)) total_area data.size # 假设我们已经定义了强度权重字典 weight_dict {0:0, 1:1, 2:2, 4:4, 8:8} # 计算分子部分 numerator sum(area_dict.get(k,0)*weight_dict.get(k,0) for k in area_dict) # 计算LUI lui (numerator / total_area) * 100 return lui对于需要分区统计的情况如按行政区计算可以结合zonal statistics# 分区统计示例 from rasterstats import zonal_stats def zonal_lui(zone_shp, intensity_raster): stats zonal_stats(zone_shp, intensity_raster, statssum count, categoricalTrue) # 后续处理逻辑... return zone_results3. 空间分析与进阶处理基础LUI计算完成后通常需要进行以下空间分析移动窗口分析计算局部区域的LUI变化时间序列比较分析多年份LUI变化趋势空间自相关检验评估LUI的空间聚集特征移动窗口分析的实现示例import numpy as np from scipy.ndimage import uniform_filter def moving_window_lui(raster_path, window_size3): with rasterio.open(raster_path) as src: data src.read(1) # 转换权重值 weight_matrix np.where(data0,0, np.where(data1,1, np.where(data2,2, np.where(data4,4,8)))) # 计算移动平均 window (window_size, window_size) sum_arr uniform_filter(weight_matrix, sizewindow, modeconstant) * (window_size**2) count_arr uniform_filter(np.ones_like(weight_matrix), sizewindow, modeconstant) * (window_size**2) # 计算局部LUI local_lui (sum_arr / count_arr) * 100 return local_lui对于时间序列分析建议使用xarray库处理多维数据import xarray as xr def multiyear_lui_analysis(file_list, year_list): # 创建空数据集 lui_ds xr.Dataset() for f, year in zip(file_list, year_list): lui calculate_lui(f) lui_ds[flui_{year}] xr.DataArray(lui) # 计算变化率 lui_ds[change_rate] (lui_ds[flui_{year_list[-1]}] - lui_ds[flui_{year_list[0]}]) / len(year_list) return lui_ds4. 结果可视化与解读高质量的可视化能显著提升研究成果的传达效果。推荐以下几种可视化方式1. 分级渲染地图在ArcGIS中创建LUI分级图右键点击图层 → 属性 → 符号系统选择分类渲染方式设置适当的分级数和色带添加图例和比例尺2. 动态交互图表使用matplotlib创建组合图表import matplotlib.pyplot as plt def plot_lui_trend(lui_data, years): fig, (ax1, ax2) plt.subplots(2, 1, figsize(10,8)) # 折线图 ax1.plot(years, lui_data, bo-) ax1.set_title(LUI Temporal Trend) ax1.set_ylabel(LUI Index) # 柱状图 changes np.diff(lui_data) ax2.bar(years[1:], changes, colornp.where(changes0,r,g)) ax2.set_title(Annual Change) ax2.axhline(0, colorgray) plt.tight_layout() return fig3. 空间变化热点图结合空间统计结果绘制热点图from esda.moran import Moran_Local import libpysal as lps def hotspot_analysis(lui_raster): # 转换为空间权重矩阵 w lps.weights.Queen.from_dataframe(gdf) # 计算局部Morans I moran_loc Moran_Local(lui_values, w) # 分类热点冷点 hotspots (moran_loc.q1) * (moran_loc.p_sim 0.05) coldspots (moran_loc.q3) * (moran_loc.p_sim 0.05) return hotspots, coldspots5. 完整工作流与质量控制为确保分析结果可靠建议遵循以下质量控制流程数据校验阶段检查原始数据的投影一致性验证分类数据的精度混淆矩阵确认重分类规则的科学性计算过程检查中间结果抽样验证边界效应处理异常值检测与处理结果验证方法与统计年鉴数据对比实地抽样验证不同分辨率结果交叉检验完整工作流示例代码def full_lui_workflow(landuse_raster, weight_dict, output_dir): 完整的LUI计算工作流 # 1. 重分类 intensity_raster os.path.join(output_dir, intensity.tif) reclassify_landuse(landuse_raster, weight_dict, intensity_raster) # 2. 计算全局LUI global_lui calculate_lui(intensity_raster) # 3. 分区统计 if zone_shapefile: zone_stats zonal_lui(zone_shapefile, intensity_raster) # 4. 移动窗口分析 window_lui moving_window_lui(intensity_raster, window_size5) # 5. 保存结果 with rasterio.open(intensity_raster) as src: profile src.profile profile.update(dtyperasterio.float32) with rasterio.open(os.path.join(output_dir, lui_window.tif), w, **profile) as dst: dst.write(window_lui.astype(rasterio.float32), 1) return { global_lui: global_lui, zone_stats: zone_stats, window_result: os.path.join(output_dir, lui_window.tif) }实际项目中遇到的典型问题及解决方案边缘效应扩大研究区范围后裁剪或使用反射边界条件数据缺失采用时空插值或邻近区域替代尺度效应进行多尺度分析验证结论稳定性权重敏感度进行敏感性分析测试不同权重方案的影响在最近的一个城市群生态评估项目中我们发现使用5×5km的移动窗口能最佳平衡局部细节与区域趋势。当LUI超过65时生物多样性指数通常会出现显著下降这一阈值可作为生态预警的重要参考。