Python批量计算NDVI
PYTHON批量计算NDVI
做了少量修改，剔除了异常值，执⾏代价时需要更换影像对应波段及⽂件存储位置
import os
import numpy as np
from osgeo import gdal
import glob
import time
list_tif = glob.glob("F://2020TreeClassification//NDVIPY//MASKGF//*.tif")
out_path ='F://2020TreeClassification//NDVIPY//NDVIGF'
start = time.perf_counter()
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:
# 将MODIS原始数据类型转化为反射率
red = in_ds.GetRasterBand(3).ReadAsArray()*0.0001
nir = in_ds.GetRasterBand(4).ReadAsArray()*0.0001
np.seterr(divide='ignore', invalid='ignore')
ndvi =(nir - red)/(nir + red)
ndvi[ndvi >1]=0# 过滤异常值
ndvi[ndvi <0]=0# 过滤异常值
# 将NAN转化为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()
elapsed =(time.perf_counter()- start)
print("影像输出完毕")
print("计算ndvi耗时:", elapsed)
如果原影像是.dat⽂件
import os
import numpy as np
from osgeo import gdal
import glob
import time
list_dat = glob.glob("F://2020TreeClassification//NDVIPY//MASKGF//*.dat")
out_path ='F://2020TreeClassification//NDVIPY//NDVI//NDVI_'
start = time.perf_counter()
for dat in list_dat:
in_ds = gdal.Open(dat)
# 获取⽂件所在路径以及不带后缀的⽂件名
(filepath, fullname)= os.path.split(dat)
(prename, suffix)= os.path.splitext(fullname)
if in_ds is None:
print('Could not open the file '+ dat)
else:
# 将MODIS原始数据类型转化为反射率
red = in_ds.GetRasterBand(3).ReadAsArray()*0.0001
nir = in_ds.GetRasterBand(4).ReadAsArray()*0.0001
np.seterr(divide='ignore', invalid='ignore')
ndvi =(nir - red)/(nir + red)
ndvi[ndvi >1]=0# 过滤异常值
python怎么读取dat文件ndvi[ndvi <0]=0# 过滤异常值
# 将NAN转化为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()
elapsed =(time.perf_counter()- start)
print("影像输出完毕")
print("计算ndvi耗时:", elapsed)

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。