根据不同土地类型提取对应地表温度，并统计对应最小值、最大值、均值和方差

# -*- coding: utf-8 -*-
"""
Created on Tue Sep  6 09:35:57 2022

@author: DR
"""

##根据不同土地类型提取对应地表温度，并统计对应最小值、最大值、均值和方差

import glob
import os
import shutil
from arcpy.sa import *
from arcpy.sa import Con

##提取单个土地类型
path2 = r"I:\DR\test\2lc_HF"
path3 = "I:/DR/test"
if not os.path.exists(path3):##path5是否存在，不存在先建立。
    os.mkdir(path3)

rasters2 = glob.glob(os.path.join(path2, "*.tif"))
for ras2 in rasters2:
    nameT2 = ras2[18:22] + ras2[27:31] + "_1Cro.tif"#cl：croplands农田
    print(nameT2)
    newFileDirs2 = path3 + "/" + ras2[18:23] + ras2[27:31]
    # print(newFileDirs2)
    if not os.path.exists(newFileDirs2):
        os.mkdir(newFileDirs2)
    outname2 = os.path.join(newFileDirs2, nameT2)
    attExtract2 = ExtractByAttributes(ras2, "VALUE = 1") 
    attExtract2.save(outname2)
   
    nameT3 = ras2[18:22] + ras2[27:31] + "_2For.tif"#forest森林
    print(nameT3)
    outname3 = os.path.join(newFileDirs2, nameT3)
    attExtract3 = ExtractByAttributes(ras2, "VALUE = 2") 
    attExtract3.save(outname3)
    
    nameT4 = ras2[18:22] + ras2[27:31] + "_5Wat.tif"#water水体
    print(nameT4)
    outname4 = os.path.join(newFileDirs2, nameT4)
    attExtract4 = ExtractByAttributes(ras2, "VALUE = 5") 
    attExtract4.save(outname4)
    
    nameT5 = ras2[18:22] + ras2[27:31] + "_8Imp.tif"#不透水面
    print(nameT5)
    outname5 = os.path.join(newFileDirs2, nameT5)
    attExtract5 = ExtractByAttributes(ras2, "VALUE = 8") 
    attExtract5.save(outname5)
    print "Fire in the hole!"
    
##path1重采样1km
path1 = r"I:\DR\test"
path2 = r"I:\DR\CLCD_Rec"
path_list3 = os.listdir(path1)
# print(path_list3)
for i in range(len(path_list3)):
    inpath3 = path1 + "\\" + path_list3[i]
    # print(inpath3)
    rasters3 = glob.glob(os.path.join(inpath3, "*.tif"))
    # print(rasters4)
    for ras3 in rasters3:
        nameT3 = os.path.basename(ras3).split(".")[0] + "_Res.tif"
        # print(nameT3)
        outname3 = os.path.join(path2, nameT3)
        # print(outname3)
        arcpy.Resample_management(ras3, outname3, 0.010942416, "NEAREST")
        print os.path.basename(ras3) + " OK!"
    
##path1*path2,分别提取不同年份不同土地利用对应的地表温度
path1 = r"I:\DR\DAY_Avg_Year"
path2 = r"I:\DR\CLCD_Rec"
path3 = "I:/DR/LST_CLCD/"

rasters1 = glob.glob(os.path.join(path1, "*.tif"))
for i in range(2003,2022):
    for ras1 in rasters1:
        if int(ras1[19:23])==i:
            rasters2 = glob.glob(os.path.join(path2, "*.tif"))
            for ras2 in rasters2:
                if int(ras2[19:23])==i:
                    outname = path3 + "LST_" + ras2[15:]
                    if int(ras2[24])==1:
                        arcpy.gp.RasterCalculator_sa("'" + ras1 + "'*'" + ras2 + "'/1", outname)
                    ##(75条消息) arcPython细节汇总-setll函数和地图计算Raster Calculator函数_忠言睿长的博客-CSDN博客_arcgis setnull函数  https://blog.csdn.net/liyanzhong/article/details/46270757
                    elif int(ras2[24])==2:
                        arcpy.gp.RasterCalculator_sa("'" + ras1 + "'*'" + ras2 + "'/2", outname)#缩小2倍
                    elif int(ras2[24])==5:
                        arcpy.gp.RasterCalculator_sa("'" + ras1 + "'*'" + ras2 + "'/5", outname)#缩小5倍
                    else:
                        arcpy.gp.RasterCalculator_sa("'" + ras1 + "'*'" + ras2 + "'/8", outname)#缩小8倍
                        
##统计不同年份不同土地利用对应的地表温度的最小值，最大值、均值和方差
outputTxtFile =  r"I:\DR\result1.txt"
arcpy.env.workspace = r"I:\DR\DAY_Avg_Year"
rasters = arcpy.ListRasters("*","tif")
# print(rasters)
for raster in rasters:
    minResult = arcpy.GetRasterProperties_management(raster, "MINIMUM")
    maxResult = arcpy.GetRasterProperties_management(raster, "MAXIMUM")
    meanResult = arcpy.GetRasterProperties_management(raster, "MEAN")
    stdResult = arcpy.GetRasterProperties_management(raster, "STD")
    img_SR = arcpy.Describe(raster).spatialReference.name
    minRes = minResult.getOutput(0)
    maxRes = maxResult.getOutput(0)
    meanRes = meanResult.getOutput(0)
    stdRes = stdResult.getOutput(0)
    # minRes=str(minResult)
    # maxRes = str(maxResult)
    # meanRes = str(meanResult)
 
 
    f = open(outputTxtFile, 'a+')#打开txt文件，并追加信息
    f.write(raster + "," + minRes + "," + maxRes + "," + meanRes + "," + stdRes + "\n")
 
 
print ("finish")