农作物数据构建和统计（基于实测数据、GEE及python）

一、技术流程概览

本文通过 Google Earth Engine (GEE) 获取遥感影像数据，结合 Python 实现空间数据处理、分类模型训练及统计可视化，完整流程如下：

1. GEE 数据获取：下载预处理后的 Sentinel-2 影像

2. Python 数据对齐：融合矢量数据与栅格影像，生成标签矩阵

3. 随机森林分类：训练模型并预测农作物类别

4. 后处理与统计：平滑分类结果，计算面积并可视化

二、详细步骤与代码实现

1. GEE 数据获取

代码文件：gee.txt
功能：下载指定时间范围和区域的 Sentinel-2 影像，并进行云掩膜处理。

// 加载矢量边界  
var shp_asset = 'projects/ee-crop/assets/merged';  
var roi = ee.FeatureCollection(shp_asset).geometry();  

// 设置时间范围  
var start_date = '2025-03-01';  
var end_date = '2025-03-31';  

// 获取 Sentinel-2 数据并预处理  
var s2 = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')  
    .filterBounds(roi)  
    .filterDate(start_date, end_date)  
    .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 15))  
    .select(['B2', 'B3', 'B4', 'B8', 'B11', 'B12', 'QA60'])  
    .mosaic()  
    .clip(roi);  

// 云掩膜处理  
var s2_cloudmask = function(img) {  
    var qa60 = img.select('QA60');  
    var cloudBitMask = 1 << 10;  
    var cirrusBitMask = 1 << 11;  
    var mask = qa60.bitwiseAnd(cloudBitMask).eq(0)  
               .and(qa60.bitwiseAnd(cirrusBitMask).eq(0));  
    return img.updateMask(mask)  
              .select(['B2', 'B3', 'B4', 'B8', 'B11', 'B12']);  
};  
s2 = s2_cloudmask(s2);  

// 导出数据到 Google Drive  
Export.image.toDrive({  
    image: s2.toUint16(),  
    description: 'Sentinel2_5m_Bands',  
    folder: 'GEE_Exports',  
    region: roi,  
    scale: 10,  
    crs: 'EPSG:4326'  
});  

2. Python 数据对齐与标签生成

代码文件：GENERATE.py
功能：对齐栅格影像与矢量数据，生成训练标签。

import os  
import numpy as np  
import geopandas as gpd  
import rasterio  
from rasterio.mask import mask  

# 配置参数  
RASTER_PATH = "path/to/Sentinel2_5m_Bands.tif"  
VECTOR_PATHS = ["耕地数据1.shp", "耕地数据2.shp"]  
OUTPUT_DIR = "dataset_output"  
LABEL_MAPPING = {"小麦":1, "未耕":2, "油菜":3, "草地":4, "毛豆":5, "大棚":6, "玉米":7, "水稻":8}  

# 空间对齐函数  
def align_data(raster_path, vector_paths):  
    gdfs = []  
    for path in vector_paths:  
        gdf = gpd.read_file(path).to_crs(rasterio.open(raster_path).crs)  
        gdf['class_id'] = gdf["ZZZWMC"].map(LABEL_MAPPING).fillna(9)  # 未映射类别标记为9  
        gdfs.append(gdf)  
    combined_gdf = gpd.pd.concat(gdfs)  
    # 裁剪栅格并保存  
    with rasterio.open(raster_path) as src:  
        out_image, out_transform = mask(src, [combined_gdf.unary_union.__geo_interface__], crop=True)  
        out_meta = src.meta.copy()  
    # 保存对齐后的栅格  
    aligned_raster = os.path.join(OUTPUT_DIR, "aligned.tif")  
    with rasterio.open(aligned_raster, "w", **out_meta) as dest:  
        dest.write(out_image)  
    return aligned_raster, combined_gdf  

# 主流程  
if __name__ == "__main__":  
    os.makedirs(OUTPUT_DIR, exist_ok=True)  
    aligned_raster, gdf = align_data(RASTER_PATH, VECTOR_PATHS)  
    # 生成标签矩阵（详见原文件）  

生成数据

数据结构

3. 随机森林分类与后处理

代码文件：RF2_filter.py
功能：训练随机森林模型，预测农作物类别并平滑结果。

from sklearn.ensemble import RandomForestClassifier  
from skimage.morphology import opening, closing, square  

def postprocess(pred_matrix):  
    # 中值滤波 + 形态学操作  
    med_footprint = square(2)  
    morpho_footprint = square(2)  
    smoothed = median(pred_matrix, footprint=med_footprint)  
    smoothed = opening(smoothed, morpho_footprint)  
    smoothed = closing(smoothed, morpho_footprint)  
    return smoothed  

def main():  
    # 加载数据  
    X = np.load("dataset_output/image_cube.npy")  
    y = np.load("dataset_output/label_matrix.npy")  
    # 训练模型  
    rf = RandomForestClassifier(n_estimators=100, class_weight='balanced')  
    rf.fit(X_train, y_train)  
    # 全图预测  
    pred_matrix = rf.predict(X.reshape(-1, X.shape[-1]))  
    pred_matrix = pred_matrix.reshape(y.shape)  
    # 平滑处理  
    smoothed_pred = postprocess(pred_matrix)  
    # 保存结果（详见原文件）  

4. 统计与可视化

输出结果：

• 面积统计：生成 Excel 表格（area_results.xlsx）

• 分类报告：模型精度评估（classification_report.xlsx）

• 可视化图：原始影像、真实标签、预测结果对比图

最终结果