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

dairui1240

已于 2022-09-09 16:45:17 修改

阅读量1k

点赞数

文章标签： servlet java 服务器

于 2022-09-09 15:57:10 首次发布

本文链接：https://blog.youkuaiyun.com/jnbvcyn/article/details/126784595

版权

这段代码展示了如何使用Python的arcpy模块处理遥感数据。它首先根据土地类型提取地表温度，然后对数据进行重采样，并计算不同年份不同土地利用的地表温度的统计指标，包括最小值、最大值、均值和方差。最后将结果写入文本文件。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

# -*- 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函数_忠言睿长的博客-优快云博客_arcgis setnull函数 https://blog.youkuaiyun.com/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")