python NDVI栅格数据提取特定经纬度以及NDVI值

我有2000-2024年每个月的NDVI数据，想提取2022-2024年5-9月的，特定三个经纬度点的NDVI值

将提取结果写入excel

 

上代码：

#!usr/bin/env python
# -*- coding:utf-8 -*-
"""
@author: Suyue
@file: NDVItiqu.py
@time: 2025/08/08
@desc:
"""
import os
import rasterio
from rasterio.transform import from_origin
import pandas as pd
from pyproj import Transformer
import numpy as np
import warnings

warnings.filterwarnings("ignore")

# 配置参数
input_folder = "G:/标准地图/NDVI/内蒙古自治区"  # 替换为你的文件夹路径
output_excel = "G:/标准地图/NDVI/内蒙古自治区/NDVI_extracted.xlsx"  # 输出Excel文件名

# 定义要提取的经纬度点 (lat, lon)
target_points = [
    (38.18, 107.48),
    (39.81, 108.71),
    (38.6, 108.83)
]


def get_tif_files(folder, start_year=2002, end_year=2022):
    """获取2002-2022年5-9月的TIFF文件"""
    files = []
    for year in range(start_year, end_year + 1):
        for month in range(5, 10):  # 5-9月
            # 构建文件名模式，如NDVI_200205.tif
            pattern = f"NDVI_{year}{month:02d}.tif"
            # 查找匹配的文件
            for f in os.listdir(folder):
                if f.startswith(pattern):
                    files.append(os.path.join(folder, f))
                    break
    return sorted(files)


def latlon_to_pixel(lat, lon, transform, crs):
    """将经纬度坐标转换为像素坐标"""
    # 创建坐标转换器 (WGS84 to raster CRS)
    transformer = Transformer.from_crs("EPSG:4326", crs, always_xy=True)
    # 转换坐标
    x, y = transformer.transform(lon, lat)
    # 将地理坐标转换为像素坐标
    col, row = ~transform * (x, y)
    return int(row), int(col)


def extract_ndvi_value(tif_path, points):
    """从TIFF文件中提取指定点的NDVI值（修正版）"""
    results = []
    with rasterio.open(tif_path) as src:
        # 读取NDVI数据 (假设是第一个波段)
        ndvi_data = src.read(1)

        # 获取年月信息从文件名
        filename = os.path.basename(tif_path)
        year = int(filename.split('_')[1][:4])
        month = int(filename.split('_')[1][4:6])

        for lat, lon in points:
            try:
                # 转换经纬度到像素坐标
                row, col = latlon_to_pixel(lat, lon, src.transform, src.crs)

                # 检查坐标是否在图像范围内
                if 0 <= row < src.height and 0 <= col < src.width:
                    ndvi = ndvi_data[row, col]

                    # 检查是否是无效值 (MODIS NDVI通常用-3000表示无效)
                    if ndvi == src.nodata:
                        ndvi = None
                    else:
                        # 先检查原始值范围
                        if -1 <= ndvi <= 1:  # 如果已经在合理范围内
                            pass  # 不需要缩放
                        elif -3000 <= ndvi <= 10000:  # 如果是原始MODIS值
                            ndvi = ndvi * 0.0001  # 应用标准缩放因子
                        else:
                            print(f"警告: 文件 {filename} 中的值 {ndvi} 超出预期范围")
                            ndvi = None
                else:
                    ndvi = None

                results.append({
                    "year": year,
                    "month": month,
                    "lat": lat,
                    "lon": lon,
                    "NDVI": ndvi
                })

            except Exception as e:
                print(f"Error processing point ({lat}, {lon}) in {filename}: {e}")
                results.append({
                    "year": year,
                    "month": month,
                    "lat": lat,
                    "lon": lon,
                    "NDVI": None
                })

    return results


def main():
    # 获取所有符合条件的TIFF文件
    tif_files = get_tif_files(input_folder)
    print(f"找到 {len(tif_files)} 个符合条件的文件")

    # 创建空的DataFrame存储结果
    df = pd.DataFrame(columns=["year", "month", "lat", "lon", "NDVI"])

    # 处理每个文件
    for tif_file in tif_files:
        print(f"正在处理: {os.path.basename(tif_file)}")
        results = extract_ndvi_value(tif_file, target_points)
        df = pd.concat([df, pd.DataFrame(results)], ignore_index=True)

    # 保存到Excel
    df.to_excel(output_excel, index=False)
    print(f"数据已保存到 {output_excel}")

    # 打印一些样本数据检查
    print("\n样本数据预览:")
    print(df.head(10))


if __name__ == "__main__":
    main()