extract_Index.zip

  • AI_Stu
    了解作者
  • Python
    开发工具
  • 2KB
    文件大小
  • zip
    文件格式
  • 0
    收藏次数
  • 1 积分
    下载积分
  • 0
    下载次数
  • 2021-04-27 18:55
    上传日期
批处理计算一个文件夹下所有.TIF文件的NDVI、DVI、SAVI、RVI、EVI指数
extract_Index.zip
  • extract_Index.py
    8.3KB
内容介绍
# -*- 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)
评论
    相关推荐
    • arcpy.rar
      可以进行批量裁剪TIF格式文件,适用于大量数据预处理
    • rest-query-to-fgdb-arcpy
      休息查询到 fgdb-arcpy 这是一个脚本,用于下载通过 ESRI 其余服务公开的数据并将其保存为地理数据库要素类。 如果 REST 服务有 n 个特征的限制,脚本将一次下载 n/2 个特征,然后合并结果。 此脚本使用 ArcPy 并...
    • arcpy共举包
      arcpy ArcGIS 10 引入了 ArcPy,这是一个 Python 站点包,它涵盖并进一步加强了 ArcGIS 9.2 中所采用的 arcgisscripting 模块的功能。ArcPy 提供了一种用于开发 Python 脚本的功能丰富的动态环境,同时提供每个函数...
    • arcpy 工具包.zip
      arcpy ArcGIS 10 引入了 ArcPy,这是一个 Python 站点包,它涵盖并进一步加强了 ArcGIS 9.2 中所采用的 arcgisscripting 模块的功能。ArcPy 提供了一种用于开发 Python 脚本的功能丰富的动态环境,同时提供每个函数...
    • ArcPy and ArcGIS - Second Edition
      The arcpy.mapping Moduleb'Chapter 6: The arcpy.mapping Module'b'Using Arc Py with map documents \xc3\x82\xc2\xa0'b'Summary'7: Advanced Analysis Topicsb'Chapter 7: Advanced Analysis Topics'b'Using ...
    • ArcGIS开发arcpy教程
      arcgis数据处理过程中使用到的arcpy脚本,基础教程。有助于gis从业人员开发arcpy脚本便捷处理数据。提供基础与案例应用结合方式,讲解arcpy知识。 1.python基础 2.基本图形创建 3.shapefile数据操作 4.常见数据txt....
    • geopublisher:使用 Esri 的 Arcpy 发布 GIS 数据的工具
      使用 Esri 的 Arcpy 发布 GIS 数据的工具。 免费软件:Apache License 2.0 文档: : 。 特征 向/从 shapefile、文件地理数据库、企业级地理数据库 (SDE) 发布 可以通过将已发布要素类(如有必要)导出到文件名中带...
    • arcpy 工具包
      arcpy ArcGIS 10 引入了 ArcPy,这是一个 Python 站点包,它涵盖并进一步加强了 ArcGIS 9.2 中所采用的 arcgisscripting 模块的功能。ArcPy 提供了一种用于开发 Python 脚本的功能丰富的动态环境,同时提供每个函数...
    • arcpy-extensions:ArcGIS python界面arcpy的扩展和帮助程序
      arcpy扩展 正在进行的python程序包,其中包含arcpy的辅助函数和类。 当前的开发是由其他项目的需求驱动的,而不是统一的开发重点,但是我愿意帮助,指导和/或提出要求。 安装 要安装此软件包,请使用pip pip ...
    • describer:一个简单的python模块,帮助在arcpy中使用Describe函数
      描述者 描述器试图通过根据输入数据类型显示您有权访问的所有有效... Describer 还会为您进行描述,因此在将输入输入到类中之前无需使用 arcpy 来描述输入。 例子 import describer desc = ' \\ Example.gdb' D = desc