# -*- 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
没有合适的资源?快使用搜索试试~ 我知道了~
SC LSTR.py.zip
共1个文件
py:1个
1.该资源内容由用户上传,如若侵权请联系客服进行举报
2.虚拟产品一经售出概不退款(资源遇到问题,请及时私信上传者)
2.虚拟产品一经售出概不退款(资源遇到问题,请及时私信上传者)
版权申诉
0 下载量 186 浏览量
2024-04-16
00:10:54
上传
评论
收藏 4KB ZIP 举报
温馨提示
SC LSTR.py.zip
资源推荐
资源详情
资源评论
收起资源包目录
SC LSTR.py.zip (1个子文件)
java
SC LSTR.py 4KB
共 1 条
- 1
资源评论
手把手教你学AI
- 粉丝: 7972
- 资源: 4768
下载权益
C知道特权
VIP文章
课程特权
开通VIP
上传资源 快速赚钱
- 我的内容管理 展开
- 我的资源 快来上传第一个资源
- 我的收益 登录查看自己的收益
- 我的积分 登录查看自己的积分
- 我的C币 登录后查看C币余额
- 我的收藏
- 我的下载
- 下载帮助
安全验证
文档复制为VIP权益,开通VIP直接复制
信息提交成功