用GEE和Landsat 8数据，5步搞定城市生态健康“体检报告”（附完整代码）

城市生态健康体检实战用GEE和Landsat 8生成可视化评估报告城市规划师和环保工作者常常需要快速评估城市生态状况但传统方法耗时费力。Google Earth EngineGEE平台结合Landsat 8数据为我们提供了一种高效解决方案。本文将手把手教你如何通过5个步骤完成城市生态健康体检并生成专业级的可视化报告。1. 准备工作与环境搭建在开始之前我们需要准备好GEE账号和基础开发环境。GEE是一个强大的地理空间分析平台集成了PB级卫星影像数据特别适合处理大范围、长时间的遥感分析任务。基础环境配置步骤访问Google Earth Engine官网注册开发者账号登录后进入代码编辑器界面新建一个空白脚本文件提示GEE目前对学术和非商业用途免费开放但商业使用需要申请企业版权限对于初次接触GEE的用户建议先熟悉平台的基本操作左侧面板脚本管理、资产管理和搜索功能中间区域代码编辑窗口右侧面板地图显示和图层控制上方工具栏运行、保存和分享按钮// 测试GEE环境是否正常工作 print(Hello, Google Earth Engine!); Map.setOptions(HYBRID); // 设置底图为卫星混合地图2. 定义研究区域与数据准备确定研究区域是整个分析的第一步。在GEE中我们可以通过多种方式定义研究区边界最常用的是绘制多边形或导入已有的Shapefile文件。研究区定义技巧城市尺度分析建议范围在50×50公里以内区域尺度分析可适当扩大范围但需考虑计算资源时间范围选择Landsat 8数据从2013年至今建议选择完整年度的数据// 示例定义青岛市黄岛区为研究区域 var roi ee.Geometry.Polygon( [[[120.121, 35.975], [120.121, 35.886], [120.257, 35.886], [120.257, 35.975]]], null, false ); Map.centerObject(roi, 12); // 将地图中心定位到研究区缩放级别12Landsat 8数据筛选与预处理Landsat 8数据在GEE中有多个版本我们选择地表反射率产品(TOA)以保证数据质量。云覆盖是影响遥感分析的主要因素需要特别处理。// 定义去云函数 function removeCloud(image) { var qa image.select(BQA); var cloudMask qa.bitwiseAnd(1 4).eq(0); var cloudShadowMask qa.bitwiseAnd(1 8).eq(0); var valid cloudMask.and(cloudShadowMask); return image.updateMask(valid); } // 加载Landsat 8数据并进行预处理 var L8 ee.ImageCollection(LANDSAT/LC08/C01/T1_TOA) .filterBounds(roi) .filterDate(2018-01-01, 2019-12-31) .filterMetadata(CLOUD_COVER, less_than, 50) .map(function(image) { return image.set(year, ee.Image(image).date().get(year)); }) .map(removeCloud);3. 生态指标计算与标准化遥感生态指数(RSEI)由四个核心指标构成分别反映城市生态的不同方面。我们需要依次计算这些指标并进行标准化处理确保它们具有可比性。四大生态指标说明指标名称计算公式生态意义理想值范围绿度(NDVI)(NIR-Red)/(NIRRed)植被覆盖状况0.2-0.8湿度(Wet)波段线性组合地表湿度状况-热度(LST)热红外波段反演地表温度状况-干度(NDBSI)(IBISI)/2建筑和裸地状况-// 计算四大生态指标 var L8imgCol L8imgCol.map(function(img){ // 湿度指标计算 var Wet img.expression( B*(0.1509) G*(0.1973) R*(0.3279) NIR*(0.3406) SWIR1*(-0.7112) SWIR2*(-0.4572),{ B: img.select([B2]), G: img.select([B3]), R: img.select([B4]), NIR: img.select([B5]), SWIR1: img.select([B6]), SWIR2: img.select([B7]) }); img img.addBands(Wet.rename(WET)); // 绿度指标(NDVI)计算 var ndvi img.normalizedDifference([B5, B4]); img img.addBands(ndvi.rename(NDVI)); // 干度指标(NDBSI)计算 var ibi img.expression( (2 * SWIR1 / (SWIR1 NIR) - (NIR / (NIR RED) GREEN / (GREEN SWIR1))) / (2 * SWIR1 / (SWIR1 NIR) (NIR / (NIR RED) GREEN / (GREEN SWIR1))), { SWIR1: img.select(B6), NIR: img.select(B5), RED: img.select(B4), GREEN: img.select(B3) }); var si img.expression( ((SWIR1 RED) - (NIR BLUE)) / ((SWIR1 RED) (NIR BLUE)), { SWIR1: img.select(B6), NIR: img.select(B5), RED: img.select(B4), BLUE: img.select(B2) }); var ndbsi (ibi.add(si)).divide(2); return img.addBands(ndbsi.rename(NDBSI)); }); // 数据标准化处理 var img_normalize function(img){ var minMax img.reduceRegion({ reducer: ee.Reducer.minMax(), geometry: roi, scale: 1000, maxPixels: 10e13, }); var year img.get(year); var normalize ee.ImageCollection.fromImages( img.bandNames().map(function(name){ name ee.String(name); var band img.select(name); return band.unitScale( ee.Number(minMax.get(name.cat(_min))), ee.Number(minMax.get(name.cat(_max))) ); }) ).toBands().rename(img.bandNames()).set(year, year); return normalize; }; var imgNorcol L8imgCol.map(img_normalize);4. 主成分分析与RSEI计算主成分分析(PCA)是RSEI计算的核心步骤它能将多个相关指标转换为少数几个不相关的主成分保留原始数据的大部分信息。PCA在生态评价中的优势消除指标间的冗余信息自动确定各指标权重简化复杂生态系统评估过程结果具有明确的数学意义// 主成分分析实现 var pca function(img){ var bandNames img.bandNames(); var region roi; var year img.get(year); // 数据中心化 var meanDict img.reduceRegion({ reducer: ee.Reducer.mean(), geometry: region, scale: 1000, maxPixels: 10e13 }); var means ee.Image.constant(meanDict.values(bandNames)); var centered img.subtract(means).set(year, year); // 计算协方差矩阵和特征值分解 var arrays centered.toArray(); var covar arrays.reduceRegion({ reducer: ee.Reducer.centeredCovariance(), geometry: region, scale: 1000, bestEffort: true, maxPixels: 10e13 }); var covarArray ee.Array(covar.get(array)); var eigens covarArray.eigen(); // 计算主成分 var eigenVectors eigens.slice(1, 1); var arrayImage arrays.toArray(1); var principalComponents ee.Image(eigenVectors).matrixMultiply(arrayImage); return principalComponents .arrayProject([0]) .arrayFlatten([[PC1, PC2, PC3, PC4]]) .set(year, year); }; var PCA_imgcol imgNorcol.map(pca); Map.addLayer(PCA_imgcol.first().select(PC1), {min: -2, max: 2}, First Principal Component);RSEI计算与可视化第一主成分(PC1)通常包含了原始数据的大部分信息我们可以用它来计算最终的遥感生态指数。// RSEI计算与归一化 var RSEI_imgcol PCA_imgcol.map(function(img){ img img.addBands(ee.Image(1).rename(constant)); var rsei img.expression(constant - pc1, { constant: img.select(constant), pc1: img.select(PC1) }); rsei img_normalize(rsei); return img.addBands(rsei.rename(RSEI)); }); // 可视化参数设置 var visParam { min: 0, max: 1, palette: [ d7191c, fdae61, ffffbf, a6d96a, 1a9641 ] }; Map.addLayer(RSEI_imgcol.first().select(RSEI), visParam, RSEI);5. 结果解读与应用实践获得RSEI结果后如何正确解读并将其应用于实际工作中是关键。RSEI值范围在0-1之间值越大表示生态状况越好。RSEI分级标准参考等级RSEI范围生态状况建议措施优0.8-1.0生态系统健康保护现有生态格局良0.6-0.8生态系统较健康优化生态空间结构中0.4-0.6生态系统一般加强生态修复差0.2-0.4生态系统较差重点治理区域极差0-0.2生态系统恶劣优先整治区域结果应用场景城市规划识别生态敏感区优化城市空间布局环境评估监测生态变化趋势评估政策效果生态修复确定优先治理区域指导修复工程学术研究支持城市生态系统的长期监测研究// 统计不同生态等级的面积比例 var areaStats function(image) { var area ee.Image.pixelArea().addBands(image); var stats area.reduceRegion({ reducer: ee.Reducer.sum().group({ groupField: 1, groupName: level, }), geometry: roi, scale: 30, maxPixels: 1e13 }); return stats; }; // 定义分级函数 var classifyRSEI function(image) { return image.select(RSEI).gt(0.8).multiply(1) .add(image.select(RSEI).gt(0.6).multiply(1)) .add(image.select(RSEI).gt(0.4).multiply(1)) .add(image.select(RSEI).gt(0.2).multiply(1)) .rename(RSEI_class); }; var classified RSEI_imgcol.map(classifyRSEI); var stats classified.map(areaStats); print(Area statistics by RSEI level:, stats);报告生成技巧使用GEE的导出功能将结果保存为GeoTIFF或CSV格式在QGIS或ArcGIS中进行进一步制图结合行政区划数据生成分区统计报告制作时间序列动画展示生态变化过程