从森林到城市：用GEE分析土地利用变化的保姆级教程（含避坑指南）

从森林到城市用GEE分析土地利用变化的保姆级教程含避坑指南当一片森林在五年内消失取而代之的是钢筋混凝土的楼群这种变化背后隐藏着怎样的数据故事Google Earth EngineGEE提供的土地利用转移矩阵分析工具就像一台时间机器能让我们精确量化每一寸土地的命运转折。本教程将手把手带您完成从数据准备到可视化输出的全流程特别针对生态保护和城市规划领域的常见需求设计。1. 数据准备选择适合的分析武器在开始分析前选择合适的数据集至关重要。GEE平台上有多个主流土地利用数据集各有特点数据集名称空间分辨率时间跨度覆盖类型适用场景CGLS-LC100100米2015-至今21类高精度变化检测MODIS Land Cover500米2001-202117类大范围趋势分析ESA WorldCover10米2020-202111类城市精细规划USGS NLCD30米2001-201916类北美区域研究以CGLS-LC100为例加载数据// 定义研究区域以昆明为例 var roi ee.Geometry.Rectangle([102.3, 24.8, 103.0, 25.2]); // 加载2018和2022年数据 var landcover2018 ee.Image(COPERNICUS/Landcover/100m/Proba-V-C3/Global/2018) .select(discrete_classification) .clip(roi); var landcover2022 ee.Image(COPERNICUS/Landcover/100m/Proba-V-C3/Global/2022) .select(discrete_classification) .clip(roi);注意不同数据集使用的分类体系可能不同跨数据集比较时需要重分类统一标准2. 转移矩阵计算核心步骤解析2.1 土地类型编码处理GEE中的分类数据通常用数字编码表示不同土地类型。我们需要先了解这些编码的含义// 获取CGLS-LC100分类体系 var classNames [ 未知, 森林, 灌木, 草地, 耕地, 建筑, 裸地, 雪冰, 水体, 湿地 ]; // 创建分类字典 var classDict ee.Dictionary.fromLists( ee.List.sequence(0, 9), ee.List(classNames) );2.2 构建转移组合图像转移矩阵的核心是将两个时期的分类图像合并为一个组合图像// 将2018年编码放大1000倍后与2022年编码相加 var transferImage landcover2018 .multiply(1000) .add(landcover2022); // 计算每个像素的面积平方米 var areaImage ee.Image.pixelArea();2.3 面积统计与矩阵构建使用分组统计方法计算各类转移的面积// 统计各转移组合的面积 var areaStats areaImage.addBands(transferImage) .reduceRegion({ reducer: ee.Reducer.sum().group({ groupField: 1, groupName: transfer_code }), geometry: roi, scale: 100, maxPixels: 1e13 }); // 将统计结果转换为特征集合 var transferFeatures ee.FeatureCollection( ee.List(areaStats.get(groups)).map(function(group) { group ee.Dictionary(group); var code ee.Number(group.get(transfer_code)); return ee.Feature(null, { from: code.divide(1000).int(), to: code.mod(1000), area_km2: ee.Number(group.get(sum)).divide(1e6) }); }) );3. 可视化技巧与常见问题解决3.1 动态转移矩阵可视化将统计结果转换为交互式热力图// 生成转移矩阵表格数据 var matrixData ee.FeatureCollection( classNames.map(function(fromName, fromIdx) { return ee.Feature(null, ee.Dictionary.fromLists( classNames, classNames.map(function(toName, toIdx) { var match transferFeatures .filter(ee.Filter.eq(from, fromIdx)) .filter(ee.Filter.eq(to, toIdx)); return ee.Number( ee.Algorithms.If( match.size(), match.first().get(area_km2), 0 ) ); }) )); }) ); // 导出到Google Drive Export.table.toDrive({ collection: matrixData, description: LandTransferMatrix, fileFormat: CSV });3.2 典型错误排查指南问题1统计结果面积异常大原因未设置合适的scale参数解决明确指定与数据分辨率匹配的scale值问题2部分转移组合缺失检查项确认两个时期使用相同的分类体系验证研究区域是否包含该土地类型检查filter条件是否准确问题3导出数据乱码解决方案使用CSV而非Excel格式在本地用文本编辑器检查原始数据确认字符编码为UTF-84. 进阶应用场景案例4.1 森林退化预警系统通过时序分析建立早期预警模型// 计算森林转其他类型的年变化率 var forestLoss transferFeatures .filter(ee.Filter.eq(from, 1)) // 1森林 .filter(ee.Filter.neq(to, 1)) .aggregate_sum(area_km2);4.2 城市扩张模拟预测基于历史趋势预测未来城市用地// 计算城市扩张速率 var urbanGrowthRate ee.Number( transferFeatures .filter(ee.Filter.eq(to, 5)) // 5城市 .aggregate_sum(area_km2) ).divide(4); // 4年变化量 // 预测5年后城市面积 var currentUrban landcover2022.eq(5).multiply(ee.Image.pixelArea()) .reduceRegion({reducer: ee.Reducer.sum(), geometry: roi, scale: 100}) .get(constant); var predictedUrban ee.Number(currentUrban).divide(1e6) .add(urbanGrowthRate.multiply(5));4.3 生态保护成效评估建立保护效果量化指标// 计算保护区内自然植被净变化 var protectedArea ee.FeatureCollection(path/to/protected_areas) .filterBounds(roi); var natureChange transferFeatures .filterBounds(protectedArea) .map(function(f) { var from ee.Number(f.get(from)); var to ee.Number(f.get(to)); var isDegradation from.lte(4) to.gt(4); // 自然→人工 var isRestoration from.gt(4) to.lte(4); // 人工→自然 return f.set({ change_type: ee.Algorithms.If(isDegradation, -1, ee.Algorithms.If(isRestoration, 1, 0)) }); });