Google Earth Engine导出的MODIS月度NDVI栅格全为NaN值，但可视化正常的问题求助

大家好，我最近在使用Google Earth Engine（GEE）处理MODIS MOD13Q1数据集，想要导出特定火灾多边形区域内的月度中位数NDVI栅格，但遇到了一个奇怪的问题：

• 在GEE控制台里，NDVI的可视化完全正常，数值也符合预期；
• 但导出的GeoTIFF文件却全是NaN值，而且日志里没有任何报错信息；
• 脚本能成功启动导出任务，可不管用QGIS还是Python的rasterio打开，文件里都找不到有效的NDVI数据。

我已经尝试过的操作 #

1. 在GEE控制台验证NDVI数据（成功） #

为了确认目标区域存在MODIS NDVI数据，我先写了一段JS代码展示2020年8月的中位数NDVI，结果显示完全正常：

var fireRegion = ee.Geometry.Polygon([
  [[-121.974, 37.634], [-121.974, 36.972], [-120.700, 36.972], [-120.700, 37.634]]
]);

var modis = ee.ImageCollection("MODIS/061/MOD13Q1")
    .select("NDVI")
    .filterBounds(fireRegion)
    .filterDate("2020-08-01", "2020-08-31")
    .median()
    .multiply(0.0001)  // 正确缩放NDVI至标准范围
    .clip(fireRegion);

var ndviVis = {min: -0.2, max: 1.0, palette: ["red", "yellow", "green"]};

Map.centerObject(fireRegion, 8);
Map.addLayer(modis, ndviVis, "MODIS NDVI (August 2020)");
print("NDVI Image:", modis);

这段代码运行后，NDVI在GEE地图上渲染完全正常，说明数据本身是存在的。

2. Python批量导出脚本（失败） #

接下来我写了Python脚本，想要批量导出2019-2022年每个月的NDVI，代码如下：

import ee
import geopandas as gpd
import logging
from datetime import datetime, timedelta

# 设置日志输出
logging.basicConfig(level=logging.INFO, format="%(asctime)s - %(levelname)s - %(message)s")

# 初始化Earth Engine
try:
    ee.Initialize(project="sound-of-resiliency")
except Exception as e:
    ee.Authenticate()
    ee.Initialize(project="sound-of-resiliency")

# 加载火灾多边形shapefile
shapefile = gpd.read_file("fire_poly.shp")

# 确保shapefile转换为WGS84投影（EPSG:4326）
if shapefile.crs != "EPSG:4326":
    shapefile = shapefile.to_crs("EPSG:4326")

# 将shapefile转换为GEE几何对象
fire_region_geom = shapefile.geometry.union_all()
fire_region = ee.Geometry(fire_region_geom.__geo_interface__)

# 加载MODIS NDVI数据集
modis = ee.ImageCollection("MODIS/061/MOD13Q1").select("NDVI").filterBounds(fire_region)

def get_monthly_ndvi(year, month):
    """计算指定月份的中位数NDVI"""
    start = datetime(year, month, 1)
    end = (start + timedelta(days=32)).replace(day=1)  # 取下个月第一天作为结束日期

    monthly_collection = modis.filterDate(start.strftime("%Y-%m-%d"), end.strftime("%Y-%m-%d"))

    # 调试：打印当月找到的影像数量
    image_count = monthly_collection.size().getInfo()
    logging.info(f"{year}-{month} 找到的MODIS影像数量: {image_count}")

    if image_count == 0:
        logging.warning(f"{year}-{month} 没有找到MODIS NDVI影像，跳过该月")
        return None  

    monthly_ndvi = (
        monthly_collection.median()
        .multiply(0.0001)  # 转换为NDVI标准范围(-1到1)
        .clip(fire_region)
    )

    return monthly_ndvi

def export_ndvi_raster(image, year, month):
    """将NDVI影像导出到Google Drive"""
    if image is None:
        logging.warning(f"跳过 {year}-{str(month).zfill(2)} 的导出（无NDVI数据）")
        return  

    # 重投影到EPSG:4326，分辨率250m
    image = image.reproject(crs="EPSG:4326", scale=250)

    # 计算区域内NDVI均值，验证是否有有效像素
    stats = image.reduceRegion(
        reducer=ee.Reducer.mean(),
        geometry=fire_region,
        scale=250,
        maxPixels=1e13
    )

    ndvi_value = stats.get("NDVI")

    if ndvi_value is not None:
        task = ee.batch.Export.image.toDrive(
            image=image,
            description=f"NDVI_{year}_{str(month).zfill(2)}",
            folder="GEE_Exports",
            fileNamePrefix=f"NDVI_{year}_{str(month).zfill(2)}",
            scale=250,  
            region=fire_region,
            fileFormat="GeoTIFF",
            maxPixels=1e13
        )
        task.start()
        logging.info(f"{year}-{str(month).zfill(2)} 的NDVI导出任务已启动")
    else:
        logging.warning(f"跳过 {year}-{str(month).zfill(2)} 的导出（无有效NDVI像素）")

# 遍历2019-2022年每个月
for year in range(2019, 2023):  
    for month in range(1, 13):  
        logging.info(f"正在处理 {year}-{str(month).zfill(2)} 的NDVI")
        monthly_ndvi = get_monthly_ndvi(year, month)
        export_ndvi_raster(monthly_ndvi, year, month)

logging.info("所有NDVI导出任务已启动，请查看GEE任务面板")

3. 检查导出文件 #

从Google Drive下载导出的.tif文件后，我用Python的rasterio做了分析：

import rasterio
import numpy as np

with rasterio.open("NDVI_2022_12.tif") as src:
    data = src.read(1)
    print(f"最小值: {np.nanmin(data)}, 最大值: {np.nanmax(data)}, 平均值: {np.nanmean(data)}")

结果输出全是NaN：

最小值: nan, 最大值: nan, 平均值: nan

我已经排查过的点 #

• 火灾多边形（fire_region）的定义是正确的，因为控制台可视化正常；
• 每个月都能找到对应的MODIS NDVI影像（通过size().getInfo()验证）；
• NDVI已经做了正确的缩放（乘以0.0001）；
• 导出前已经将影像重投影到EPSG:4326；
• 导出任务在GEE的“任务”面板里显示正常完成。

但导出的栅格还是全是NaN，我现在有几个疑问想请教大家：

1. 为什么NDVI在GEE控制台渲染正常，但用.median()计算后导出就变成全NaN？
2. 有没有更可靠的方法从MOD13Q1计算并导出月度NDVI？
3. .median()操作会不会导致像素缺失的问题？
4. 导出时的scale=250需要调整吗？还是MODIS的投影存在兼容性问题？

备注：内容来源于stack exchange，提问作者Derek Corcoran