SC LSTR.py.zip

# -*- coding: utf-8 -*- # Landsat 8 LST Retrieval according to SC algorithm #LSE determinded according to empirical NDVI 0.7 and 0.05 from urllib2 import urlopen from arcpy import env import arcpy from dbfread import DBF import numpy as np import os,re,glob,math,json import pandas as pd arcpy.CheckOutExtension('spatial') env.overwriteOutput = True def RC(): B10 = arcpy.sa.Raster(''.join(glob.glob('*_B10.TIF'))) #L=MQ+A M is RADIANCE_MULT_BAND_10;A is RADIANCE_ADD_BAND_10 txt = ''.join(glob.glob('*.txt')) txt = open(txt).read() M = float(re.findall(r'RADIANCE_MULT_BAND_10 = (.*?)\n',txt)[0]) A = float(re.findall(r'RADIANCE_ADD_BAND_10 = (.*?)\n',txt)[0]) TOA = B10*M + A TOA.save('TOA.img') return 'TOA.img' #raster type def LSER(): ndvi = arcpy.RasterToNumPyArray('newNDVI.TIF') ndvi = list(ndvi.flatten()) ndvi.sort() NDVIVeg = 0.7 NDVISoil = 0.05 NDVI = arcpy.sa.Raster('newNDVI.TIF') #1 is Built-up; 2 is water; 3 is nature surface LUCC = arcpy.sa.Raster('LUCC.IMG') Pv = ((NDVI-NDVISoil)/(NDVIVeg-NDVISoil))**2 Pv = arcpy.sa.Con(NDVI>NDVIVeg,1,Pv) Pv = arcpy.sa.Con(NDVI<NDVISoil,0,Pv) outCon = arcpy.sa.Con(LUCC == 2,0.9908,0) #water outCon = arcpy.sa.Con((LUCC == 3)&(NDVI<NDVISoil),0.9722,outCon) #soil outCon = arcpy.sa.Con(NDVI>NDVIVeg,0.982,outCon)#vegetated outCon = arcpy.sa.Con((LUCC == 3)&(NDVI<NDVIVeg)&(NDVI>NDVISoil), 0.982*Pv+0.9722*(1-Pv)+0.01,outCon) #soil-veg-mixed outCon = arcpy.sa.Con((LUCC == 1)&(NDVI<NDVIVeg),0.9212,outCon) #built outCon = arcpy.sa.Con((LUCC == 1)&(NDVI<NDVIVeg)&(NDVI>NDVISoil), 0.982*Pv+0.9212*(1-Pv)+0.01,outCon) #built-veg-mixed outCon.save('LSE.img') return 'LSE.img' def CalcParam(waterImage,AT,UP,DOWN): water = arcpy.sa.Raster(waterImage) param1 = 0.04019 * water**2 + 0.02916 * water + 1.01523 param2 = -0.38333 * water**2 - 1.50294 * water + 0.20324 param3 = 0.00918 * water**2 + 1.36072 * water - 0.27514 param1 = arcpy.sa.Con(water<=3,param1,1.0/AT) param2 = arcpy.sa.Con(water<=3,param2,-DOWN-UP*1.0/AT) param3 = arcpy.sa.Con(water<=3,param3,DOWN) param1.save('param1.img') param2.save('param2.img') param3.save('param3.img') return 'param1.img','param2.img','param3.img' def LSTR(LSE,TOA,BT,param1,param2,param3): TOA = arcpy.sa.Raster(TOA) LSE = arcpy.sa.Raster(LSE) BT = arcpy.sa.Raster(BT) param1 = arcpy.sa.Raster(param1) param2 = arcpy.sa.Raster(param2) param3 = arcpy.sa.Raster(param3) br = 1324 r = BT**2/(br*TOA*1.0) LST = r*((param1*TOA+param2)/LSE+param3)+BT-BT**2/br LST = LST - 273 LST.save('scLST.img') return LST def BTR(TOA): TOA = arcpy.sa.Raster(TOA) BT = 1321.08/arcpy.sa.Ln(774.89/TOA+1) BT.save('BT.img') return 'BT.img' AI = pd.read_csv('e:/research/NDVI_LST/data/Atmospheric Index.csv') for i in os.listdir('e:/research/NDVI_LST/data/subLandsat8'): workspace=os.path.join('e:/research/NDVI_LST/data/subLandsat8',i) os.chdir(workspace) env.workspace = os.getcwd() index = AI.ix[AI.city==i[:-8],1:] AT = float(index['AT']); UP = float(index.UP); DOWN = float(index.DOWN) TOA = RC() BT = BTR(TOA) LSE = LSER() waterImage = 'water5.img' param1,param2,param3 = CalcParam(waterImage,AT,UP,DOWN) lst = LSTR(LSE,TOA,BT,param1,param2,param3) print i,lst.mean