马山半岛生态环境分析 | 马山半岛生态与旅游分析

第 1 步：获取“矢量边界”数据

你需要一个马山半岛的精确轮廓。有了它，我们才能把庞大的地球影像裁剪成只保留研究区的范围。

我们可以使用 Overpass Turbo 网站来收集相关的边界数据。

javascript

[out:json][timeout:30];

// nwr = node + way + relation (所有类型)
nwr({{bbox}});

out geom;

这段代码使用的是 Overpass QL (Overpass Query Language)，这是专门用来查询 OpenStreetMap 数据的查询语言。但是要注意，这段代码会加载屏幕地图上显示的所有数据，数据量可能会很大。因此，建议先放大到你想查询的区域，然后再点击运行。

研究区域的北部边界太湖边界

或者，如果你知道研究区域的边界编码，可以直接基于 ID 进行精确查询：

javascript

[out:json][timeout:30];

(
  way(4559718);    // 道路网界线
  way(71306905);   // “无锡”河界线
);

out geom;

研究区边界导出 GeoJSON

从 Overpass Turbo 导出为标准格式：

看网页最上方的菜单栏，点击 “导出” (Export) 按钮。
在弹出的菜单中，找到 “数据” (Data) 这一栏。
点击下载 GeoJSON 格式（这是目前最通用、且 ArcGIS 完美支持的现代地理数据格式）。

浏览器会下载一个文件，默认名称为 /Users/ousin/Downloads/export.geojson。

:::{.fold} 新建项目的文件夹结构是怎样的？

ArcGIS 中新建项目的文件夹

当你新建一个名为 MyProject 的项目时，ArcGIS 只是在指定路径下创建了一个 MyProject 文件夹，并在其中打包了 4 个核心项目：

Index 文件夹（用于加速本地搜索索引）
MyProject.aprx（项目排版文件，记录你打开了哪些地图）
MyProject.atbx（该项目的默认工具箱）
MyProject.gdb（该项目的默认地理数据库）

:::

:::{.fold} 如何在 ArcGIS 中删除项目？

删除项目

这是一个非常好的问题！许多初学者在这个界面都会感到困惑：为什么右键菜单里只有“从列表中移除 (Remove Project From List)”，而没有“删除 (Delete)”？

因为 ArcGIS Pro 的逻辑是：项目本质上就是硬盘上的一个文件夹。软件内部不提供一键删除功能，是为了防止你不小心误删了里面重要的地理数据库（.gdb）。

在 ArcGIS 中删除项目的步骤：

从欢迎界面清理：右键点击你想删除的项目（如图中的 MyProject）并选择“从列表中移除”。这只会删除快捷方式。
定位文件夹：如果你还没有点击移除，可以点击“在文件资源管理器中显示 (Show In File Explorer)”来打开项目所在的文件夹。
[关闭 ArcGIS Pro：非常重要！]{style="background-color: red; color: white;"} 如果你在打开软件的情况下尝试删除文件夹，Windows 会报错提示“文件被占用”，或者导致 .gdb 文件被锁定并损坏。
删除硬盘文件：关闭软件后，将该项目对应的整个文件夹（例如 MyProject1）删除到回收站。

💡 极客删除法：因为你的文件都共享在 Mac 上，你甚至不需要打开 Windows 的文件资源管理器。你可以直接在 Mac 终端中输入命令，瞬间将它们彻底删除（例如，我们要删除 MyProject 和 MyProject1）：

bash

rm -rf ~/Documents/ArcGIS/Projects/MyProject
rm -rf ~/Documents/ArcGIS/Projects/MyProject1

(警告：rm -rf 是不可逆的删除命令，千万小心不要把 mashan 给删了！)

:::

:::{.fold} QGIS 要在 APP 中删除项目文件夹，而 ArcGIS 可以直接删除文件夹。

这是一个非常敏锐的对比！你对 QGIS 的记忆很准确，但这主要是因为 QGIS 和 ArcGIS Pro 在“项目管理”的底层哲学上完全不同。我们可以简单对比一下：

QGIS：散养模式（指针链接）

QGIS 没有强制的“项目文件夹”概念。当你保存一个 QGIS 项目时，它只是生成了一个极其小巧的 .qgz（或 .qgs）单文件。
这个文件里面只保存了路径指针和样式设置。比如它记录了：“马山边界在下载文件夹里，颜色是红色”。
为什么不能直接删文件夹：因为你的数据（Shapefile、TIFF）通常是散落在电脑各个地方的。如果你只删了 .qgz 所在的文件夹，你的数据并没有被删除；反过来，如果你不小心清理了其他文件夹（比如删了下载目录），你的 QGIS 项目再打开时就会全屏飘红，提示“图层丢失 (Bad Layers)”。

ArcGIS Pro：全家桶模式（高内聚）

ArcGIS Pro 强制推行“工程 (Project)”思维。只要你新建项目，它就强制在硬盘上圈出一块地（建一个文件夹）。
它强制为你在这个文件夹里建好专属的数据库（.gdb）。它的理念是：把你所有处理过的数据、中间产物、出图排版，全部“收容”在这个文件夹里。
为什么可以直接删：正因为这种“打包”机制，只要你遵循它的规范（把数据都存进 .gdb），当你想要销毁这个项目时，只要把外层的“全家桶”文件夹扔进垃圾桶，就干干净净，绝不拖泥带水。

这也是为什么在企业级应用和大型科研项目中，很多人觉得 ArcGIS Pro 的文件管理更让人安心，不容易出现把工程拷给同事后，对方打开发现“图层全部丢失”的尴尬局面（除非你还是坚持把数据存在系统桌面而不是 .gdb 里）。

:::

使用命令删除项目只剩下我们的 mashan 项目文件

先在你的 ArcGIS 工程目录下，建立一个专门用来存放原始数据的结构（包含矢量和栅格分类）：

bash

mkdir -p ~/Documents/ArcGIS/Projects/mashan/01_Raw_Data/Vector
mkdir -p ~/Documents/ArcGIS/Projects/mashan/01_Raw_Data/Raster

使用命令文件夹效果

把刚才下载的 export.geojson 移动到新建的 Vector 文件夹中，并顺手改个能看懂的名字（Mashan_Boundary）：

bash

mv ~/Downloads/export.geojson ~/Documents/ArcGIS/Projects/mashan/01_Raw_Data/Vector/Mashan_Boundary.geojson

查看确认：

bash

ls -l ~/Documents/ArcGIS/Projects/mashan/01_Raw_Data/Vector/

text

📁 /Users/ousin/Documents/ArcGIS/Projects/mashan/
 ├── 📁 01_Raw_Data           <-- [新建的] 存放所有未经处理的原始下载数据
 │    ├── 📁 Raster           <-- 之后下载的卫星影像(.tif)放这里
 │    └── 📁 Vector           <-- 矢量边界放这里
 │         └── Mashan_Boundary.geojson  <-- [我刚把你的文件移动过来并改名了！]
 ├── mashan.aprx              <-- 你的工程文件
 └── mashan.gdb               <-- 你的数据库（处理过后的数据、最终图层都会存进这里）

我将外部下载的原始数据（GeoJSON、TIFF）通常放在单独的 Raw_Data 普通文件夹里保管。

而经过转换、裁剪、真正用于出图的内部数据（比如你马上要生成的 Mashan_ROI 多边形面），都会全部收纳进那个 .gdb（地理数据库）文件中统一管理。

第 2 步：将 GeoJSON 导入 ArcGIS 并生成多边形 (ROI)

在准备好了矢量线段后，我们需要在 ArcGIS Pro 中将其转换为实心的多边形（Polygon），作为后续裁剪卫星影像的“模具”。

2.1 JSON To Features (导入线条)

首先，使用 JSON To Features 工具将 GeoJSON 导入到我们的 .gdb 数据库中。

JSON To Features 工具

Input JSON (输入 JSON)：点击旁边的文件夹图标 📂，找到存好的 Mashan_Boundary.geojson。
Output Feature Class (输出要素类)：默认保存在 mashan.gdb 里，我们将其命名为 Mashan_Lines_Raw（或者如图所示的其他名字）。
Geometry Type (几何类型)：[非常重要]{style="color: red;"} 由于我们在 Overpass 中提取的是道路和河流（线段），这里必须选择 Polyline（折线）或保持默认。

⚠️ 避坑指南：如果你强行将 Geometry Type 选为 Polygon，由于输入数据本身没有面要素，ArcGIS 会弹出一个黄色的警告 WARNING 000117: Warning empty output generated. 并且生成一个没有任何数据的空图层。

报错示例

正确选择 Polyline 后，运行成功，地图上会加载出半岛的线框：

成功显示线框

2.2 Feature To Polygon (线转面，清理悬挂线段)

从上面的截图中可以看到，提取的路网向北部延伸出了两条多余的“天线”（悬挂线段）。我们不需要手动去删除它们，ArcGIS 的拓扑算法可以自动搞定。

使用 Feature To Polygon (要素转面) 工具：

Feature To Polygon 配置

Input Features：选择刚才生成的线图层。
Output Feature Class：命名为 Mashan_ROI (Region of Interest)。

工具会自动寻找线条交叉形成的闭合包围圈，生成实心的半岛多边形，同时自动抛弃那些没有闭合的无用线段！

至此，我们的矢量边界预处理阶段全部完成。有了这个完美的粉色多边形，我们接下来就可以去天上“扒”卫星图了。

第 3 步：使用 Google Earth Engine (GEE) 获取无云卫星影像

为了进行生态分析（如计算 NDVI），我们需要一张覆盖马山半岛的高清多光谱卫星图。由于从传统的 USGS 官网下载原始影像体积巨大且需要繁琐的大气校正，我们选择使用极客方案：Google Earth Engine (GEE)。

GEE 的超级计算机可以在云端直接过滤云层、拼接影像，并裁剪出刚好覆盖研究区的一小块，直接导出为轻量级的 .tif 格式。

3.1 在 GEE 中框选研究区

登录 GEE Code Editor。
在地图上方的搜索框输入 Wuxi 定位到无锡，并放大找到太湖边的马山半岛。

GEE Code Editor 界面定位无锡马山

点击地图左上角的正方形（矩形工具） ⬛。

矩形工具

在地图上画一个刚好能框住整个马山半岛的矩形。此时，代码顶部会自动生成一行 var geometry = ...。

拖动框选

3.2 运行 GEE 代码截取 Sentinel-2 影像

我们选择欧洲空间局的 Sentinel-2 卫星，它的分辨率高达 10 米，非常适合区县级的精细生态分析。

在代码框第一行粘贴以下代码：

javascript

// 1. 设置时间范围 (选取植被最茂盛的夏季，云量少的时间段)
var startDate = '2023-06-01';
var endDate = '2023-10-31';

// 2. 调用 Sentinel-2 表面反射率数据集 (已做过大气校正)
var s2 = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
  .filterBounds(geometry)
  .filterDate(startDate, endDate)
  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20));

// 3. 叠加求中位数，生成无云合成图
var composite = s2.median().clip(geometry);

// 4. 在地图上预览真彩色 (RGB) 效果
var visParams = {bands: ['B4', 'B3', 'B2'], min: 0, max: 3000};
Map.addLayer(composite, visParams, 'Mashan Sentinel-2 Preview');

// 5. 核心命令：将提取的 4 个关键波段导出到 Google 云盘
Export.image.toDrive({
  image: composite.select(['B2', 'B3', 'B4', 'B8']), 
  description: 'Mashan_Sentinel2_10m',
  folder: 'GEE_Exports',
  scale: 10,
  region: geometry,
  maxPixels: 1e10
});

:::{.fold} 我们用到了 Sentinel-2，使用了它的哪些波段？

好的，让我帮你梳理清楚。在你笔记的第 3 步中，我们在 GEE 里导出卫星影像时，代码里有这一行关键语句：

javascript

image: composite.select(['B2', 'B3', 'B4', 'B8']),

我们从 Sentinel-2 卫星一共选了 4 个波段导出，它们分别是：

导出时的波段名	导入 ArcGIS 后的名称	光的类型	波长	人眼可见？	用途
B2	Band_1	🔵 蓝光	~490 nm	✅	水体识别、大气校正
B3	Band_2	🟢 绿光	~560 nm	✅	植被反射峰
B4	Band_3	🔴 红光	~665 nm	✅	植物大量吸收（光合作用）
B8	Band_4	🟤 近红外 (NIR)	~842 nm	❌	植物强烈反射

其中前三个（蓝、绿、红）组合起来就是你在地图上看到的真彩色卫星图（和人眼看到的一样）。

第四个（近红外 B8）是人眼看不见的，但它是计算 NDVI 的灵魂波段。这也就是为什么在第 5 步计算 NDVI 时，我们设置 NIR Band = 4（对应 B8）、Red Band = 3（对应 B4）。

顺带提一个有趣的细节

Sentinel-2 卫星其实一共有 13 个波段，我们只挑了 4 个最核心的。如果当时多导出几个短波红外波段（B11、B12），我们就可以直接在 ArcGIS 里计算 MNDWI（水体指数）和 NDBI（建筑指数）——这恰恰就是你之前 GEE 脚本 lucc_classification.py 中用来做全自动 LUCC 分类的那三个指数！

:::

点击顶部的 Run 按钮执行代码。随后，右侧的 Tasks 面板会亮起黄灯，提示有导出任务。点击 Tasks 里的 Run，确认参数无误后，再次点击蓝色 RUN 开始导出。

导出设置任务运行中

等待约 1~3 分钟，当齿轮变成绿色对勾 ✅，说明影像已经成功导出到了 Google Drive！

3.3 将影像移入本地 ArcGIS 目录

访问 Google Drive，在 GEE_Exports 文件夹中找到我们导出的 Mashan_Sentinel2_10m.tif。由于 GEE 在云端已经做了精准裁剪和波段筛选，这个原本应该有几个 GB 的原始文件，被精简到了只有 8MB 左右。

Google Drive 中的成果

下载文件后，我们在 Mac 终端中使用 mv 命令将其移动到我们事先规划好的 ArcGIS 工程 Raster 目录下：

bash

mv '/Users/ousin/Library/CloudStorage/GoogleDrive-susangharabaghipii9@gmail.com/My Drive/GEE_Exports/Mashan_Sentinel2_10m.tif' ~/Documents/ArcGIS/Projects/mashan/01_Raw_Data/Raster/

执行 mv 命令

现在，你的 ArcGIS 项目原材料已经全部归位！矢量边界和栅格影像已准备就绪，可以开始接下来的 ArcGIS 掩膜提取与分析了。

完美的文件夹结构

第 4 步：在 ArcGIS 中加载影像并按掩膜提取

影像准备好后，我们回到 Windows 虚拟机中的 ArcGIS Pro 进行最终的“修边”工作。

4.1 避坑：跨系统共享路径 Bug 与火星坐标系偏移

在这个过程中，你很可能会遇到两个非常经典的 GIS 跨平台实战“暗坑”：

踩坑 1：虚拟机网络路径导致 Invalid Path 在 ArcGIS Pro 中，如果使用顶部菜单的 Add Data -> Data From Path 尝试加载存放在 Mac 共享文件夹（如 \\Mac\Home\Documents\...）中的 .tif 文件，GDAL 渲染引擎经常会报错 unsupported data type 或 Invalid Path。

尝试使用 Add Data 从路径加载 Invalid Path 报错

解决方案：直接避开 Add Data 菜单。使用右侧的 Catalog (目录) 面板找到文件，或者直接在 Windows 资源管理器中按住鼠标左键拖拽 .tif 文件到地图画布中，即可完美绕过路径解析 Bug！

直接拖拽导入成功

踩坑 2：中国大陆特有的“火星坐标系 (GCJ-02)”偏移 当你把第一次导出的卫星图拖进 ArcGIS 并盖上粉色多边形时，你可能会发现：**为什么半岛最西边有一小块凸出了卫星图的范围？**明明在 GEE 里画框的时候把整个半岛都包住了呀？

多边形左侧超出卫星图边界在 GEE 街道地图中看似已经包住

原理解析：这是因为 GEE 默认的街道地图/地形图在中国大陆区域强制使用了加密偏移的火星坐标系 (GCJ-02)。但我们通过代码提取的 Sentinel-2 卫星影像使用的是绝对真实的物理坐标 (WGS-84)。
因此，你看着偏移后的“假地图”画的框，套在“真卫星图”上时，就会出现几百米的错位！

街道图与实际卫星图的偏移对比

终极解法：在 GEE 画框之前，必须将右上角的底图切换为 “卫星检视 (Satellite)”。Google 的卫星底图是没有偏移的 WGS-84 坐标，与 Sentinel-2 完全吻合。
在真实的卫星底图上，重新大方地画一个把半岛包裹得严严实实的矩形，重新运行代码导出 _V2 版本。

切换为卫星视图并重新画框导出 V2 版本

4.2 加载 V2 影像与生成附属文件 (.aux)

将重新下载的 _V2.tif 拖进 ArcGIS Pro 后，软件会弹出一个提示框：Build Pyramids and Calculate Statistics (构建金字塔和计算统计数据)。

构建金字塔与统计数据弹窗

为什么会弹这个框？ 这说明你的原始遥感图像是一张“原片”，ArcGIS 需要扫描图像内所有像素的最亮/最暗值（Statistics 统计数据），从而决定如何正确地给你呈现色彩。同时，因为图像较大，它建议为你生成不同缩放比例下的模糊预览图（Pyramids 金字塔），以保证你在缩放地图时不卡顿。
直接保持默认勾选，点击 OK。

点击确定后，你如果去 Windows 资源管理器里查看，会发现 .tif 原文件旁边多出了同名的 .aux 或 .aux.xml 文件。这就是 ArcGIS 的 附属文件 (Sidecar files)，它把刚才计算的金字塔和统计数据记在了里面，这样下次再加载这张图时就能瞬间读取，不会再弹框了。

自动生成的 .aux 附属文件

4.3 执行 Extract by Mask（按掩膜提取）

现在，卫星原片已经完美覆盖，并且没有了边界残缺的问题。为了不让周围大面积的太湖水体干扰我们后续的植被指数 (NDVI) 平均值计算，我们需要将多余的水面裁掉。

在顶部菜单栏点击 Analysis (分析) -> Tools (工具)，打开右侧的 Geoprocessing 搜索面板。
搜索并打开工具：Extract by Mask (按掩膜提取)。
Input Raster (输入栅格)：选择我们的卫星原片 Mashan_Sentinel2_10m_V2.tif。
Input raster or feature mask data (输入掩膜数据)：选择我们的粉色多边形“刻刀” Mashan_ROI。
Output Raster (输出栅格)：命名为 Mashan_Image_Masked（注意它会被默认存入 mashan.gdb 地理数据库中）。

Extract by Mask

点击底部的 Run。

运行完成后，在左侧的 Contents 面板中隐藏掉原来的 _V2.tif 图层，你就会得到一张边缘像刀切一样整齐，完美贴合马山半岛轮廓的卫星图片！

马山半岛卫星图

第 5 步：提取大自然的“绿色指纹” (计算 NDVI)

在生态地理学中，NDVI（归一化植被指数） 是衡量植被健康度最核心的指标。它的基本公式是： $NDVI = \frac{NIR - Red}{NIR + Red}$ 。由于健康的绿色植物会大量吸收红光并强烈反射近红外光，因此这两者的差值越大，代表植被越茂密。

:::{.fold} NDVI 是什么？

NDVI（Normalized Difference Vegetation Index，归一化植被指数） 本质上是一个介于 -1 到 1 之间的连续数值，用来回答一个问题："这块地上的植物有多绿、多健康？"

健康的绿色植物有一个独特的物理特性：它会疯狂[吸收红光（用于光合作用），同时强烈反射近红外光]{style="background-color: red; color: white;"}。卫星从太空拍下来后，通过对比这两个波段的差异，就能精确判断地面的"绿"的程度：

NDVI 范围	含义	对应地物
0.8 ~ 1.0	极高植被覆盖	茂密森林
0.3 ~ 0.5	中等植被覆盖	农田、草地
0 ~ 0.1	几乎无植被	水面、建筑、裸地
< 0	非植被	云、雪、水体

在论文中的作用：NDVI 相当于一张"植被体检报告"。通过对比不同年份的 NDVI 图，可以直观地看到马山半岛的绿色在 20 年间是变多了还是变少了，从而量化生态环境的演变趋势。 :::

:::{.fold} 红光是什么？近红外光是什么？

好问题！这涉及遥感的基础物理知识。

光的本质是电磁波

我们看到的"光"只是电磁波谱中极其狭窄的一小段。按照波长从短到长排列：

text

紫外线 → [紫  蓝  绿  黄  橙  红] → 近红外 → 短波红外 → 热红外 → 微波
        ← 人眼可见光（380~700nm）→

名称	波长范围	人眼能看到吗？	遥感中的用途
红光 (Red)	~620–700 nm	✅ 能看到（就是红色）	植物会大量吸收它来做光合作用
近红外 (NIR)	~700–1300 nm	❌ 看不到	植物叶片会强烈反射它（像镜子一样）
短波红外 (SWIR)	~1300–2500 nm	❌ 看不到	对水分和土壤极敏感
热红外 (TIR)	~8000–14000 nm	❌ 看不到	测量地表温度（热力学）

所以回答你的三个问题：

红光：就是你肉眼能看到的红色光，波长约 620~700 nm。Sentinel-2 的 B4 波段就是专门拍它的。
近红外光 (NIR)：紧挨着红光"右边"（波长稍长），人眼已经看不见了，但卫星的传感器可以"看到"。Sentinel-2 的 B8 波段拍的就是它。
红外光：是的，"红外"是一个很大的家族，按波长从短到长分为近红外 → 短波红外 → 中红外 → 热红外。我们 NDVI 用的只是最"近"的那一段（近红外）。

为什么植物对这两种光的反应如此不同？

这背后是植物细胞的物理结构决定的：

叶绿素会贪婪地吸收红光和蓝光（用于光合作用），所以植物在红光波段反射极低
叶肉细胞的内部结构（海绵组织）会像一面镜子一样把近红外光强烈反射回去

正因如此，NDVI 公式 $NDVI = \frac{NIR - Red}{NIR + Red}$ 的精妙之处在于：

健康植物 → NIR 极高、Red 极低 → 差值大 → NDVI 接近 1
水面/建筑 → NIR 和 Red 都很低或差不多 → 差值小 → NDVI 接近 0

:::

5.1 认清波段序号并一键计算

我们在 GEE 中导出 Sentinel-2 影像时，波段的对应关系如下：

Band_1 = B2 (蓝光)
Band_2 = B3 (绿光)
Band_3 = B4 (红光 Red)
Band_4 = B8 (近红外 NIR)

在 ArcGIS Pro 中，只需点击顶部 Imagery 选项卡，在 Indices 菜单中选择 NDVI；在弹出的窗口中，手动将 NIR Band 设为 4，将 Red Band 设为 3：

正确的

点击确定后，我们得到了一张极其精准的 NDVI 黑白图像（白色代表高值森林，黑灰代表低值水体和建筑）：

未配色的NDVI图像

5.2 避坑：临时图层断连与“转正大法”

注意，通过 Indices 工具一键生成的 NDVI 默认是一个存在内存里的临时图层。当你直接给它赋予复杂的渲染色带时，图层很容易因为无法及时计算后台的统计数据（Statistics）而崩溃断连，表现为图层前出现红色感叹号 ❗️，右侧面板弹出 Failed to set statistics 警告。

完美解法（转正大法）：遇到临时图层罢工，直接删掉重做。趁新图层刚生成还没报错，立刻右键图层 -> Data -> Export Raster (导出栅格)。将输出格式设为 TIFF，Pixel Type 保持 32 Bit float (32位浮点型) 以保留小数精度，将其导出为物理硬盘上的永久文件。

Export Raster 导出成功

5.3 生态学配色与掩膜去黑底

拥有永久的 NDVI 文件后，我们在 Symbology (符号系统) 中给它换上学术级配色，并去除由于裁切导致的冗余黑底。

Symbology(符号系统) Condition Number

设置色带：在 Color scheme 搜索栏输入 Condition Number，选中这条从红到绿的连续渐变色带。记得

💡 混用两种卫星在学术上是完全可以的，大量发表在《Remote Sensing of Environment》、《ISPRS Journal》等顶刊的论文都采用 Landsat + Sentinel-2 联合方案。只需在方法论章节中注明数据来源和分辨率差异即可。

但经过进一步调研，我们决定采用更稳妥的方案：全部 6 个年份统一使用 Landsat 30m 数据（详见下方"多年份数据方案"一节）。

:::{.fold} 各个卫星的分辨率 10 米和 30 米是什么意思？

分辨率 = 每个像素代表地面上多大一块地

想象你用马赛克来拼一幅画：

10 米分辨率（Sentinel-2）：每一块马赛克小方块代表地面上 10m × 10m = 100 m² 的区域。相当于用小颗粒拼图，细节丰富。
30 米分辨率（Landsat）：每一块马赛克代表 30m × 30m = 900 m² 的区域。相当于用大颗粒拼图，比较粗糙。

也就是说，同样覆盖 1 km² 的区域：

	Sentinel-2 (10m)	Landsat (30m)
每个像素覆盖面积	100 m²	900 m²
1 km² 需要多少像素	10,000 个	约 1,111 个
精细度	能分辨一栋楼	只能分辨一个街区

实际效果差异

举个例子，马山半岛上有一条 20 米宽的公路：

10 米分辨率：公路占了 2 个像素的宽度，卫星能"看见"这条路 ✅
30 米分辨率：公路不到 1 个像素宽，它会和旁边的农田、树木"混"在同一个像素里，被平均掉 ❌

对论文的影响

本论文分析的是半岛尺度的大面积土地变化（森林 vs 建设用地），不需要精确到单栋建筑。30 米分辨率对于这个研究尺度完全够用，这也是绝大多数 LUCC 论文都用 Landsat 30m 数据的原因。

:::

5.1 MNDWI / NDBI 可视化验证

将计算好的 MNDWI 和 NDBI 进行配色，用于视觉验证计算结果的正确性。

MNDWI 水体识别效果：使用 Classify 分类模式，设 2 个类，断点设为 0（MNDWI > 0 标蓝色为水体，≤ 0 标灰色为陆地）：

MNDWI 配色操作 Mask 去底

MNDWI 水体识别结果

可以看到研究区内的池塘、水库等小型水体被成功标蓝，说明 MNDWI 计算正确。

:::{.fold} 阈值法（如 MNDWI > 0 就判定为水）可靠吗？

不完全可靠。 MNDWI > 0 只是一个粗略的分界线：

MNDWI 值	通常对应	但也可能是...
> 0.3	✅ 几乎确定是水	极少误判
0 ~ 0.3	⚠️ 可能是水	也可能是潮湿泥土、阴影、深色屋顶
< 0	通常不是水	植被、建筑、农田

这就是为什么我们不能仅靠阈值法做最终分类。随机森林分类会同时综合 NDVI、MNDWI、NDBI + 6 个原始波段，从多个维度判断每个像素的真实类型。例如：

某像素	MNDWI	NDBI	NDVI	阈值法判断	随机森林判断
池塘	0.5	-0.3	-0.1	✅ 水	✅ 水
深色屋顶	0.1	0.4	-0.2	❌ 误判为水	✅ 建筑
湿泥地	0.05	0.1	0.2	❌ 误判为水	✅ 农田

:::

5.2 指数的定位：中间产品，而非最终结论

⚠️ 重要理解：NDVI、MNDWI、NDBI 这三个指数不需要提前分类。它们保持原始的连续数值（-1 到 +1），作为随机森林分类的输入特征使用。

随机森林的输入（每个像素有 9 个数值）：

  Band_1 ~ Band_6    ← 6 个原始波段（光谱信息）
  NDVI               ← 连续数值（植被信号）
  MNDWI              ← 连续数值（水体信号）
  NDBI               ← 连续数值（建筑信号）

随机森林的输出：→ "这个像素是 植被 / 水体 / 建设用地 / 农田"

上面做的配色和二分类只是可视化展示，用于验证计算结果、制作论文插图。真正的 LUCC 分类将由随机森林一步完成。

第 6 步：多年份数据方案与技术路线

6.1 为什么需要多年份数据？

目前我们的所有指数都基于 2023 年一个时间点的 Sentinel-2 影像。但论文的核心论点是"旅游开发在 20 年间蚕食了生态"，这需要多年份的对比数据才能得出结论。

6.2 分辨率统一方案

经过调研，学术界对长时序 LUCC 分析的主流做法是：

全部年份统一使用 Landsat 30m 数据集，保持数据源的一致性，避免因卫星更换产生的伪变化。

如果前面用 Landsat (30m)、后面换 Sentinel-2 (10m)，分辨率的突然跳变可能导致误判——比如把"分辨率提高后多识别出的小路"当成"新增的建设用地"。

确定采用的方案：

年份	卫星	分辨率	数据集
2000	Landsat 5	30m	`LANDSAT/LT05/C02/T1_L2`
2005	Landsat 5	30m	`LANDSAT/LT05/C02/T1_L2`
2010	Landsat 5	30m	`LANDSAT/LT05/C02/T1_L2`
2015	Landsat 8	30m	`LANDSAT/LC08/C02/T1_L2`
2020	Landsat 8	30m	`LANDSAT/LC08/C02/T1_L2`
2025	Landsat 9	30m	`LANDSAT/LC09/C02/T1_L2`

我们之前用 Sentinel-2 (10m) 做的 NDVI / MNDWI / NDBI 可以作为 2023 年现状的高精度补充分析，与 Landsat 的长时序分析并行使用。

6.3 分类方法：随机森林监督分类

	阈值法	随机森林监督分类
原理	人工设定规则（如 NDVI > 0.4 → 植被）	标注样本，让机器学习自己找规律
精度	⚠️ 较低，边界模糊	✅ 高，能识别混合像素
工作量	快，5 分钟搞定	每个年份需手动标注 50-100 个训练样本
学术认可度	⚠️ 低，审稿人会质疑	✅ 高，当前顶刊标准方法

6.4 完整技术路线

GEE 导出 Landsat 6 波段（每年一张，统一 30m）
    ↓
ArcGIS 导入 + Extract by Mask（裁剪研究区）
    ↓
计算 NDVI / MNDWI / NDBI（每年各 3 个指数）
    ↓
Classification Wizard → 随机森林监督分类 → 6 张 LUCC 分类图
    ↓
面积统计 + 转移矩阵 → 论文核心数据与结论

6.5 论文中需要的分析与结论

📊 面积统计（基础数据）

用 6 期 LUCC 数据，统计每个年份中各类土地面积：

年份	水体(km²)	植被(km²)	建设用地(km²)	农田(km²)
2000	?	?	?	?
2005	?	?	?	?
2010	?	?	?	?
2015	?	?	?	?
2020	?	?	?	?
2025	?	?	?	?

🔄 土地转移矩阵（核心证据）

对比两个年份（如 2000 → 2025），计算"谁变成了谁"，量化生态被侵占的面积。

🏔️ DEM 地形叠加分析

将建设用地扩张与海拔/坡度进行交叉分析，揭示开发集中在什么地形条件下。的研究，30 米完全满足学术要求。

回到 GEE，在之前的脚本中追加以下代码（复用已有的 geometry 矩形区域）：

javascript

var dem = ee.Image("USGS/SRTMGL1_003").clip(geometry);

Export.image.toDrive({
  image: dem,
  description: 'Mashan_DEM_30m',
  folder: 'GEE_Exports',
  scale: 30,
  region: geometry,
  maxPixels: 1e10
});

⚠️ 注意：必须在之前保存的旧脚本中追加代码（里面已有画好的 geometry），而不是新建脚本。新脚本中 geometry 未定义会报错。同时，确保底图仍为卫星视图，避免 GCJ-02 偏移。

追加代码下载 DEM

导出完成后，从 Google Drive 下载 Mashan_DEM_30m.tif，使用终端命令移入项目目录：

bash

mv '/Users/ousin/Library/CloudStorage/GoogleDrive-susangharabaghipii9@gmail.com/My Drive/GEE_Exports/Mashan_DEM_30m.tif' ~/Documents/ArcGIS/Projects/mashan/01_Raw_Data/Raster/

下载完成移动命令 mv 成功

在 ArcGIS Pro 中，通过 View -> Catalog Pane 打开目录面板，展开 Folders -> mashan -> 01_Raw_Data -> Raster。如果看不到新文件，右键 Raster 文件夹选择 Refresh (刷新)，然后将 Mashan_DEM_30m.tif 拖拽到地图画布中加载。

Catalog 面板 Refresh 后出现 DEM

加载成功后的 DEM 显示马山半岛海拔范围为 -13m ~ 258m，白色为山峰，黑色为低洼水面区域：

DEM 加载成功

5.5 按掩膜提取 DEM 与视觉误区避坑

加载好的 DEM 是一个巨大的矩形，我们需要像之前处理卫星影像一样，用 Mashan_ROI 多边形将其裁剪成半岛的形状。

打开 Analysis -> Tools，搜索并运行 Extract by Mask (按掩膜提取)。
Input raster: Mashan_DEM_30m.tif
Feature mask data: Mashan_ROI
Output raster: Mashan_DEM_Masked

按掩膜提取 DEM

WARNING

经典视觉误区与避坑

运行完毕后，如果你发现在多边形内部大面积呈现黑色，千万不要误以为那是“没裁干净的黑底”！实际上，Extract by Mask 已经完美成功，多边形外部已经完全透明。你看到的黑色，是半岛内部真实存在的低海拔地形（平地与水面）。在 DEM 默认的黑白拉伸配色中，低海拔区域自然会显示为黑色。

绝对不要去 Symbology 的 Mask 标签页里勾选 Display background value = 0 试图“去黑底”！因为马山半岛紧邻太湖，边缘水面和海岸线的真实海拔刚好就是 0 米。如果强制把 0 设为透明，就会在真实数据上“挖孔”，导致半岛边缘破损，出现边缘参差不齐、露出底层彩色影像的错误效果。

正确的DEM提取错误的去底操作导致半岛边缘穿孔破损

正确的做法是：直接保留默认渲染。不需要任何额外去底操作，在 Contents 面板中将其取消勾选隐藏备用即可。

第 7 步：土地利用分类 (LUCC)

:::{.fold} LUCC 是什么？它和 NDVI 有什么关系？

LUCC（Land Use / Land Cover Change，土地利用/覆盖变化） 是一张将地面上每一个像素都打上具体标签的分类地图——"这是水"、"这是森林"、"这是房子"、"这是农田"。

NDVI 只能告诉你"这块地有多绿"，但它分不清"绿色的森林"和"绿色的农田"。LUCC 则更进一步，综合多个光谱指标把每个像素归类为具体的土地类型。

	NDVI	LUCC
是什么	一个连续的数值（0~1）	一个离散的分类标签（1/2/3/4）
回答的问题	"植被有多茂密？"	"这块地到底是什么？"
类比	体温计（告诉你温度是多少）	诊断书（告诉你是感冒还是发烧）
论文中的角色	辅助指标，量化绿色变化趋势	核心数据，支撑土地转型分析

在论文中的作用：LUCC 回答的是一个更深层的问题——"马山半岛的土地用途在 20 年间发生了什么变化？"例如：

2000 年这块地是农田，到 2020 年变成了酒店 → 旅游开发侵占了农业用地
2000 年这块地是森林，到 2020 年变成了停车场 → 生态被破坏了

简而言之：NDVI 是体检数据，LUCC 是诊断结论。 NDVI 可以作为 LUCC 分类的输入特征之一，而 LUCC 的最终分类结果才是做面积统计、画转移矩阵、得出"旅游开发蚕食生态"结论的核心依据。 :::

Classification Wizard 分类类型与方案划分四个类别

训练水域

手动标注太麻烦，遂放弃。

第 8 步：补全波段，计算全套光谱指数

在第 5 步中，我们仅使用 4 个波段（B2、B3、B4、B8）成功计算了 NDVI。但为了进行更全面的生态分析（水体提取、建筑识别等），我们需要补充导出 短波红外（SWIR） 波段，构建完整的 6 波段底图。

8.1 在 GEE 中重新导出 6 波段影像（V3）

搜索原来保存的 GEE 脚本 Mashan_Sentinel2_10m_V2，在其基础上修改导出代码，沿用 V2 中已经校正过的 geometry 划定区域。

run run

javascript

// 1. 设置时间范围 (选取植被最茂盛的夏季，云量少的时间段)
var startDate = '2023-06-01';
var endDate = '2023-10-31';

// 2. 调用 Sentinel-2 表面反射率数据集 (已做过大气校正)
var s2 = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
  .filterBounds(geometry)
  .filterDate(startDate, endDate)
  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20));

// 3. 叠加求中位数，生成无云合成图
var composite = s2.median().clip(geometry);

// 4. 在地图上预览真彩色 (RGB) 效果
var visParams = {bands: ['B4', 'B3', 'B2'], min: 0, max: 3000};
Map.addLayer(composite, visParams, 'Mashan Sentinel-2 Preview');

// ===== 完整版：导出 6 个关键波段（覆盖所有生态指数）=====
Export.image.toDrive({
  image: composite.select(['B2', 'B3', 'B4', 'B8', 'B11', 'B12']),
  description: 'Mashan_Sentinel2_10m_V3',
  folder: 'GEE_Exports',
  scale: 10,
  region: geometry,
  maxPixels: 1e10
});

⚠️ 注意：B11 和 B12 在 Sentinel-2 上原生分辨率是 20 米（不是 10 米），但 GEE 在导出时会自动重采样到指定的 scale: 10，所以文件里所有波段都是统一的 10 米像素。

导出完成后，从 Google Drive 下载并移入本地项目目录：

bash

mv '/Users/ousin/Library/CloudStorage/GoogleDrive-susangharabaghipii9@gmail.com/My Drive/GEE_Exports/Mashan_Sentinel2_10m_V3.tif' ~/Documents/ArcGIS/Projects/mashan/01_Raw_Data/Raster/

mv 命令 mv 成功

8.2 导入 ArcGIS 并按掩膜裁剪

在 Catalog 面板中找到 V3 文件，直接拖拽到地图画布中加载：

直接从目录中拖拽进来导入成功

然后使用 Extract by Mask 将矩形大图裁剪为半岛轮廓：

顶部菜单 Analysis → Tools → 搜索 Extract by Mask
Input Raster：Mashan_Sentinel2_10m_V3.tif
Feature mask data：Mashan_ROI
Output Raster：Extract_Mash_V3
点击 Run

Extract_Mash_V3

裁剪后的 Extract_Mash_V3 包含 6 个波段，对应关系如下：

ArcGIS 中的名称	Sentinel-2 波段	光的类型	主要用途
Band_1	B2	🔵 蓝光	BSI 计算
Band_2	B3	🟢 绿光	MNDWI 计算
Band_3	B4	🔴 红光	NDVI 计算
Band_4	B8	🟤 近红外 (NIR)	NDVI / NDBI
Band_5	B11	🟫 短波红外1 (SWIR1)	MNDWI / NDBI / BSI
Band_6	B12	🟫 短波红外2 (SWIR2)	土壤湿度（备用）

8.3 基于 V3 重新计算 NDVI

因为底图从 4 波段升级为 6 波段，我们需要基于新的 Extract_Mash_V3 重新计算 NDVI：

在左侧 Contents 面板中点击选中 Extract_Mash_V3（让它成为活跃图层）
顶部菜单 → Imagery → Indices → NDVI
设置：NIR Band = 4，Red Band = 3
点击确定
生成后立刻右键图层 → Data → Export Raster，导出为永久 TIFF 文件（避免临时图层崩溃）

配色步骤：

右键 NDVI_Extract_Mash_V3.tif → Symbology
在 Color scheme 搜索栏输入 Condition Number，选中从红到绿的色带
勾选 Invert（反转）：绿色 = 高植被，红色 = 建设用地
底部 Mask 标签页 → 勾选 Display background value，值设为 0 → No Color（去黑底）

NDVI 配色效果

8.4 计算 MNDWI（水体指数）

NDVI 只能告诉我们植被的情况，但论文还需要精准识别水体。MNDWI（Modified Normalized Difference Water Index，改进型归一化水体指数）利用水的物理特性来实现这一目标。

公式：MNDWI = (绿光 - 短波红外) / (绿光 + 短波红外)

:::{.fold} 为什么这个公式能识别水体？

水有一个独特的物理特性：

强烈反射绿光（所以水看起来发蓝绿色）→ Band_2 值高
几乎完全吸收短波红外（红外穿不透水）→ Band_5 值极低

所以对于水面：(高 - 低) / (高 + 低) → MNDWI 接近 +1（地图上显示为亮白色）

而对于建筑和植被：短波红外的反射比绿光更强 → MNDWI 接近 -1 或 0（暗黑色） :::

使用 Raster Calculator 计算（Analysis → Tools → 搜索 Raster Calculator）：

搜索 Raster Calculator

:::{.fold} 为什么会出现两个 Raster Calculator？

两个 Raster Calculator

这两个功能几乎完全一样，只是来自不同的扩展模块：

Raster Calculator (Image Analyst Tools) → 来自 Image Analyst 扩展
Raster Calculator (Spatial Analyst Tools) → 来自 Spatial Analyst 扩展

公式语法和计算结果完全相同，随便选一个即可。 :::

在 Raster Calculator 的表达式框中输入：

python

Float(Raster("Extract_Mash_V3/Band_2") - Raster("Extract_Mash_V3/Band_5")) / Float(Raster("Extract_Mash_V3/Band_2") + Raster("Extract_Mash_V3/Band_5"))

:::{.fold} 这个公式每一部分是什么意思？

语法部分	含义
`Raster("Extract_Mash_V3/Band_2")`	从裁剪后的 6 波段栅格中取出第 2 波段（🟢 绿光）
`Raster("Extract_Mash_V3/Band_5")`	取出第 5 波段（🟫 短波红外 SWIR1）
`Band_2 - Band_5`	绿光减去短波红外（差值）
`Band_2 + Band_5`	绿光加上短波红外（总和）
`差值 / 总和`	归一化，结果限定在 -1 到 1 之间
`Float(...)`	把整数结果转成小数（浮点型），保留精度
:::

Output raster：MNDWI_V3
点击 Run

MNDWI 计算结果

计算成功，MNDWI 值范围 -0.70 ~ 0.71。地图上左侧太湖水面呈亮白色（MNDWI 接近 1），半岛内部陆地呈暗黑色（MNDWI 接近 0 或负值），水体被完美识别。

8.5 计算 NDBI（建筑指数）

最后一个指数——NDBI（Normalized Difference Built-up Index，归一化建筑指数），用于识别建设用地（房屋、道路、硬化地面）。

公式：NDBI = (短波红外 - 近红外) / (短波红外 + 近红外)

原理：建筑材料（水泥、沥青）在短波红外波段的反射率高于近红外，而植被恰好相反。

在 Raster Calculator 中输入：

python

Float(Raster("Extract_Mash_V3/Band_5") - Raster("Extract_Mash_V3/Band_4")) / Float(Raster("Extract_Mash_V3/Band_5") + Raster("Extract_Mash_V3/Band_4"))

Output raster：NDBI_V3
点击 Run

NDBI 计算结果

计算成功，NDBI 值范围 -0.59 ~ 0.55。城区和建筑密集区域偏亮，森林和水面偏暗，建筑用地被正确识别。

8.6 全套光谱指数总结

至此，三大生态分析核心指数已全部计算完毕：

指数	公式	数值范围	识别目标	存储位置
NDVI	$\frac{NIR - Red}{NIR + Red}$	-0.33 ~ 0.89	🌿 植被健康度	`NDVI_Extract_Mash_V3.tif`
MNDWI	$\frac{Green - SWIR1}{Green + SWIR1}$	-0.70 ~ 0.71	💧 水体	`mashan.gdb/MNDWI_V3`
NDBI	$\frac{SWIR1 - NIR}{SWIR1 + NIR}$	-0.59 ~ 0.55	🏢 建设用地	`mashan.gdb/NDBI_V3`

💡 MNDWI 和 NDBI 通过 Raster Calculator 直接输出到了 mashan.gdb，本身就是永久文件，无需像 NDVI 那样额外执行 Export Raster。

:::{.fold} 各个卫星都有什么不同？有哪些主流的卫星可以获取数据？

好问题！在做生态分析之前，搞清楚"天上有哪些卫星可以用"是最基本的知识。以下是免费开放的主流对地观测卫星：

🛰️ Landsat 系列（美国 NASA/USGS）

Landsat 是人类历史上运行时间最长的对地观测计划，从 1972 年至今，积累了超过 50 年的地表影像档案。

卫星	运行时间	分辨率	重访周期	在 GEE 中的数据集 ID
Landsat 5 (TM)	1984–2012	30m	16 天	`LANDSAT/LT05/C02/T1_L2`
Landsat 7 (ETM+)	1999–至今	30m	16 天	`LANDSAT/LE07/C02/T1_L2`
Landsat 8 (OLI)	2013–至今	30m	16 天	`LANDSAT/LC08/C02/T1_L2`
Landsat 9 (OLI-2)	2021–至今	30m	16 天	`LANDSAT/LC09/C02/T1_L2`

⚠️ Landsat 7 的致命缺陷：2003 年 5 月传感器发生硬件故障（SLC-off），此后所有影像都带有条纹状黑色条带（约 22% 数据丢失）。虽然可以通过插值算法修补，但一般能避免就避免。

适用场景：需要追溯 2000 年以前的历史变化时，Landsat 是唯一的选择。

🛰️ Sentinel-2（欧洲 ESA）

欧洲空间局的"哨兵"计划，2015 年发射，是目前免费卫星里分辨率最高的。

卫星	运行时间	分辨率	重访周期	在 GEE 中的数据集 ID
Sentinel-2A	2015–至今	10m	5 天（双星）	`COPERNICUS/S2_SR_HARMONIZED`
Sentinel-2B	2017–至今	10m	5 天（双星）	同上

优势：分辨率是 Landsat 的 3 倍（10m vs 30m），有 13 个波段，重访周期更短（5 天 vs 16 天）。

局限：只有 2015 年以后的数据，无法追溯更早的历史。

🛰️ 两者的波段对比

这是你论文中最关键的对照表——同一个指数，在不同卫星上要用不同的波段编号：

光谱	Sentinel-2	Landsat 8/9	Landsat 5/7
🔵 蓝光	B2 (490nm)	B2 (482nm)	B1 (485nm)
🟢 绿光	B3 (560nm)	B3 (562nm)	B2 (560nm)
🔴 红光	B4 (665nm)	B4 (655nm)	B3 (660nm)
🟤 近红外 (NIR)	B8 (842nm)	B5 (865nm)	B4 (830nm)
🟫 短波红外 (SWIR1)	B11 (1610nm)	B6 (1609nm)	B5 (1650nm)
🟫 短波红外 (SWIR2)	B12 (2190nm)	B7 (2201nm)	B7 (2220nm)

可以看到，虽然波段编号不同，但它们捕获的光的物理波长几乎完全一致。所以公式逻辑不变，只是编号要对应替换。

🛰️ 其他你可能听说过的卫星

卫星	国家	分辨率	特点
MODIS	美国	250m~1km	超大范围、每天覆盖，适合全球/大洲尺度，但太粗糙不适合半岛级别分析
高分系列 (GF)	中国	1m~16m	分辨率极高，但数据获取渠道复杂，不在 GEE 上
Planet	美国（商业）	3m	每天覆盖全球，但要付费

🎯 对于你的论文，最优方案

经过进一步调研，为了保证长达 20 多年的时序趋势分析的稳定性，学术界的主流共识是：

统一使用 Landsat 30m 数据集，是目前做 2000-2020 年长时序 LUCC 分析最稳妥的方案。

如果前期用 Landsat (30m)、后期换 Sentinel-2 (10m)，分辨率的突然跳变可能会导致伪变化（例如把"分辨率提高后多识别出的小路"误判为"新增的建设用地"）。

因此，我们最终确定的全统一 Landsat 方案如下：

年份	卫星	分辨率	数据集
2000	Landsat 5	30m	`LANDSAT/LT05/C02/T1_L2`
2005	Landsat 5	30m	`LANDSAT/LT05/C02/T1_L2`
2010	Landsat 5	30m	`LANDSAT/LT05/C02/T1_L2`
2015	Landsat 8	30m	`LANDSAT/LC08/C02/T1_L2`
2020	Landsat 8	30m	`LANDSAT/LC08/C02/T1_L2`
2025	Landsat 9	30m	`LANDSAT/LC09/C02/T1_L2`

💡 之前我们用 Sentinel-2 (10m) 做的 NDVI / MNDWI / NDBI 可以作为2023 年现状的高精度补充分析，与 Landsat 的长时序分析并行使用。

:::

:::{.fold} 这些数据是什么处理级别？相当于 Sentinel-2 的 L2A 吗？

数据集 ID 中的 T1_L2 揭示了答案：

C02 = Collection 2（USGS 最新的数据处理版本）
T1 = Tier 1（最高质量等级，几何精度最好的影像）
L2 = Level-2（已完成大气校正的地表反射率 Surface Reflectance）

因此，全部 6 年的数据都是 Level-2 地表反射率产品，已经完成了以下所有预处理步骤：

✅ 辐射校正 — DN 值已转换为物理量
✅ 几何精校正 — Tier 1 是最高定位精度
✅ 大气校正 — 已从 TOA（大气层顶部反射率）校正为 BOA（地表反射率）

Landsat 与 Sentinel-2 的级别命名对照

两个卫星体系的处理级别含义相同，但命名方式不同：

处理内容	Sentinel-2 命名	Landsat 命名	我们的数据
大气层顶部反射率（TOA）	L1C	T1 (Level-1)	—
地表反射率（BOA，已大气校正）	L2A	T1_L2 (Level-2)	✅ 就是这个

一句话总结：我们的 Landsat T1_L2 数据等同于 Sentinel-2 的 L2A 级别，是可以直接用于光谱分析和分类的最高质量产品，无需再做任何预处理。

:::

:::{.fold} 各个卫星的分辨率 10 米和 30 米是什么意思？

你的理解完全正确！让我用一个直观的类比来解释：

分辨率 = 每个像素代表地面上多大一块地

想象你用马赛克来拼一幅画：

10 米分辨率（Sentinel-2）：每一块马赛克小方块代表地面上 10m × 10m = 100 m² 的区域。相当于用小颗粒拼图，细节丰富。
30 米分辨率（Landsat）：每一块马赛克代表 30m × 30m = 900 m² 的区域。相当于用大颗粒拼图，比较粗糙。

也就是说，同样覆盖 1 km² 的区域：

	Sentinel-2 (10m)	Landsat (30m)
每个像素覆盖面积	100 m²	900 m²
1 km² 需要多少像素	10,000 个	约 1,111 个
精细度	能分辨一栋楼	只能分辨一个街区

实际效果差异

举个例子，马山半岛上有一条 20 米宽的公路：

10 米分辨率：公路占了 2 个像素的宽度，卫星能"看见"这条路 ✅
30 米分辨率：公路不到 1 个像素宽，它会和旁边的农田、树木"混"在同一个像素里，被平均掉 ❌

对你论文的影响

好消息是：你的论文分析的是半岛尺度的大面积土地变化（森林 vs 建设用地），不是精确到单栋建筑。30 米分辨率对于这个研究尺度完全够用，这也是为什么绝大多数 LUCC 论文都用 Landsat 30m 数据。

:::

8.7 MNDWI 可视化验证与分类方法讨论

为了直观验证 MNDWI 的计算结果，我们可以在 ArcGIS Pro 中对其进行配色：

模式：Classify（分类模式），设置为 2 个类
断点：设为 0
配色：≤0（陆地）设为透明或浅灰，>0（水体）设为蓝色

MNDWI 二分法配色设置去除背景值

MNDWI 水体识别结果

通过验证，研究区内的池塘、太湖水面等被成功识别为蓝色，证明公式计算正确。

:::{.fold} 质疑：这种阈值法（>0 就是水）真的可靠吗？

确实不够可靠！ MNDWI > 0 只是一个粗略的分界线：

MNDWI 值	通常对应	但也可能是...
> 0.3	✅ 几乎确定是水	极少误判
0 ~ 0.3	⚠️ 可能是水	也可能是潮湿泥土、阴影、深色屋顶
< 0	通常不是水	植被、建筑、农田

如果只依靠阈值法做分类，会产生大量误判。这就是为什么我们在论文的正式分类环节必须使用随机森林（Random Forest）等监督分类算法。

在后续的随机森林模型中，NDVI、MNDWI、NDBI 都不需要提前分类，而是保持其连续的原始数值（-1 到 +1），和 6 个原始光谱波段一起喂给模型：

某像素的 9 个特征输入 → 随机森林模型自动学习权重规则 → 输出类别（如：建筑）

例如：一个像素的 MNDWI = 0.1（看似像水），但它的 NDBI = 0.4（典型建筑特征），随机森林就会聪明地将其判断为"深色建筑"而不是水体。 :::

第 9 步：多年份时序分析路线图

在完成了 2023 年高精度现状分析（技术打样）后，我们将正式进入论文核心——时序演变分析。

研究目的：通过生成 2000-2025 年间 6 期 LUCC 分类图，计算各类土地的面积增减和相互转化（转移矩阵），量化论证"旅游开发如何蚕食生态环境"。

整体技术路线：

GEE 批量导出：导出 2000-2025 年共 6 期的 Landsat 6 波段原始影像（统一 30m 分辨率）。
ArcGIS 预处理：导入数据并使用 Extract by Mask 裁剪出马山半岛范围。
特征工程：利用 Raster Calculator 为每一年分别计算出 NDVI、MNDWI、NDBI 三大生态指数层。
随机森林监督分类：
- 将原始波段 + 3 个指数合并为多波段特征层
- 手动选择并标注训练样本（水体、植被、建设用地、农田）
- 运行 Classification Wizard 进行分类
空间统计与制图：
- 计算各年份土地分类面积表
- 计算土地利用转移矩阵
- （可选）叠合 DEM 数据，分析建设用地的地形扩张规律

9.1 导出 ROI 矢量边界

为了在 GEE 中精确裁剪研究区，需要将 ArcGIS Pro 中的 Mashan_ROI 矢量边界导出为 Shapefile 格式，然后上传到 GEE Assets。

从 ArcGIS Pro 导出 Shapefile

右键 Mashan_ROI 图层 → Data → Export Features

Export Features

默认导出路径指向 .gdb 数据库，GEE 不支持这种格式：

默认导出到 .gdb 数据库里

点击右侧 📁 文件夹图标，将路径改为项目文件夹下的 .shp 文件：

Mashan_ROI.shp

导出后会生成一组文件：

导出的文件

:::{.fold} Shapefile 为什么有这么多文件？

Shapefile 是 1990 年代 ESRI 设计的格式，一个"文件"实际上是一组文件：

文件	必需？	作用
`.shp`	✅ 必需	几何数据（多边形的形状坐标）
`.dbf`	✅ 必需	属性表（字段数据）
`.shx`	✅ 必需	索引文件（加速读取）
`.prj`	✅ 必需	坐标参考系信息
`.cpg`	可选	字符编码声明
`.sbn` / `.sbx`	可选	空间索引（加速查询）
`.shp.xml`	可选	元数据（ArcGIS 自动生成的描述信息）
`.sr.lock`	❌ 忽略	ArcGIS 的临时锁文件，关闭项目后会消失

上传到 GEE Assets 时，只需选中 .shp、.dbf、.shx、.prj 这 4 个必需文件即可。 :::

上传到 GEE Assets

打开 GEE Code Editor
左侧面板 → Assets 标签 → 点 NEW → Shape files
选中 .shp、.dbf、.shx、.prj 四个文件上传
等待 Ingesting 完成，获得 Asset ID（如 users/你的用户名/Mashan_ROI）

Assets -> New -> Shape files Upload a new shapefile asset

选择 4 个文件上传

等待 Ingesting 完成 View asset 预览一下

GEE 脚本中使用 ROI

上传完成后，Asset ID 为：

projects/gen-lang-client-0147864203/assets/Mashan_ROI

9.2 GEE 批量导出 Landsat 多年份影像

将以下脚本粘贴到 GEE Code Editor 中运行，即可批量导出 6 个年份的 Landsat 地表反射率数据。

脚本核心逻辑：

数据源：根据年份自动选择 Landsat 5 / 8 / 9
时间窗口：每年 6-9 月（植被生长季），云覆盖 <30%
处理：云掩膜 → 中值合成 → 应用缩放因子（转为真实反射率 0~1）
波段统一：Landsat 5 的 B1-B5,B7 和 Landsat 8/9 的 B2-B7 统一重命名为 Blue/Green/Red/NIR/SWIR1/SWIR2
输出：6 个 GeoTIFF，30m 分辨率，UTM 51N 投影，存入 Google Drive Mashan_Landsat 文件夹

:::{.fold} 完整 GEE 导出脚本

javascript

// =============================================================
// 马山半岛 LUCC 分析 —— Landsat 多年份批量导出脚本
// 用途：导出 2000/2005/2010/2015/2020/2025 六期 Landsat 地表反射率数据
// 输出：6 波段 GeoTIFF（蓝/绿/红/NIR/SWIR1/SWIR2），统一 30m 分辨率
// =============================================================

// ---------------------- 1. 定义研究区 ----------------------
// 使用上传到 GEE Assets 的精确 ROI 矢量边界
var roi = ee.FeatureCollection('projects/gen-lang-client-0147864203/assets/Mashan_ROI');

// 在地图上显示研究区
Map.centerObject(roi, 13);
Map.addLayer(roi, {color: 'red'}, '研究区范围');

// ---------------------- 2. 云掩膜函数 ----------------------

// Landsat 5 云掩膜（使用 QA_PIXEL 波段）
function maskL5(image) {
  var qa = image.select('QA_PIXEL');
  var cloudShadowBitMask = 1 << 3;
  var cloudBitMask = 1 << 4;
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
    .and(qa.bitwiseAnd(cloudBitMask).eq(0));
  
  // Landsat 5 L2 缩放因子：scale = 0.0000275, offset = -0.2
  return image.updateMask(mask)
    .select(['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B7'])
    .multiply(0.0000275).add(-0.2)
    .rename(['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2'])
    .copyProperties(image, ['system:time_start']);
}

// Landsat 8/9 云掩膜
function maskL89(image) {
  var qa = image.select('QA_PIXEL');
  var cloudShadowBitMask = 1 << 3;
  var cloudBitMask = 1 << 4;
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
    .and(qa.bitwiseAnd(cloudBitMask).eq(0));
  
  return image.updateMask(mask)
    .select(['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7'])
    .multiply(0.0000275).add(-0.2)
    .rename(['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2'])
    .copyProperties(image, ['system:time_start']);
}

// ---------------------- 3. 年份配置 ----------------------
var yearConfigs = [
  { year: 2000, collection: 'LANDSAT/LT05/C02/T1_L2', maskFn: maskL5, satellite: 'Landsat5' },
  { year: 2005, collection: 'LANDSAT/LT05/C02/T1_L2', maskFn: maskL5, satellite: 'Landsat5' },
  { year: 2010, collection: 'LANDSAT/LT05/C02/T1_L2', maskFn: maskL5, satellite: 'Landsat5' },
  { year: 2015, collection: 'LANDSAT/LC08/C02/T1_L2', maskFn: maskL89, satellite: 'Landsat8' },
  { year: 2020, collection: 'LANDSAT/LC08/C02/T1_L2', maskFn: maskL89, satellite: 'Landsat8' },
  { year: 2025, collection: 'LANDSAT/LC09/C02/T1_L2', maskFn: maskL89, satellite: 'Landsat9' },
];

// ---------------------- 4. 批量处理与导出 ----------------------
yearConfigs.forEach(function(config) {
  var startDate = config.year + '-06-01';
  var endDate   = config.year + '-09-30';
  
  var composite = ee.ImageCollection(config.collection)
    .filterBounds(roi)
    .filterDate(startDate, endDate)
    .filter(ee.Filter.lt('CLOUD_COVER', 30))
    .map(config.maskFn)
    .median()
    .clip(roi);
  
  composite = composite.clamp(0, 1);
  
  var visParams = {bands: ['Red', 'Green', 'Blue'], min: 0, max: 0.3};
  Map.addLayer(composite, visParams, config.satellite + '_' + config.year);
  
  print(config.year + ' (' + config.satellite + ') composite:', composite);
  
  Export.image.toDrive({
    image: composite,
    description: 'Mashan_' + config.satellite + '_' + config.year,
    folder: 'Mashan_Landsat',
    fileNamePrefix: 'Mashan_' + config.year + '_SR',
    region: roi,
    scale: 30,
    crs: 'EPSG:32651',
    maxPixels: 1e9,
    fileFormat: 'GeoTIFF'
  });
});

:::

:::{.fold} Landsat 5 和 Landsat 8/9 的波段编号为什么不同？

这是因为 Landsat 8/9 增加了一个海岸气溶胶波段（Coastal Aerosol），插在了最前面，导致后续所有波段编号都往后移了一位：

光谱	Landsat 5 (TM)	Landsat 8/9 (OLI)
🔵 蓝光	SR_B1	SR_B2
🟢 绿光	SR_B2	SR_B3
🔴 红光	SR_B3	SR_B4
🟤 近红外	SR_B4	SR_B5
🟫 SWIR1	SR_B5	SR_B6
🟫 SWIR2	SR_B7	SR_B7

脚本中通过 .rename() 统一命名为 Blue/Green/Red/NIR/SWIR1/SWIR2，确保后续处理代码通用。 :::

运行脚本后，在右侧 Tasks 面板中逐个点击 Run 启动导出任务，导出完成后在 Google Drive 的 Mashan_Landsat 文件夹中下载 6 个 .tif 文件。

运行导出任务

将脚本粘贴到 GEE Code Editor 后点击 Run，右侧 Tasks 面板会出现 6 个待导出的任务。逐个点击 Run 启动导出：

Mashan_Landsat5_20xx run

2010 年导出失败与修复

2010 年的导出任务报错：Error: Can't get band number 0. Image has no bands.

Error: Can't get band number 0. Image has no bands.

原因是 Landsat 5 在 2010 年已接近退役（2011 年底故障），当年 6-9 月覆盖马山半岛且云量 <30% 的影像数量为 0。

修复方案：将时间窗口从 6-9 月扩大到 4-10 月，云量限制从 30% 放宽到 50%。修复后找到 3 张可用影像，成功完成导出：

可用数量为3 2010_fix

下载并归档

6 个年份全部导出完成后，在 Google Drive Mashan_Landsat 文件夹中下载 .tif 文件，并复制到 ArcGIS 项目目录：

bash

mkdir -p ~/Documents/ArcGIS/Projects/mashan/01_Raw_Data/Raster/Landsat

bash

cp ~/Library/CloudStorage/GoogleDrive-*/My\ Drive/Mashan_Landsat/Mashan_*.tif \
   ~/Documents/ArcGIS/Projects/mashan/01_Raw_Data/Raster/Landsat/

下载完成 6 个图像 cp 命令复制完成

文件名	年份	卫星	大小
`Mashan_2000_SR.tif`	2000	Landsat 5	1.6 MB
`Mashan_2005_SR.tif`	2005	Landsat 5	1.5 MB
`Mashan_2010_SR.tif`	2010	Landsat 5	1.8 MB
`Mashan_2015_SR.tif`	2015	Landsat 8	1.9 MB
`Mashan_2020_SR.tif`	2020	Landsat 8	1.9 MB
`Mashan_2025_SR.tif`	2025	Landsat 9	1.9 MB

💡 GEE 导出时已使用 .clip(roi) 裁剪到马山半岛边界，后续无需再做 Extract by Mask。

在 ArcGIS Pro 的 Catalog 面板中刷新后，可以看到这 6 个 Landsat 影像已经成功出现在项目中。

拖动或右键添加到当前地图 Landsat 影像已存入项目文件夹

TIP

加载到地图中 此时这些文件仅存在于硬盘文件夹里。要真正在地图上使用它们，请全选这 6 个 .tif 文件，右键选择 Add To Current Map（添加到当前地图），或者用鼠标直接拖拽到左侧的 Contents 面板中。拖入时如果提示 Build Pyramids and Calculate Statistics，直接点击 OK 即可。

:::{.fold} 图层太多，Contents 显示太杂乱，如何处理？

为了让图层列表保持整洁，我们可以将这 6 张原始地表反射率底图合并到一个图层组 (Group Layer) 中。

在左侧 Contents 面板中，按住 Shift 全选这 6 个新加载的图层。
右键选择 Group。
双击新建的 Group Layer，将其重命名为 Landsat_2000_2025，以区分后续即将生成的特征指数。

影像太多，创建一个 Group

注：你可以取消勾选该组来隐藏这些图层。这完全不影响后续 Python 脚本在后台读取并处理它们。 :::

:::{.fold} 导入的 6 个影像如何去掉黑色的区域？

在进行批量计算之前，你可能会注意到一个“碍眼”的细节：马山半岛以外的背景区域呈现出大片的黑色。这不仅影响视觉，更让人担心它是否会干扰后续的分析。

为了隐藏这片黑色，我们的第一反应通常是去 Symbology (符号系统) 面板的 Mask 选项卡里，勾选 Display background value 0,0,0 as Transparent。但奇怪的是，打勾之后竟然毫无变化，黑边依然存在！

这是一个极其关键的警告信号！🚨

如果勾选了 0 依然没有变透明，这说明这片黑色的像元值根本不是 0！ 如果在 ArcGIS 眼里它是一个具体的数字，那么待会儿跑 Python 代码算 NDVI 时，背景就会被强行算出 0.0 这种假数据，直接毁掉我们后面的随机森林分类。

为了查明真相，我们使用了 ArcGIS 顶部的 Explore (浏览) 工具，直接点击地图上的黑块。真相大白：

探索出黑块的真面目是 nan

弹窗显示各个波段的值全是 nan (Not a Number)。这是浮点型数据表示 NoData（无数据）的最高级标准。

这就完美解释了刚才发生的所有悬案：

为什么勾选背景 0 无效？ 因为那片黑色的真实身份是 nan，而不是 0。ArcGIS 的输入框只能匹配纯数字，匹配不到 nan，所以无法变透明。它只是一层渲染 nan 时的默认黑色皮肤罢了。
为什么完全不需要处理它？ 因为 GEE 导出时非常严谨地把非研究区标记为了 nan。后续 ArcPy 在进行像素加减乘除时，任何数字遭遇 nan 其结果依然是 nan。

结论：这 6 张影像是极高规格的科研级数据。千万不要试图把它强制设为 0（那代表真实的极深水体）。直接在左侧取消勾选隐藏这些原始图层，眼不见为净，然后放心地去跑下一步的 Python 代码吧！算出的新图层自动就是透明边界。 :::

:::{.fold} 为什么导入的 TIF 图像下方有三个颜色？RGB 是什么意思？

仔细观察 Contents 面板中刚刚加载的图层，你会发现每个 .tif 文件展开后，都带有红、绿、蓝（RGB）三个颜色的通道。

RGB

这其实是遥感软件在进行**“多波段色彩合成”（RGB Composite）**。

这不是一张普通的“照片”：我们从 GEE 导出的 .tif 文件并不是一张可以直接双击查看的 JPG 图片，而是一个包含 6 层数据的“数字三明治”（里面包含了 Blue, Green, Red, NIR, SWIR1, SWIR2 共 6 个波段的数值）。
屏幕显色的局限性：我们的电脑显示器只能通过红 (Red)、绿 (Green)、蓝 (Blue) 三原色来发光。因此，为了让你能“看见”这 6 层数据，ArcGIS 必须从中挑出 3 个波段，分别塞进屏幕的红、绿、蓝三个发光通道里。
为什么默认的颜色看起来很诡异？：默认情况下，ArcGIS 会按顺序把文件的第 1、2、3 个波段，强制分配给屏幕的 R、G、B 通道。但回忆一下我们的 GEE 代码，我们的波段排序是 [Blue, Green, Red...]。这就导致 ArcGIS 搞了一个大乌龙：
- 屏幕的 红光 (Red) = 卫星的 蓝光 (Blue) 波段
- 屏幕的 绿光 (Green) = 卫星的 绿光 (Green) 波段
- 屏幕的 蓝光 (Blue) = 卫星的 红光 (Red) 波段
这就叫作假彩色合成 (False Color Composite)。如果你想看到符合人类肉眼认知的“真彩色地图”，只需要在右侧的 Symbology (符号系统) 中，把 Red 通道设为 Red 波段，Blue 通道设为 Blue 波段（把被颠倒的红蓝通道换回来），地图的颜色瞬间就会变得极其自然啦！ :::

第 10 步：使用 ArcPy 批量计算特征指数

我们现在有 6 个年份的原始遥感影像。根据路线图，我们需要为每一年计算三个核心生态指数：

NDVI (植被，需要 Red 和 NIR 波段)
MNDWI (水体，需要 Green 和 SWIR1 波段)
NDBI (建设用地，需要 SWIR1 和 NIR 波段)

如果手动在 Raster Calculator 里重复点 18 次，既费时又容易选错波段。作为进阶用法，我们将使用 ArcGIS Pro 内置的 Python 窗口 一键批量完成计算。

10.1 打开 Python 窗口

在 ArcGIS Pro 顶部菜单栏点击 Analysis 选项卡
点击 Python 图标，打开 Python 交互窗口

10.2 运行批量计算脚本

Python Window 输入代码

将以下代码复制并粘贴到 Python 窗口中，然后按回车运行。这段代码会自动遍历地图中加载的 6 个年份的 Landsat 影像，提取对应的波段，利用公式计算出 18 个指数图层，并将结果自动保存到默认数据库（mashan.gdb）中。

python

import arcpy
from arcpy.sa import *
import os

# 签出 Spatial Analyst 扩展权限
arcpy.CheckOutExtension("Spatial")

# 获取当前项目、当前地图以及默认的 .gdb 数据库
aprx = arcpy.mp.ArcGISProject("CURRENT")
m = aprx.activeMap
out_gdb = aprx.defaultGeodatabase

print(f"输出数据库为: {out_gdb}")

# 遍历地图中所有以 "Mashan_" 开头、以 "_SR.tif" 结尾的图层
for lyr in m.listLayers("Mashan_*_SR.tif"):
    # 从图层名提取年份，例如 "Mashan_2000_SR.tif" -> "2000"
    year = lyr.name.split("_")[1]
    print(f"\n开始处理 {year} 年影像...")
    
    # 获取原始 .tif 文件的硬盘绝对路径
    tif_path = lyr.dataSource
    
    # 提取所需波段 
    # (我们在 GEE 中使用了 .rename 重新命名了波段)
    green = Raster(os.path.join(tif_path, "Green"))
    red   = Raster(os.path.join(tif_path, "Red"))
    nir   = Raster(os.path.join(tif_path, "NIR"))
    swir1 = Raster(os.path.join(tif_path, "SWIR1"))
    
    # 1. 计算 NDVI = (NIR - Red) / (NIR + Red)
    ndvi = Float(nir - red) / Float(nir + red)
    ndvi_name = f"NDVI_{year}"
    ndvi.save(os.path.join(out_gdb, ndvi_name))
    print(f"  ✅ {ndvi_name} 已保存")
    
    # 2. 计算 MNDWI = (Green - SWIR1) / (Green + SWIR1)
    mndwi = Float(green - swir1) / Float(green + swir1)
    mndwi_name = f"MNDWI_{year}"
    mndwi.save(os.path.join(out_gdb, mndwi_name))
    print(f"  ✅ {mndwi_name} 已保存")
    
    # 3. 计算 NDBI = (SWIR1 - NIR) / (SWIR1 + NIR)
    ndbi = Float(swir1 - nir) / Float(swir1 + nir)
    ndbi_name = f"NDBI_{year}"
    ndbi.save(os.path.join(out_gdb, ndbi_name))
    print(f"  ✅ {ndbi_name} 已保存")

print("\n🎉 全部年份的 NDVI、MNDWI、NDBI 计算完成！请在 Catalog 的 mashan.gdb 中查看结果。")

10.3 检查结果

代码运行大约需要几分钟。运行完成后：

打开右侧 Catalog 面板，展开 mashan.gdb 数据库
如果看不到新生成的文件，右键点击 mashan.gdb 选择 Refresh (刷新)
你会看到 18 个崭新的图层（如 NDVI_2000、MNDWI_2025 等）已经整整齐齐地保存在数据库中了！

Output

第 11 步：使用 ArcPy 批量合成多波段特征 (Composite Bands)

在 ArcGIS Pro 的图像分类向导（Image Classification Wizard）中，分类器每次只能输入一个多波段图层。而目前我们每个年份拥有 4 个分散的图层（1 个原始 6 波段图像 + 3 个特征指数图像）。

为了让随机森林能同时学习到光谱信息和生态指数特征，我们需要把这 4 个图层“像三明治一样”合并成一个完整的**【9波段超级特征图像 (Composite Bands)】**，然后把这个特征图像输入给机器学习分类器。

在 Python 窗口中，运行以下代码进行批量合成：

python

import arcpy
import os

aprx = arcpy.mp.ArcGISProject("CURRENT")
m = aprx.activeMap
out_gdb = aprx.defaultGeodatabase

arcpy.env.overwriteOutput = True
print(f"Output GDB: {out_gdb}")

years = ["2000", "2005", "2010", "2015", "2020", "2025"]

for year in years:
    print(f"\nCompositing {year} features...")
    
    # 1. 找到该年份的原始多波段影像
    raw_lyr = m.listLayers(f"Mashan_{year}_SR.tif")[0]
    raw_path = raw_lyr.dataSource
    
    # 2. 找到对应的三个指数图层路径 (在 gdb 中)
    ndvi_path = os.path.join(out_gdb, f"NDVI_{year}")
    mndwi_path = os.path.join(out_gdb, f"MNDWI_{year}")
    ndbi_path = os.path.join(out_gdb, f"NDBI_{year}")
    
    # 3. 组合输入列表 (6 个原始波段 + 3 个指数波段 = 9 波段)
    in_rasters = [raw_path, ndvi_path, mndwi_path, ndbi_path]
    
    # 4. 输出路径
    out_name = f"RF_Features_{year}"
    out_path = os.path.join(out_gdb, out_name)
    
    # 5. 执行 Composite Bands
    arcpy.management.CompositeBands(in_rasters, out_path)
    
    print(f"  [OK] {out_name} generated with 9 bands")

print("\nAll composite bands ready for Random Forest!")

Output

执行完毕后，mashan.gdb 中将生成 6 个名为 RF_Features_20xx 的 9 波段栅格图层。这正是机器学习分类器的完美输入源！

共 9 个图像分为两个 Group

第 12 步：阈值自动标样本 + 随机森林分类

本步骤是整个 LUCC 分析的核心环节。我们采用"混合法"：先用严格的指数阈值自动生成高纯度训练样本，再将这些样本输入随机森林分类器，让机器自动学习并分类整张影像。

:::{.fold} 为什么不直接用阈值法分类，而要多此一举训练随机森林？

阈值法的致命缺陷：处理不了"模糊地带"。

比如你设 MNDWI > 0.3 = 水体。那 MNDWI = 0.28 的像素呢？它可能是浅水、湿地、或者刚下完雨的泥地。阈值法只能"一刀切"，要么全算水体，要么全不算。

而随机森林的做法是：你只需要给它那些 MNDWI > 0.4 的"绝对是水"的像素作为样本，它会自己去学习"水体在 9 个波段上的综合特征模式"。然后它回过头来处理 MNDWI = 0.28 的像素时，不只看 MNDWI 这一个值，而是同时看这个像素在 Blue、Green、Red、NIR、SWIR1、SWIR2、NDVI、MNDWI、NDBI 这 9 个维度上的综合表现，做出更准确的判断。

阈值法是"用一把尺子量"，随机森林是"用九把尺子同时量"。

:::

12.1 定义分类体系

我们将马山半岛的土地覆盖分为 4 个一级类别：

类别 ID	名称	含义	阈值规则（用于自动提取纯净样本）
1	🌊 水体	海域、湖泊、河流、水库	`MNDWI > 0.3`
2	🌿 植被	林地、草地、农田	`NDVI > 0.5`
3	🏘️ 建设用地	城镇、道路、工业区	`NDBI > 0.05` 且 `NDVI < 0.2`
4	🟫 裸地	未利用地、荒地、裸露岩石	`NDVI < 0.1` 且 `MNDWI < 0` 且 `NDBI < 0`

:::{.fold} 这些阈值是怎么确定的？

这些阈值的设定原则是**"宁缺毋滥"——我们不追求覆盖所有该类地物，而是只提取那些百分之百确定**的纯净像素。

水体 MNDWI > 0.3：MNDWI 大于 0.3 的像素，在物理上几乎只可能是开阔水面。浅水、湿地的 MNDWI 通常在 0.1-0.3 之间，我们故意不要它们。
植被 NDVI > 0.5：NDVI 大于 0.5 意味着非常茂密的绿色植被（通常是林冠中心）。稀疏草地的 NDVI 一般在 0.2-0.4，我们也不要。
建设用地 NDBI > 0.05 且 NDVI < 0.2：NDBI 为正说明地表有较多的建筑材料（水泥、沥青），同时 NDVI 很低排除了城市绿化带。
裸地：三个指数都很低，说明既不是水、也不是植被、也不是建筑。

这些阈值在实际运行后，如果某个类别的样本点太少（不足 100 个），可以适当放宽阈值；如果太多（超过 5000 个），可以收紧。

:::

12.2 自动生成训练样本

将以下脚本粘贴到 ArcGIS Pro 的 Python 窗口中运行。脚本会自动：

对每个年份，用阈值从 NDVI/MNDWI/NDBI 中提取"纯净区域"
将纯净区域转为多边形
在多边形内随机撒点
为每个点标注类别标签

python

import arcpy
import os
from arcpy.sa import *

arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = True

# ---- 配置 ----
out_gdb = r"C:\Mac\Home\Documents\ArcGIS\Projects\mashan\mashan.gdb"
arcpy.env.workspace = out_gdb

years = [2000, 2005, 2010, 2015, 2020, 2025]

# 阈值配置（可根据实际效果微调）
thresholds = {
    "Water":    {"rule": "MNDWI > 0.3", "class_id": 1},
    "Veg":      {"rule": "NDVI > 0.5",  "class_id": 2},
    "Built":    {"rule": "NDBI > 0.05 AND NDVI < 0.2", "class_id": 3},
    "Bare":     {"rule": "NDVI < 0.1 AND MNDWI < 0 AND NDBI < 0", "class_id": 4},
}

samples_per_class = 500  # 每个类别每年生成的随机点数

print(f"Output GDB: {out_gdb}")

for year in years:
    print(f"\n{'='*50}")
    print(f"Processing {year}...")

    ndvi  = Raster(f"NDVI_{year}")
    mndwi = Raster(f"MNDWI_{year}")
    ndbi  = Raster(f"NDBI_{year}")

    # ---- 创建各类别的纯净区域掩膜 ----
    masks = {}
    masks["Water"] = Con(mndwi > 0.3, 1)
    masks["Veg"]   = Con(ndvi > 0.5, 2)
    masks["Built"] = Con((ndbi > 0.05) & (ndvi < 0.2), 3)
    masks["Bare"]  = Con((ndvi < 0.1) & (mndwi < 0) & (ndbi < 0), 4)

    # ---- 合并所有训练点 ----
    all_points = None

    for class_name, mask_raster in masks.items():
        class_id = thresholds[class_name]["class_id"]
        temp_mask = f"temp_mask_{class_name}_{year}"
        temp_poly = f"temp_poly_{class_name}_{year}"
        temp_pts  = f"temp_pts_{class_name}_{year}"

        try:
            # 保存掩膜栅格
            mask_raster.save(os.path.join(out_gdb, temp_mask))

            # 栅格转多边形
            arcpy.conversion.RasterToPolygon(
                in_raster=temp_mask,
                out_polygon_features=temp_poly,
                simplify="NO_SIMPLIFY"
            )

            # 在多边形内随机撒点
            arcpy.management.CreateRandomPoints(
                out_path=out_gdb,
                out_name=temp_pts,
                constraining_feature_class=temp_poly,
                number_of_points_or_field=samples_per_class
            )

            # 添加类别字段并赋值
            arcpy.management.AddField(temp_pts, "ClassID", "SHORT")
            arcpy.management.AddField(temp_pts, "ClassName", "TEXT", field_length=20)
            with arcpy.da.UpdateCursor(temp_pts, ["ClassID", "ClassName"]) as cursor:
                for row in cursor:
                    row[0] = class_id
                    row[1] = class_name
                    cursor.updateRow(row)

            # 合并到总点集
            if all_points is None:
                arcpy.management.CopyFeatures(temp_pts, f"TrainingSamples_{year}")
                all_points = f"TrainingSamples_{year}"
            else:
                arcpy.management.Append(temp_pts, all_points, "NO_TEST")

            count = int(arcpy.management.GetCount(temp_pts)[0])
            print(f"  [OK] {class_name}: {count} points")

        except Exception as e:
            print(f"  [WARN] {class_name}: {e}")

        finally:
            # 清理临时数据
            for temp in [temp_mask, temp_poly, temp_pts]:
                try:
                    arcpy.management.Delete(temp)
                except:
                    pass

    total = int(arcpy.management.GetCount(f"TrainingSamples_{year}")[0])
    print(f"  >> Total training samples for {year}: {total}")

print("\n" + "="*50)
print("All training samples generated!")

运行完毕后，mashan.gdb 中将生成 6 个训练样本点图层：TrainingSamples_2000 ~ TrainingSamples_2025，每个包含约 2000 个点（4 类 × 500 点/类）。

:::{.fold} 可以验证一下生成的样本点吗？

在 ArcGIS Pro 中，右键点击 TrainingSamples_2020 → Open Attribute Table，你应该会看到类似这样的结构：

OID	Shape	ClassID	ClassName
1	Point	1	Water
2	Point	1	Water
...	...	...	...
501	Point	2	Veg
...	...	...	...

你也可以将这些点叠加到原始影像上，目视检查它们是否确实落在了正确的地物上。如果发现某些点明显标错（比如"水体"点落在了陆地上），说明阈值需要收紧。

:::

12.3 训练随机森林分类器

用生成的样本点训练随机森林模型：

python

import arcpy
from arcpy.sa import *

arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = True

out_gdb = r"C:\Mac\Home\Documents\ArcGIS\Projects\mashan\mashan.gdb"
ecd_dir = r"C:\Mac\Home\Documents\ArcGIS\Projects\mashan"
years = [2000, 2005, 2010, 2015, 2020, 2025]

for year in years:
    print(f"\nTraining RF classifier for {year}...")

    in_raster  = f"{out_gdb}\\RF_Features_{year}"
    in_samples = f"{out_gdb}\\TrainingSamples_{year}"
    out_ecd    = f"{ecd_dir}\\RF_Model_{year}.ecd"

    arcpy.sa.TrainRandomTreesClassifier(
        in_raster=in_raster,
        in_training_features=in_samples,
        out_classifier_definition=out_ecd,
        in_additional_raster=None,
        max_num_trees=100,
        max_tree_depth=30,
        max_samples_per_class=1000
    )
    print(f"  [OK] Model saved: RF_Model_{year}.ecd")

print("\nAll models trained!")

:::{.fold} 这些参数是什么意思？

参数	值	含义
`max_num_trees`	100	森林里种 100 棵决策树。越多越稳定，但越慢。100 是学术界最常用的默认值
`max_tree_depth`	30	每棵树最多分裂 30 层。越深 → 模型越复杂 → 越容易过拟合。30 对 4 个类别绰绰有余
`max_samples_per_class`	1000	每个类别最多取 1000 个样本参与训练。防止某个类别样本过多导致偏向
输出 `.ecd` 文件	—	Esri Classifier Definition 文件，保存了训练好的模型。可以复用到其他影像上

:::

12.4 执行分类

用训练好的模型对每年的 9 波段影像进行逐像素分类：

python

import arcpy
from arcpy.sa import *

arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = True

out_gdb = r"C:\Mac\Home\Documents\ArcGIS\Projects\mashan\mashan.gdb"
ecd_dir = r"C:\Mac\Home\Documents\ArcGIS\Projects\mashan"
years = [2000, 2005, 2010, 2015, 2020, 2025]

for year in years:
    print(f"\nClassifying {year}...")

    in_raster = f"{out_gdb}\\RF_Features_{year}"
    in_ecd    = f"{ecd_dir}\\RF_Model_{year}.ecd"
    out_class = f"{out_gdb}\\LULC_{year}"

    classified = arcpy.sa.ClassifyRaster(
        in_raster=in_raster,
        in_classifier_definition=in_ecd
    )
    classified.save(out_class)
    print(f"  [OK] Classification saved: LULC_{year}")

print("\nAll years classified!")
print("Result layers: LULC_2000, LULC_2005, ..., LULC_2025")

运行完毕后，mashan.gdb 中将生成 6 张分类结果图：LULC_2000 ~ LULC_2025。每个像素的值为 1-4，分别对应水体、植被、建设用地、裸地。

12.5 分类后处理

原始分类结果通常包含零散的错分像素（"椒盐噪声"），需要平滑处理：

python

from arcpy.sa import *

arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = True

out_gdb = r"C:\Mac\Home\Documents\ArcGIS\Projects\mashan\mashan.gdb"
years = [2000, 2005, 2010, 2015, 2020, 2025]

for year in years:
    print(f"\nPost-processing {year}...")

    raw = Raster(f"{out_gdb}\\LULC_{year}")

    # 众数滤波：将孤立像素替换为周围 8 邻域的众数
    filtered = MajorityFilter(raw, "EIGHT", "MAJORITY")

    # 边界清理：平滑类别边界
    cleaned = BoundaryClean(filtered, "NO_SORT", "TWO_WAY")

    cleaned.save(f"{out_gdb}\\LULC_{year}_clean")
    print(f"  [OK] LULC_{year}_clean saved")

print("\nPost-processing complete!")

:::{.fold} 众数滤波和边界清理在做什么？

众数滤波（Majority Filter）

想象分类结果是一幅马赛克图。如果一个"植被"像素被 8 个"建设用地"像素包围，它很可能是分错了。众数滤波会检查每个像素的 8 个邻居，如果多数邻居都是另一个类别，就把这个"孤独的像素"替换掉。

边界清理（Boundary Clean）

两个类别之间的边界通常是锯齿状的（因为像素是方块）。边界清理会让较大的地物区域"吞噬"边界处的零星像素，使得边界更加平滑自然。

这两步不会改变大面积区域的分类结果，只会清理掉边缘的零散噪点。

:::

12.6 设置符号化（配色）

为分类结果设置直观的颜色方案：

python

import arcpy

aprx = arcpy.mp.ArcGISProject("CURRENT")
m = aprx.activeMap

# 颜色方案
colors = {
    1: {"name": "Water",   "rgb": [0, 112, 255]},     # 蓝色
    2: {"name": "Veg",     "rgb": [56, 168, 0]},       # 绿色
    3: {"name": "Built",   "rgb": [255, 0, 0]},        # 红色
    4: {"name": "Bare",    "rgb": [255, 235, 175]},     # 浅黄色
}

for year in [2000, 2005, 2010, 2015, 2020, 2025]:
    lyr_name = f"LULC_{year}_clean"
    layers = m.listLayers(lyr_name)
    if layers:
        lyr = layers[0]
        sym = lyr.symbology
        if hasattr(sym, 'updateColorizer'):
            sym.updateColorizer('RasterUniqueValueColorizer')
            lyr.symbology = sym
            print(f"  [OK] Symbology updated for {lyr_name}")

print("Done! Manually adjust colors if needed via Layer Properties > Symbology.")

如果 ArcPy 自动配色效果不理想，可以手动操作：右键图层 → Symbology → Unique Values → 为 1/2/3/4 分别设置蓝/绿/红/浅黄色。

第 1 步：获取“矢量边界”数据 ​

从 Overpass Turbo 导出为标准格式： ​

在 ArcGIS 中删除项目的步骤： ​

QGIS：散养模式（指针链接） ​

ArcGIS Pro：全家桶模式（高内聚） ​

第 2 步：将 GeoJSON 导入 ArcGIS 并生成多边形 (ROI) ​

2.1 JSON To Features (导入线条) ​

2.2 Feature To Polygon (线转面，清理悬挂线段) ​

第 3 步：使用 Google Earth Engine (GEE) 获取无云卫星影像 ​

3.1 在 GEE 中框选研究区 ​

3.2 运行 GEE 代码截取 Sentinel-2 影像 ​

顺带提一个有趣的细节 ​

3.3 将影像移入本地 ArcGIS 目录 ​

第 4 步：在 ArcGIS 中加载影像并按掩膜提取 ​

4.1 避坑：跨系统共享路径 Bug 与火星坐标系偏移 ​

4.2 加载 V2 影像与生成附属文件 (.aux) ​

4.3 执行 Extract by Mask（按掩膜提取） ​

第 5 步：提取大自然的“绿色指纹” (计算 NDVI) ​

光的本质是电磁波 ​

为什么植物对这两种光的反应如此不同？ ​

5.1 认清波段序号并一键计算 ​

5.2 避坑：临时图层断连与“转正大法” ​

5.3 生态学配色与掩膜去黑底 ​

分辨率 = 每个像素代表地面上多大一块地 ​

实际效果差异 ​

对论文的影响 ​

5.1 MNDWI / NDBI 可视化验证 ​

5.2 指数的定位：中间产品，而非最终结论 ​

第 6 步：多年份数据方案与技术路线 ​

6.1 为什么需要多年份数据？ ​

6.2 分辨率统一方案 ​

6.3 分类方法：随机森林监督分类 ​

6.4 完整技术路线 ​

6.5 论文中需要的分析与结论 ​

📊 面积统计（基础数据） ​

🔄 土地转移矩阵（核心证据） ​

🏔️ DEM 地形叠加分析 ​

5.5 按掩膜提取 DEM 与视觉误区避坑 ​

第 7 步：土地利用分类 (LUCC) ​

第 8 步：补全波段，计算全套光谱指数 ​

8.1 在 GEE 中重新导出 6 波段影像（V3） ​

8.2 导入 ArcGIS 并按掩膜裁剪 ​

8.3 基于 V3 重新计算 NDVI ​

8.4 计算 MNDWI（水体指数） ​

8.5 计算 NDBI（建筑指数） ​

8.6 全套光谱指数总结 ​

🛰️ Landsat 系列（美国 NASA/USGS） ​

🛰️ Sentinel-2（欧洲 ESA） ​

🛰️ 两者的波段对比 ​

🛰️ 其他你可能听说过的卫星 ​

🎯 对于你的论文，最优方案 ​

Landsat 与 Sentinel-2 的级别命名对照 ​

分辨率 = 每个像素代表地面上多大一块地 ​

实际效果差异 ​

对你论文的影响 ​

8.7 MNDWI 可视化验证与分类方法讨论 ​

第 9 步：多年份时序分析路线图 ​

9.1 导出 ROI 矢量边界 ​

从 ArcGIS Pro 导出 Shapefile ​

上传到 GEE Assets ​

GEE 脚本中使用 ROI ​

9.2 GEE 批量导出 Landsat 多年份影像 ​

运行导出任务 ​

2010 年导出失败与修复 ​

下载并归档 ​

第 10 步：使用 ArcPy 批量计算特征指数 ​

10.1 打开 Python 窗口 ​

10.2 运行批量计算脚本 ​

10.3 检查结果 ​

第 11 步：使用 ArcPy 批量合成多波段特征 (Composite Bands) ​

第 12 步：阈值自动标样本 + 随机森林分类 ​

12.1 定义分类体系 ​

12.2 自动生成训练样本 ​

12.3 训练随机森林分类器 ​

12.4 执行分类 ​

12.5 分类后处理 ​

12.6 设置符号化（配色） ​

第 1 步：获取“矢量边界”数据

从 Overpass Turbo 导出为标准格式：

在 ArcGIS 中删除项目的步骤：

QGIS：散养模式（指针链接）

ArcGIS Pro：全家桶模式（高内聚）

第 2 步：将 GeoJSON 导入 ArcGIS 并生成多边形 (ROI)

2.1 JSON To Features (导入线条)

2.2 Feature To Polygon (线转面，清理悬挂线段)

第 3 步：使用 Google Earth Engine (GEE) 获取无云卫星影像

3.1 在 GEE 中框选研究区

3.2 运行 GEE 代码截取 Sentinel-2 影像

顺带提一个有趣的细节

3.3 将影像移入本地 ArcGIS 目录

第 4 步：在 ArcGIS 中加载影像并按掩膜提取

4.1 避坑：跨系统共享路径 Bug 与火星坐标系偏移

4.2 加载 V2 影像与生成附属文件 (.aux)

4.3 执行 Extract by Mask（按掩膜提取）

第 5 步：提取大自然的“绿色指纹” (计算 NDVI)

光的本质是电磁波

为什么植物对这两种光的反应如此不同？

5.1 认清波段序号并一键计算

5.2 避坑：临时图层断连与“转正大法”

5.3 生态学配色与掩膜去黑底

分辨率 = 每个像素代表地面上多大一块地

实际效果差异

对论文的影响

5.1 MNDWI / NDBI 可视化验证

5.2 指数的定位：中间产品，而非最终结论

第 6 步：多年份数据方案与技术路线

6.1 为什么需要多年份数据？

6.2 分辨率统一方案

6.3 分类方法：随机森林监督分类

6.4 完整技术路线

6.5 论文中需要的分析与结论

📊 面积统计（基础数据）

🔄 土地转移矩阵（核心证据）

🏔️ DEM 地形叠加分析

5.5 按掩膜提取 DEM 与视觉误区避坑

第 7 步：土地利用分类 (LUCC)

第 8 步：补全波段，计算全套光谱指数

8.1 在 GEE 中重新导出 6 波段影像（V3）

8.2 导入 ArcGIS 并按掩膜裁剪

8.3 基于 V3 重新计算 NDVI

8.4 计算 MNDWI（水体指数）

8.5 计算 NDBI（建筑指数）

8.6 全套光谱指数总结

🛰️ Landsat 系列（美国 NASA/USGS）

🛰️ Sentinel-2（欧洲 ESA）

🛰️ 两者的波段对比

🛰️ 其他你可能听说过的卫星

🎯 对于你的论文，最优方案

Landsat 与 Sentinel-2 的级别命名对照

分辨率 = 每个像素代表地面上多大一块地

实际效果差异

对你论文的影响

8.7 MNDWI 可视化验证与分类方法讨论

第 9 步：多年份时序分析路线图

9.1 导出 ROI 矢量边界

从 ArcGIS Pro 导出 Shapefile

上传到 GEE Assets

GEE 脚本中使用 ROI

9.2 GEE 批量导出 Landsat 多年份影像

运行导出任务

2010 年导出失败与修复

下载并归档

第 10 步：使用 ArcPy 批量计算特征指数

10.1 打开 Python 窗口

10.2 运行批量计算脚本

10.3 检查结果

第 11 步：使用 ArcPy 批量合成多波段特征 (Composite Bands)

第 12 步：阈值自动标样本 + 随机森林分类

12.1 定义分类体系

12.2 自动生成训练样本

12.3 训练随机森林分类器

12.4 执行分类

12.5 分类后处理

12.6 设置符号化（配色）