【记录】GEE下载MOD15A2H Fpar Lai数据 Arcpy裁剪及去除填充值

var roi = ee.FeatureCollection("projects/xx/xxx/xx");
// 1. 加载 ROI  

// 2. 读取 MOD15A2H 数据(2000-2024,每 8 天)
var dataset = ee.ImageCollection('MODIS/061/MOD15A2H')
  .filterBounds(roi)
  .filterDate('2022-10-15', '2022-10-17');
print('MODIS dataset loaded:', dataset);

// 3. 质量控制函数:将质量不好的像元赋值为 255
function maskQuality(image, bandName) {
  var qc = image.select('FparLai_QC');
  
  // 解析 bit 0,值为0表示质量好
  var goodQuality = qc.bitwiseAnd(1).eq(0);
  
  // 选择目标波段
  var band = image.select(bandName);
  
  // 对质量不好的像元(goodQuality为false)赋值为255,其余保持原值
  band = band.where(goodQuality.not(), 255);
  
  // 返回裁剪到ROI范围内并保留原影像属性的影像
  return band.clip(roi).copyProperties(image, image.propertyNames());
}

// 4. 处理 Fpar 和 Lai 波段
var fparCollection = dataset.map(function(image) {
  return maskQuality(image, 'Fpar_500m');
});
print('Processed Fpar collection:', fparCollection);

var laiCollection = dataset.map(function(image) {
  return maskQuality(image, 'Lai_500m');
});
print('Processed Lai collection:', laiCollection);

// 5. 按年份遍历影像,并逐个导出
var years = ee.List.sequence(2022, 2022);
print('Years to export:', years);

// 定义导出函数
function exportImagesByYear(year) {
  var startDate = ee.Date.fromYMD(year, 1, 1);
  var endDate = ee.Date.fromYMD(year, 12, 31);
  
  var fparYear = fparCollection.filterDate(startDate, endDate);
  var laiYear = laiCollection.filterDate(startDate, endDate);
  
  print('Processing year:', year);
  
  // 导出 Fpar 影像
  fparYear.evaluate(function(fparList) {
    fparList.features.forEach(function(feature) {
      var image = ee.Image(feature.id).select('Fpar_500m');
      var dateString = ee.Date(image.get('system:time_start')).format('YYYYMMdd').getInfo();
      
      print('Exporting Fpar for date:', dateString);
      
      Export.image.toDrive({
        image: image,
        description: 'Fpar_' + dateString,
        folder: 'GEE_Data_2022',
        scale: 500,
        region: roi.geometry(),
        crs: 'EPSG:4326',
        fileFormat: 'GeoTIFF',
        maxPixels: 1e13
      });
    });
  });

  // 导出 Lai 影像
  laiYear.evaluate(function(laiList) {
    laiList.features.forEach(function(feature) {
      var laiImage = ee.Image(feature.id).select('Lai_500m');
      var dateString = ee.Date(laiImage.get('system:time_start')).format('YYYYMMdd').getInfo();
      
      print('Exporting Lai for date:', dateString);
      
      Export.image.toDrive({
        image: laiImage,
        description: 'Lai_' + dateString,
        folder: 'GEE_Data_2022',
        scale: 500,
        region: roi.geometry(),
        crs: 'EPSG:4326',
        fileFormat: 'GeoTIFF',
        maxPixels: 1e13
      });
    });
  });
}

// 遍历每年并导出
years.evaluate(function(yearsList) {
  yearsList.forEach(function(y) {
    exportImagesByYear(y);
  });
});

// 6. 可视化最近一年的 Fpar 和 Lai
var latestYear = 2022;
var latestFpar = fparCollection.filterDate(ee.Date.fromYMD(latestYear, 1, 1), ee.Date.fromYMD(latestYear, 12, 31)).median();
var latestLai = laiCollection.filterDate(ee.Date.fromYMD(latestYear, 1, 1), ee.Date.fromYMD(latestYear, 12, 31)).median();

var visParams = { palette: ['blue', 'green', 'yellow', 'red'] };

Map.centerObject(roi, 6);
Map.addLayer(latestFpar, visParams, 'Fpar (' + latestYear + ')');
Map.addLayer(latestLai, visParams, 'Lai (' + latestYear + ')');
Map.addLayer(roi, {color: 'white'}, 'ROI');

质量控制的部分希望和大家多讨论 每次在gee上质量控制都有点心里没底 不知道怎么检验我质量控制到没

Arcpy批量处理(时间根据需要自己调)

import arcpy
from arcpy import env
from arcpy.sa import *
import os

# 设置输入和输出的根目录
input_root = r"F:/xxx"
output_root = r"F:/xxx/result"
mask_file = r"F:/xxxx/xxxx/xxxx/xxxx.shp"##裁剪的掩膜 研究区shp文件

# 检查并启用 Spatial Analyst 扩展
arcpy.CheckOutExtension("Spatial")

# 设置环境参数
env.overwriteOutput = True

# 定义年份范围
start_year = 2049
end_year = 2049

# 遍历每个年份文件夹
for year in range(start_year, end_year + 1):
    input_folder = os.path.join(input_root, str(year))
    output_folder = os.path.join(output_root, "Fpar_lai_{}".format(year))
    
    # 如果输入文件夹不存在,跳过
    if not os.path.exists(input_folder):
        print("输入文件夹 {} 不存在,已跳过。".format(input_folder))
        continue
    
    # 如果输出文件夹不存在,创建输出文件夹
    if not os.path.exists(output_folder):
        os.makedirs(output_folder)
        print("创建输出文件夹:{}".format(output_folder))
    
    # 设置当前工作空间
    env.workspace = input_folder
    
    # 列出当前年份文件夹中的所有栅格文件
    raster_list = arcpy.ListRasters()
    if not raster_list:
        print("文件夹 {} 中没有找到栅格文件,已跳过。".format(input_folder))
        continue
    
    print("正在处理年份:{},共 {} 个文件".format(year, len(raster_list)))

    # 遍历每个栅格文件
    for raster in raster_list:
        try:
            print("正在处理:{}".format(raster))
            
            # Step 1: 按掩膜提取
            extracted_raster = ExtractByMask(raster, mask_file)
            
            # Step 2: 执行 SetNull 操作(将值大于 100 的像素设置为 NoData)
            out_raster = SetNull(extracted_raster > 100, extracted_raster)
            
            # Step 3: 保存输出文件(保持原始文件名)
            output_path = os.path.join(output_folder, raster)
            out_raster.save(output_path)
            
            print("已保存:{}".format(output_path))
        
        except Exception as e:
            print("处理 {} 时出错:{}".format(raster, e))

print("所有文件处理完成!")

# 释放扩展许可
arcpy.CheckInExtension("Spatial")

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值