一、写在前面

最近在做人工林遥感语义分割任务时,我构建了一个面向高分辨率人工林分割的数据集。论文中的数据集覆盖亚洲、欧洲、北美和非洲四个区域,其中北美子集使用 NAIP 航空影像,并包含原始影像、NDVI、EVI 和像素级标注;其他区域主要来自高分辨率卫星影像和航空底图场景。大场景影像统一裁剪为 500×500 像素的小样本,最终形成四个区域共 4728 组 image-mask 样本。

数据集地址:

https://huggingface.co/datasets/jkhsdfja/ArtificialForest-HRSegDataset

这篇文章主要记录一个可复现的数据构建流程:如何使用 Google Earth Engine(GEE)寻找遥感图源、下载影像、计算植被指数,并进一步裁剪成语义分割数据集。

二、整体流程

整个流程可以分成 6 步:

  1.  在 GEE 中确定研究区域 ROI
  2. 选择合适的遥感影像源
  3. 进行云掩膜、影像合成和裁剪
  4. 计算 NDVI、EVI 等植被指数
  5. 导出 GeoTIFF 影像
  6. 本地裁剪为 500×500 图像块,并制作 mask 标签

需要注意的是,GEE 可以直接导出其数据目录中的公开影像,例如 Sentinel-2、Landsat、NAIP 等。
如果使用 Google Earth 网页端底图作为人工目视参考,一般不能直接从 GEE 导出 Google Earth 底图本身,更规范的做法是使用 GEE 中公开可下载的数据源,并在 QGIS、ArcGIS 或 LabelMe 中完成标注。

三、在 GEE 中获取图源


1. 北美区域:使用 NAIP 高分辨率影像

NAIP 是美国农业部提供的航空影像,分辨率较高,通常包含 RGB 和近红外波段,适合做植被分割、林地提取和人工林边界分析。

在 GEE 中可以这样加载 NAIP:

// ===============================
// 1. 定义研究区域
// ===============================
// 可以在 GEE 左侧 Geometry Imports 中手动绘制 ROI
var roi = geometry;

// ===============================
// 2. 加载 NAIP 影像
// ===============================
var naip = ee.ImageCollection("USDA/NAIP/DOQQ")
  .filterBounds(roi)
  .filterDate("2020-01-01", "2023-12-31");

// 多景影像镶嵌
var image = naip.mosaic().clip(roi);

// NAIP 常见波段:R, G, B, N
var rgb = image.select(["R", "G", "B"]);
var nir = image.select("N");

// ===============================
// 3. 计算 NDVI
// NDVI = (NIR - Red) / (NIR + Red)
// ===============================
var ndvi = image.normalizedDifference(["N", "R"]).rename("NDVI");

// ===============================
// 4. 计算 EVI
// EVI = 2.5 * (NIR - Red) / (NIR + 6*Red - 7.5*Blue + 1)
// ===============================
var evi = image.expression(
  "2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))",
  {
    "NIR": image.select("N"),
    "RED": image.select("R"),
    "BLUE": image.select("B")
  }
).rename("EVI");

// ===============================
// 5. 合并 RGB + NDVI + EVI
// ===============================
var outImage = rgb.addBands(ndvi).addBands(evi);

// 可视化
Map.centerObject(roi, 14);
Map.addLayer(rgb, {min: 0, max: 255}, "NAIP RGB");
Map.addLayer(ndvi, {min: -1, max: 1, palette: ["brown", "white", "green"]}, "NDVI");
Map.addLayer(evi, {min: -1, max: 1, palette: ["brown", "white", "blue"]}, "EVI");

// ===============================
// 6. 导出到 Google Drive
// ===============================
Export.image.toDrive({
  image: outImage,
  description: "NAIP_RGB_NDVI_EVI",
  folder: "GEE_Plantation_Dataset",
  fileNamePrefix: "naip_rgb_ndvi_evi",
  region: roi,
  scale: 1,
  maxPixels: 1e13,
  fileFormat: "GeoTIFF"
});

如果研究区在美国,优先推荐 NAIP,因为它包含近红外波段,可以直接计算 NDVI 和 EVI。

四、非北美区域:使用 Sentinel-2 获取可复现图源

如果研究区不在美国,NAIP 无法覆盖,可以使用 Sentinel-2。
Sentinel-2 分辨率虽然没有 NAIP 那么高,但它公开、稳定、可复现,并且具有红光、绿光、蓝光、近红外等波段,可以计算 NDVI 和 EVI。

// ===============================
// 1. 定义研究区域
// ===============================
var roi = geometry;

// ===============================
// 2. Sentinel-2 云掩膜函数
// ===============================
function maskS2Clouds(image) {
  var qa = image.select("QA60");

  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;

  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
    .and(qa.bitwiseAnd(cirrusBitMask).eq(0));

  return image.updateMask(mask).divide(10000);
}

// ===============================
// 3. 加载 Sentinel-2 SR 数据
// ===============================
var s2 = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
  .filterBounds(roi)
  .filterDate("2023-01-01", "2023-12-31")
  .filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 20))
  .map(maskS2Clouds);

// 中值合成,降低云和噪声影响
var image = s2.median().clip(roi);

// Sentinel-2 波段
// B2: Blue, B3: Green, B4: Red, B8: NIR
var rgb = image.select(["B4", "B3", "B2"]).rename(["R", "G", "B"]);

// ===============================
// 4. 计算 NDVI
// ===============================
var ndvi = image.normalizedDifference(["B8", "B4"]).rename("NDVI");

// ===============================
// 5. 计算 EVI
// ===============================
var evi = image.expression(
  "2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))",
  {
    "NIR": image.select("B8"),
    "RED": image.select("B4"),
    "BLUE": image.select("B2")
  }
).rename("EVI");

// ===============================
// 6. 合并导出
// ===============================
var outImage = rgb.addBands(ndvi).addBands(evi);

Map.centerObject(roi, 12);
Map.addLayer(rgb, {min: 0, max: 0.3}, "Sentinel-2 RGB");
Map.addLayer(ndvi, {min: -1, max: 1, palette: ["brown", "white", "green"]}, "NDVI");
Map.addLayer(evi, {min: -1, max: 1, palette: ["brown", "white", "blue"]}, "EVI");

Export.image.toDrive({
  image: outImage,
  description: "S2_RGB_NDVI_EVI",
  folder: "GEE_Plantation_Dataset",
  fileNamePrefix: "s2_rgb_ndvi_evi",
  region: roi,
  scale: 10,
  maxPixels: 1e13,
  fileFormat: "GeoTIFF"
});

这里导出的 GeoTIFF 包含 5 个通道:

R
G
B
NDVI
EVI

后续可以根据任务需求只使用 RGB,也可以将 NDVI、EVI 作为额外模态输入模型。

五、为什么要计算 NDVI 和 EVI?

在人工林分割任务中,仅使用 RGB 影像有时会遇到几个问题:

1. 人工林和天然林纹理相似
2. 草地、灌木和人工林光谱相近
3. 阴影、裸土和建筑物容易造成误分割
4. 不同地区、不同卫星影像之间存在亮度差异

NDVI 和 EVI 可以提供更直接的植被信息:

NDVI:反映植被绿度和长势
EVI:对高覆盖植被、土壤背景和大气影响更稳定

在我的论文中,模型采用了光学图像与 NDVI、EVI 对齐融合的方式,用于增强人工林区域的表征能力。

六、导出大图后,本地裁剪为 500×500 样本

论文中的大场景影像约为 5000×5000 像素,并统一裁剪成 500×500 像素的小图块,用于语义分割网络训练。

下面给出一个本地裁剪 GeoTIFF 的 Python 脚本。1. 安装依赖
pip install rasterio numpy tqdm
2. 裁剪影像为 500×500

import os
import rasterio
from rasterio.windows import Window
from tqdm import tqdm
import numpy as np

def crop_geotiff_to_tiles(
    image_path,
    out_dir,
    tile_size=500,
    stride=500,
    min_valid_ratio=0.6
):
    os.makedirs(out_dir, exist_ok=True)

    with rasterio.open(image_path) as src:
        width = src.width
        height = src.height
        profile = src.profile.copy()

        tile_id = 0

        for y in tqdm(range(0, height - tile_size + 1, stride)):
            for x in range(0, width - tile_size + 1, stride):
                window = Window(x, y, tile_size, tile_size)
                tile = src.read(window=window)

                # 判断有效像素比例,避免导出大量空白图块
                valid_mask = np.any(tile != 0, axis=0)
                valid_ratio = valid_mask.sum() / (tile_size * tile_size)

                if valid_ratio < min_valid_ratio:
                    continue

                tile_profile = profile.copy()
                tile_profile.update({
                    "height": tile_size,
                    "width": tile_size,
                    "transform": src.window_transform(window)
                })

                out_path = os.path.join(out_dir, f"image_{tile_id:05d}.tif")

                with rasterio.open(out_path, "w", **tile_profile) as dst:
                    dst.write(tile)

                tile_id += 1

    print(f"裁剪完成,共生成 {tile_id} 个图块。")


if __name__ == "__main__":
    image_path = "naip_rgb_ndvi_evi.tif"
    out_dir = "tiles/images"

    crop_geotiff_to_tiles(
        image_path=image_path,
        out_dir=out_dir,
        tile_size=500,
        stride=500
    )

如果希望增加样本数量,可以把 stride 改小,例如:

stride=250

这样会产生重叠裁剪样本。

七、如何制作语义分割 mask?

GEE 下载的是影像,不会自动生成精确的人工林 mask。
语义分割数据集中的 mask 一般需要人工标注或者半自动标注。

常见方法有三种:

1. 在 QGIS / ArcGIS 中绘制人工林矢量边界
2. 使用 LabelMe 标注 polygon
3. 使用已有林地矢量数据作为初始标签,再人工修正

如果你已经在 QGIS 中绘制了人工林 polygon,可以导出为 shapefile 或 geojson,然后使用 Python 栅格化成 mask。1. 安装依赖
pip install geopandas rasterio shapely numpy
2. 矢量标注转 mask

import rasterio
import geopandas as gpd
from rasterio.features import rasterize
import numpy as np

def vector_to_mask(image_path, vector_path, out_mask_path):
    """
    image_path: 参考影像路径
    vector_path: 人工林 polygon 标注文件,例如 shp 或 geojson
    out_mask_path: 输出 mask tif
    """

    with rasterio.open(image_path) as src:
        transform = src.transform
        height = src.height
        width = src.width
        crs = src.crs
        profile = src.profile.copy()

    gdf = gpd.read_file(vector_path)

    # 保证矢量和影像坐标系一致
    if gdf.crs != crs:
        gdf = gdf.to_crs(crs)

    shapes = [(geom, 1) for geom in gdf.geometry if geom is not None]

    mask = rasterize(
        shapes=shapes,
        out_shape=(height, width),
        transform=transform,
        fill=0,
        dtype=np.uint8
    )

    profile.update({
        "count": 1,
        "dtype": "uint8"
    })

    with rasterio.open(out_mask_path, "w", **profile) as dst:
        dst.write(mask, 1)

    print("mask 已保存:", out_mask_path)


if __name__ == "__main__":
    vector_to_mask(
        image_path="naip_rgb_ndvi_evi.tif",
        vector_path="plantation_label.geojson",
        out_mask_path="plantation_mask.tif"
    )

生成的 mask 中:

0:背景
1:人工林


八、同步裁剪 image 和 mask

语义分割任务要求 image 和 mask 严格一一对应。
因此,裁剪时必须使用相同的窗口。

import os
import rasterio
from rasterio.windows import Window
import numpy as np
from tqdm import tqdm

def crop_image_and_mask(
    image_path,
    mask_path,
    image_out_dir,
    mask_out_dir,
    tile_size=500,
    stride=500,
    min_valid_ratio=0.6
):
    os.makedirs(image_out_dir, exist_ok=True)
    os.makedirs(mask_out_dir, exist_ok=True)

    with rasterio.open(image_path) as img_src, rasterio.open(mask_path) as mask_src:
        width = img_src.width
        height = img_src.height

        img_profile = img_src.profile.copy()
        mask_profile = mask_src.profile.copy()

        tile_id = 0

        for y in tqdm(range(0, height - tile_size + 1, stride)):
            for x in range(0, width - tile_size + 1, stride):
                window = Window(x, y, tile_size, tile_size)

                image_tile = img_src.read(window=window)
                mask_tile = mask_src.read(1, window=window)

                valid_mask = np.any(image_tile != 0, axis=0)
                valid_ratio = valid_mask.sum() / (tile_size * tile_size)

                if valid_ratio < min_valid_ratio:
                    continue

                img_tile_profile = img_profile.copy()
                img_tile_profile.update({
                    "height": tile_size,
                    "width": tile_size,
                    "transform": img_src.window_transform(window)
                })

                mask_tile_profile = mask_profile.copy()
                mask_tile_profile.update({
                    "height": tile_size,
                    "width": tile_size,
                    "transform": mask_src.window_transform(window),
                    "count": 1,
                    "dtype": "uint8"
                })

                image_out_path = os.path.join(image_out_dir, f"image_{tile_id:05d}.tif")
                mask_out_path = os.path.join(mask_out_dir, f"mask_{tile_id:05d}.tif")

                with rasterio.open(image_out_path, "w", **img_tile_profile) as dst:
                    dst.write(image_tile)

                with rasterio.open(mask_out_path, "w", **mask_tile_profile) as dst:
                    dst.write(mask_tile.astype(np.uint8), 1)

                tile_id += 1

    print(f"同步裁剪完成,共生成 {tile_id} 对 image-mask 样本。")


if __name__ == "__main__":
    crop_image_and_mask(
        image_path="naip_rgb_ndvi_evi.tif",
        mask_path="plantation_mask.tif",
        image_out_dir="dataset/images",
        mask_out_dir="dataset/masks",
        tile_size=500,
        stride=500
    )

裁剪完成后,可以得到类似结构:

dataset/
├── images/
│   ├── image_00000.tif
│   ├── image_00001.tif
│   └── ...
└── masks/
    ├── mask_00000.tif
    ├── mask_00001.tif
    └── ...

九、将 GeoTIFF 转成 PNG/JPG 训练格式

有些语义分割代码习惯读取 PNG/JPG。
如果只使用 RGB 三通道,可以把 GeoTIFF 转为 PNG。

import os
import rasterio
import numpy as np
from PIL import Image
from tqdm import tqdm

def tif_to_png(tif_dir, png_dir):
    os.makedirs(png_dir, exist_ok=True)

    for name in tqdm(os.listdir(tif_dir)):
        if not name.endswith(".tif"):
            continue

        tif_path = os.path.join(tif_dir, name)

        with rasterio.open(tif_path) as src:
            img = src.read([1, 2, 3])  # 读取 RGB 三通道

        img = np.transpose(img, (1, 2, 0))

        # 如果是 0-1 浮点数据,转成 0-255
        if img.max() <= 1.0:
            img = img * 255

        img = np.clip(img, 0, 255).astype(np.uint8)

        out_name = name.replace(".tif", ".png")
        out_path = os.path.join(png_dir, out_name)

        Image.fromarray(img).save(out_path)

    print("转换完成。")


if __name__ == "__main__":
    tif_to_png(
        tif_dir="dataset/images",
        png_dir="dataset/images_png"
    )

mask 也可以转为单通道 PNG:

import os
import rasterio
import numpy as np
from PIL import Image
from tqdm import tqdm

def mask_tif_to_png(mask_tif_dir, mask_png_dir):
    os.makedirs(mask_png_dir, exist_ok=True)

    for name in tqdm(os.listdir(mask_tif_dir)):
        if not name.endswith(".tif"):
            continue

        tif_path = os.path.join(mask_tif_dir, name)

        with rasterio.open(tif_path) as src:
            mask = src.read(1)

        mask = mask.astype(np.uint8)

        out_name = name.replace(".tif", ".png")
        out_path = os.path.join(mask_png_dir, out_name)

        Image.fromarray(mask).save(out_path)

    print("mask 转换完成。")


if __name__ == "__main__":
    mask_tif_to_png(
        mask_tif_dir="dataset/masks",
        mask_png_dir="dataset/masks_png"
    )


十、划分训练集、验证集和测试集

import os
import random
import shutil

def split_dataset(
    image_dir,
    mask_dir,
    out_dir,
    train_ratio=0.7,
    val_ratio=0.2,
    seed=42
):
    random.seed(seed)

    names = sorted(os.listdir(image_dir))
    names = [n for n in names if n.endswith((".png", ".jpg", ".tif"))]

    random.shuffle(names)

    total = len(names)
    train_num = int(total * train_ratio)
    val_num = int(total * val_ratio)

    train_names = names[:train_num]
    val_names = names[train_num:train_num + val_num]
    test_names = names[train_num + val_num:]

    splits = {
        "train": train_names,
        "val": val_names,
        "test": test_names
    }

    for split, split_names in splits.items():
        os.makedirs(os.path.join(out_dir, split, "images"), exist_ok=True)
        os.makedirs(os.path.join(out_dir, split, "masks"), exist_ok=True)

        for name in split_names:
            image_path = os.path.join(image_dir, name)
            mask_name = name.replace("image", "mask")
            mask_path = os.path.join(mask_dir, mask_name)

            shutil.copy(image_path, os.path.join(out_dir, split, "images", name))
            shutil.copy(mask_path, os.path.join(out_dir, split, "masks", mask_name))

    print("数据集划分完成:")
    print("train:", len(train_names))
    print("val:", len(val_names))
    print("test:", len(test_names))


if __name__ == "__main__":
    split_dataset(
        image_dir="dataset/images_png",
        mask_dir="dataset/masks_png",
        out_dir="ArtificialForest-HRSegDataset"
    )


十一、下载我的公开数据集

如果只是想直接使用整理好的人工林语义分割数据集,可以从 Hugging Face 下载:

https://huggingface.co/datasets/jkhsdfja/ArtificialForest-HRSegDataset
方式一:使用 git clone
git lfs install
git clone https://huggingface.co/datasets/jkhsdfja/ArtificialForest-HRSegDataset
方式二:使用 huggingface_hub
pip install huggingface_hub
from huggingface_hub import snapshot_download

snapshot_download(
    repo_id="jkhsdfja/ArtificialForest-HRSegDataset",
    repo_type="dataset",
    local_dir="./ArtificialForest-HRSegDataset"
)
十二、数据集构建时需要注意的问题
1. 图像和 mask 必须严格对齐

遥感语义分割中,影像和标签只要错位几个像素,模型训练就会明显变差。
因此,裁剪 image 和 mask 时必须使用同一个窗口。

2. mask 不能使用双线性插值

如果对 mask resize,必须使用最近邻插值:

Image.NEAREST

或者:

cv2.INTER_NEAREST

mask 是类别索引图,不是普通灰度图,不能对类别编号做平均。

3. 尽量保留复杂场景

如果只保留边界清晰、背景简单的样本,模型在真实场景中容易失效。
论文数据集中保留了云影、山地、林地与非林地过渡区、建筑物和裸土干扰等复杂情况,用于提高模型泛化能力。

4. 植被指数要和原图保持同尺度

如果 RGB 图像是 500×500,NDVI 和 EVI 也要裁剪成同样大小。
否则在多模态融合时会出现空间错位。

十三、总结

本文记录了一个从 GEE 获取遥感图源并构建人工林语义分割数据集的完整流程。核心步骤包括:

1. 在 GEE 中选择研究区域
2. 下载 NAIP 或 Sentinel-2 影像
3. 计算 NDVI、EVI 等植被指数
4. 导出 GeoTIFF
5. 本地裁剪为 500×500 图像块
6. 制作人工林 mask
7. 划分训练集、验证集和测试集

对于人工林语义分割任务,RGB 图像可以提供纹理和空间结构信息,NDVI 和 EVI 可以提供植被生理信息,两者结合更有利于区分人工林、天然林、草地、裸土和建筑物等复杂背景。

公开数据集地址:

https://huggingface.co/datasets/jkhsdfja/ArtificialForest-HRSegDataset

后续我也会继续整理人工林语义分割模型训练、mask 可视化、mIoU 计算和多模态输入构建等相关内容。

Logo

AtomGit 是由开放原子开源基金会联合 CSDN 等生态伙伴共同推出的新一代开源与人工智能协作平台。平台坚持“开放、中立、公益”的理念,把代码托管、模型共享、数据集托管、智能体开发体验和算力服务整合在一起,为开发者提供从开发、训练到部署的一站式体验。

更多推荐