############################################################
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @Time : 2017/12/20 15:28
# @Author : Zhang Qi
# @File : tiffcal.py
############################################################
from osgeo import gdal
import sys
import os
import numpy as np
arg=sys.argv #获得命令行参数
filenl=arg[1:len(arg)] #arg[0]为文件名,所以去除
print(str(filenl[0]))
data_all=[]
# 读文件:
dataset = gdal.Open(str(filenl[0]))
width = dataset.RasterXSize # 栅格矩阵的列数
height = dataset.RasterYSize # 栅格矩阵的行数
bands = dataset.RasterCount # 波段数
geotrans = dataset.GetGeoTransform() # 获取仿射矩阵信息
proj = dataset.GetProjection() # 获取投影信息
for i,a in enumerate(filenl):
src_ds = gdal.Open(str(a))
# im_width = src_ds.RasterXSize # 栅格矩阵的列数
# im_height = src_ds.RasterYSize # 栅格矩阵的行数
# im_bands = src_ds.RasterCount # 波段数
# im_geotrans = src_ds.GetGeoTransform() # 获取仿射矩阵信息
# im_proj = src_ds.GetProjection() # 获取投影信息
if src_ds is None:
print('Unable to open INPUT.tif')
sys.exit(1)
print("[ RASTER BAND COUNT ]: ", src_ds.RasterCount)
if src_ds.RasterCount==1: #一层band的
im_width = src_ds.RasterXSize # 栅格矩阵的列数
im_height = src_ds.RasterYSize # 栅格矩阵的行数
data=list(src_ds.ReadAsArray(0,0,im_width,im_height))
data_all.append(data)
else:
for band in range(src_ds.RasterCount):
band += 1
print("[ GETTING BAND ]: ", band)
srcband = src_ds.GetRasterBand(band)
if srcband is None:
continue
data=srcband.ReadAsArray() #存在data中
data_all=np.asarray(data_all) #将tiff读入保存。
def lst(A,B,C): # A:T4 B:T5 C:ndvi
Pv = (C - np.min(C)) / (np.max(C) - np.min(C))
Tv = A + 2.6 * (A - B) - 2.4
Tbs = A + 2.1 * (A - B) - 3.1
Ts= Pv*Tv+(1-Pv)*Tbs
return Ts
#写文件:
def writeTiff(im_data, im_width, im_height, im_bands, im_geotrans, im_proj, path):
if 'int8' in im_data.dtype.name:
datatype = gdal.GDT_Byte
elif 'int16' in im_data.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float32
if len(im_data.shape) == 3:
im_bands, im_height, im_width = im_data.shape
else:
im_bands, (im_height, im_width) = 1, im_data.shape
# 创建文件
driver = gdal.GetDriverByName("GTiff")
dataset = driver.Create(path, im_width, im_height, im_bands, datatype)
if (dataset != None):
dataset.SetGeoTransform(im_geotrans) # 写入仿射变换参数
dataset.SetProjection(im_proj) # 写入投影
dataset.GetRasterBand(im_bands).WriteArray(im_data)
path=os.path.split(os.path.realpath(__file__))[0]+ r'\\result.tif'
print(path)
lstout=lst(data_all[0],data_all[1],data_all[2])
writeTiff(lstout,width,height,bands,geotrans,proj, path)
命令行调用方式:
python *.py T4.tiff T5.tiff ndvi.tiff