批量NDVI时序
import osfrom 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 = 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, ndvi.shape, 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()
页:
[1]