extract_Index.zip

# -*- coding: utf-8 -*- import sys '''python import模块时， 是在sys.path里按顺序查找的。sys.path是一个列表，里面以字符串的形式存储了许多路径。使用db.py文件中的函数需要先将他的文件路径放到sys.path中''' sys.path.append(r'..\db') reload(sys) sys.setdefaultencoding('utf-8') import db import arcpy import os from arcpy.sa import * arcpy.CheckOutExtension("Spatial") # ArcGIS各种功能模块许可 __name__ = 'extract_Index' path = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) _inFolder = os.path.join(path, 'data', 'originalData') # 导入TIF文件路径 _outFolder = os.path.join(path, 'data', 'temp', 'index') # 输出数据文件夹 index_Math = 0.33 # SAVI指数土壤背景百分比 LS05 = "6" LS08 = "7" # 提取各项指数 def extractIndex(inFolder, outFolder): if not os.path.isdir(inFolder): print "输入的文件夹路径无效!" return print "遍历文件夹..." files = os.listdir(inFolder) # 返回文件夹下的所有文件名 arcpy.env.workspace = inFolder # 创建工作空间 arcpy.env.overwriteOutput = True # 覆盖相同文件 for list in files: if list.endswith(".tif"): inTifPath = os.path.join(inFolder, list) # .tif文件的绝对路径 bandCount = arcpy.GetRasterProperties_management(inTifPath, "BANDCOUNT").getOutput(0) # 波段总数为字符串形式 try: Band_RED = '' Band_NIR = '' Band_Blue = '' # 针对Landsat 5 和landsat 8 波段不一致 if bandCount == LS05: # 提取 红, 近红波段 Band_RED = arcpy.Raster(inTifPath + '/Band_3') # 读取红波段 Band_NIR = arcpy.Raster(inTifPath + '/Band_4') # 读取近红外波段 Band_Blue = arcpy.Raster(inTifPath + '/Band_1') # 读取蓝波段 elif bandCount == LS08: # 提取 红, 近红波段 Band_RED = arcpy.Raster(inTifPath + '/Band_4') # 读取红波段 Band_NIR = arcpy.Raster(inTifPath + '/Band_5') # 读取近红外波段 Band_Blue = arcpy.Raster(inTifPath + '/Band_2') # 读取蓝波段 # 计算各项指数 NDVI_Index(outFolder, list, bandCount, Band_NIR, Band_RED) DVI_Index(outFolder, list, bandCount, Band_NIR, Band_RED) SAVI_Index(outFolder, list, bandCount, Band_NIR, Band_RED) RVI_Index(outFolder, list, bandCount, Band_NIR, Band_RED) EVI_Index(outFolder, list, bandCount, Band_NIR, Band_RED, Band_Blue) print "提取", list, "所有指数成功!" except Exception, e: print "读取tif数据失败!", e # NDVI=(NIR-R)/(NIR+R) def NDVI_Index(outFolder, list, bandCount, Band_NIR, Band_RED): try: # landsat 系列不同命名不同 if bandCount == LS05: out_NDVI_Name = "LS_05_" + str(list[:4]) + "_NDVI" + ".tif" # Landsat 5 数据输出名称 elif bandCount == LS08: out_NDVI_Name = "LS_08_" + str(list[:4]) + "_NDVI" + ".tif" # Landsat 8 数据输出名称 out_NDVI_path = os.path.join(outFolder, "NDVI", out_NDVI_Name) # 输出文件路径 print out_NDVI_path NDVI = Float(Band_NIR - Band_RED) / Float(Band_NIR + Band_RED) # 计算NDVI outNDVI = Con(NDVI < -1.0, -1.0, Con(NDVI > 1.0, 1.0, NDVI)) # 剔除异常值 outNDVI.save(out_NDVI_path) # 导出栅格 arcpy.BuildPyramids_management(out_NDVI_path) # 构建栅格金字塔 arcpy.BatchCalculateStatistics_management(out_NDVI_path) # 构建直方图 print "提取", list, "NDVI成功!" except Exception, e: print "提取", list, "NDVI失败!", e return # DVI=NIR-R def DVI_Index(outFolder, list, bandCount, Band_NIR, Band_RED): try: # landsat 系列不同命名不同 if bandCount == LS05: out_DVI_Name = "LS_05_" + str(list[:4]) + "_DVI" + ".tif" # Landsat 5 数据输出名称 elif bandCount == LS08: out_DVI_Name = "LS_08_" + str(list[:4]) + "_DVI" + ".tif" # Landsat 8 数据输出名称 out_DVI_path = os.path.join(outFolder, "DVI", out_DVI_Name) # 输出文件路径 print out_DVI_path DVI = Float(Band_NIR) - Float(Band_RED) # 计算DVI 必须浮点型float outDVI = Con(DVI < 0, 0, Con(DVI > 16.0, 16.0, DVI)) # 剔除异常值 outDVI.save(out_DVI_path) # 导出栅格 arcpy.BuildPyramids_management(out_DVI_path) # 构建栅格金字塔 arcpy.BatchCalculateStatistics_management(out_DVI_path) # 构建直方图 print "提取", list, "DVI成功!" except Exception, e: print "提取", list, "DVI失败!", e return # SAVI=((NIR-R)/(NIR+R+L))(1+L) def SAVI_Index(outFolder, list, bandCount, Band_NIR, Band_RED): try: # landsat 系列不同命名不同 if bandCount == LS05: out_SAVI_Name = "LS_05_" + str(list[:4]) + "_SAVI" + ".tif" # Landsat 5 数据输出名称 elif bandCount == LS08: out_SAVI_Name = "LS_08_" + str(list[:4]) + "_SAVI" + ".tif" # Landsat 8 数据输出名称 out_SAVI_path = os.path.join(outFolder, "SAVI", out_SAVI_Name) # 输出文件路径 print out_SAVI_path # 计算SAVI 必须浮点型float SAVI = (Float(Band_NIR) - Float(Band_RED)) / (Float(Band_NIR) + Float(Band_RED) + index_Math) * (1 + index_Math) outSAVI = Con(SAVI < 0, 0.01, Con(SAVI > 1.3, 1.3, SAVI)) # 剔除异常值 outSAVI.save(out_SAVI_path) # 导出栅格 arcpy.BuildPyramids_management(out_SAVI_path) # 构建栅格金字塔 arcpy.BatchCalculateStatistics_management(out_SAVI_path) # 构建直方图 print "提取", list, "SAVI成功!" except Exception, e: print "提取", list, "SAVI失败!", e return # RVI=NIR/R def RVI_Index(outFolder, list, bandCount, Band_NIR, Band_RED): try: # landsat 系列不同命名不同 if bandCount == LS05: out_RVI_Name = "LS_05_" + str(list[:4]) + "_RVI" + ".tif" # Landsat 5 数据输出名称 elif bandCount == LS08: out_RVI_Name = "LS_08_" + str(list[:4]) + "_RVI" + ".tif" # Landsat 8 数据输出名称 out_RVI_path = os.path.join(outFolder, "RVI", out_RVI_Name) # 输出文件路径 print out_RVI_path RVI = Float(Band_NIR) / Float(Band_RED) # 计算RVI 必须浮点型float outRVI = Con(RVI < 0, 0, Con(RVI > 16.0, 16.0, RVI)) # 剔除异常值 outRVI.save(out_RVI_path) # 导出栅格 arcpy.BuildPyramids_management(out_RVI_path) # 构建栅格金字塔 arcpy.BatchCalculateStatistics_management(out_RVI_path) # 构建直方图 print "提取", list, "RVI成功!" except Exception, e: print "提取", list, "RVI失败!", e return # EVI 增强型植被指数 # EVI=2.5*(NIR-R)/(NIR+6.0*R-7.5*B+1) def EVI_Index(outFolder, list, bandCount, Band_NIR, Band_RED, Band_Blue): try: # landsat 系列不同命名不同 if bandCount == LS05: out_EVI_Name = "LS_05_" + str(list[:4]) + "_EVI" + ".tif" # Landsat 5 数据输出名称 elif bandCount == LS08: out_EVI_Name = "LS_08_" + str(list[:4]) + "_EVI" + ".tif" # Landsat 8 数据输出名称 out_EVI_path = os.path.join(outFolder, "EVI", out_EVI_Name) # 输出文件路径 print out_EVI_path # 计算EVI 必须浮点型float EVI = 2.5 * (Float(Band_NIR) - Float(Band_RED)) / (Float(Band_NIR) + 6 * Float(Band_RED)