批量NDVI时序

admin

[mw_shl_code=python,true]import os
from PIL import Image
import numpy as np
from osgeo import gdal
import glob
import cv2
list_tif = glob.glob('D:/BaiduNetdiskDownload/Shangyu/20210605/files/*.tif')
out_path = 'D:/NDVI/20210605/'
print(list_tif)
#记录异常值
count = 0

for tif in list_tif:
    in_ds = gdal.Open(tif)
    # 获取⽂件所在路径以及不带后缀的⽂件名
    (filepath, fullname) = os.path.split(tif)
    (prename, suffix) = os.path.splitext(fullname)
    if in_ds is None:
        print('Could not open the file ' + tif)
    else:
        # 读取波段数据
        red = in_ds.GetRasterBand(1).ReadAsArray()*0.0001
        nir = in_ds.GetRasterBand(4).ReadAsArray()*0.0001
        np.seterr(divide='ignore', invalid='ignore')
        ndvi = (nir - red) / (nir + red)
        # 将NAN转化为0值,如nodata转为0
        nan_index = np.isnan(ndvi)
        ndvi[nan_index] = 0
        ndvi = ndvi.astype(np.float32)
        #将计算好的NDVI保存为GeoTiff⽂件
        gtiff_driver = gdal.GetDriverByName('GTiff')
        # 批量处理需要注意⽂件名是变量，这⾥截取对应原始⽂件的不带后缀的⽂件名
        out_ds = gtiff_driver.Create(out_path + prename + 'NDVI.tif',
                         ndvi.shape[1], ndvi.shape[0], 1, gdal.GDT_Float32)
        # 将NDVI数据坐标投影设置为原始坐标投影
        out_ds.SetProjection(in_ds.GetProjection())
        out_ds.SetGeoTransform(in_ds.GetGeoTransform())
        out_band = out_ds.GetRasterBand(1)
        out_band.WriteArray(ndvi)
        out_band.FlushCache()
[/mw_shl_code]