Arcpy进行地图代数栅格运算

Python利用Arcpy进行栅格运算ArcGIS栅格计算器介绍Arcpy对栅格文件的读取、创建和写入Arcpy进行栅格的数学分析和逻辑运算归纳总结环境;ArcGIS、Python 2.7ArcGIS栅格计算器介绍ArcGIS中的工具箱有栅格计算器工具，它提供的时地图代数功能，可以实现很多复杂的栅格代数运算处理，但栅格计算器工具专门用于 ArcGIS for Desktop 应用程序（仅作为 GP

jing_zhong

13464人浏览 · 2020-10-07 21:39:53

jing_zhong · 2020-10-07 21:39:53 发布

Python利用Arcpy进行栅格运算

ArcGIS栅格计算器介绍
Arcpy对栅格文件的读取、创建和写入
Arcpy进行栅格数学分析的三角函数和逻辑运算归纳总结
实例及运行结果

环境;ArcGIS、Python 2.7

ArcGIS栅格计算器介绍

ArcGIS中的工具箱有栅格计算器工具，它提供的时地图代数功能，可以实现很多复杂的栅格代数运算处理，但栅格计算器工具专门用于 ArcGIS for Desktop 应用程序（仅作为 GP 工具对话框）或模型构建器。它不适用于脚本的编写，而且也不能用于 ArcPy Spatial Analyst 模块。

在这里插入图片描述
为此，通过查看ArcGIS 工具箱中的空间分析模块（Spatial Analyst简称sa），发现其中有很多数学分析运算（我们主要用到的多为三角函数和逻辑运算）（见下图）

因此，为了更好地利用Arcpy中的空间分析工具来进行栅格数据的地图代数运算，需要首先在Python的py文件中导入ArcPy和Spatial Analysis模块，见下面的代码：

#-*- coding: UTF-8 -*-
import arcpy
from arcpy import env
from arcpy.sa import *

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

Arcpy对栅格文件的读取、创建和写入

Arcpy读取栅格：

inRaster = arcpy.Raster("D:\\test\\myinputraster.tif")#读取tif文件

Arcpy可以创建常量栅格、正态栅格和随机栅格：

outConstRaster = CreateConstantRaster(12.7, "FLOAT", 2, Extent(0, 0, 250, 250))#常量栅格创建
outNormalRaster = CreateNormalRaster(2, Extent(0, 0, 150, 150))#正态栅格创建
outRandRaster = CreateRandomRaster(100, 2, Extent(0, 0, 150, 150))#随机栅格创建

在这里插入图片描述
Arcpy保存栅格：

OutputRaster.save("D:\\test\\myoutputraster.tif")#保存tif文件

Arcpy进行栅格数学分析的三角函数和逻辑运算归纳总结

三角函数运算	对应代码	作用
Sin	outSinRaster = Sin(inRaster)	计算栅格中像元的正弦
Cos	outCosRaster = Cos(inRaster)	计算栅格中像元的余弦
Tan	outTanRaster = Tan(inRaster)	计算栅格中像元的正切
SinH	outSinHRaster = SinH(inRaster)	计算栅格中像元的双曲正弦
CosH	outCosHRaster = CosH(inRaster)	计算栅格中像元的双曲余弦
TanH	outTanHRaster = TanH(inRaster)	计算栅格中像元的双曲正切
ASin	outASinRaster = ASin(inRaster)	计算栅格中各像元的反正弦
ACos	outACosRaster = ACos(inRaster)	计算栅格中像元的反余弦
ATan	outATanRaster = ATan(inRaster)	计算栅格中各像元的反正切
ATan2	outATan2Raster = ATan2(inRaster1, inRaster2)	计算栅格中各像元的反正切(基于 x、y)
ASinH	outASinHRaster = ASinH(inRaster)	计算栅格中各像元的反双曲正弦
ACosH	outACosHRaster = ACosH(inRaster)	计算栅格中各像元的反双曲余弦
ATanH	outATanHRaster = ATanH(inRaster)	计算栅格中各像元的反双曲正切

在这里插入图片描述

逻辑运算	对应代码	作用
加	outPlus = Plus(inRaster1, inRaster2)	逐个像元地将两个栅格的值相加（求和）
减	outMinus = Minus(inRaster1, inRaster2)	逐个像元地从第一个输入栅格的值中减去第二个输入栅格的值
乘	outTimes = Times(inRaster, inConstant)	将两个栅格的值逐个像元地相乘
除	outDivide = Divide(inRaster01, inRaster02)	将两个栅格的值逐个像元地相除
平方	outSquare = Square(inRaster)	计算栅格中像元值的平方值
平方根	outSQRT = SquareRoot(inRaster)	计算栅格中像元值的平方根
幂	outPower = Power(inRaster1, inRaster2)	对栅格中的像元值进行乘方运算，结果为在另一个栅格中找到的值
取反	outNegate = Negate(inRaster)	逐个像元地更改输入栅格的像元值符号（乘以 -1）
绝对值	outAbs = Abs(inRaster)	计算栅格中像元的绝对值
Exp	outExp = Exp(inRaster)	计算栅格中各像元以 e 为底的指数
Exp2	outExp2 = Exp2(inRaster)	计算栅格中各像元以 2 为底的指数
Exp10	outExp10 = Exp10(inRaster)	计算栅格中各像元以 10 为底的指数
Ln	outLn = Ln(inRaster)	计算栅格中各像元的自然对数(以 e 为底)
Log2	outLog2 = Log2(inRaster)	计算栅格中各像元以 2 为底的对数
Log 10	outLog10 = Log10(inRaster)	计算栅格中各像元以 10 为底的对数
上舍入	outRoundURaster = RoundUp(inRaster)	返回栅格中每个像元的最近的较大整数，但表示为浮点型
下舍入	outRoundDRaster = RoundDown(inRaster)	返回栅格中每个像元的最近的较小整数，但表示为浮点型
求模	outMod = Mod(inRaster1, inRaster2)	逐个像元地求出第一个栅格数据除以第二个栅格数据的余数（模）
转为整型	outInt = Int(inRaster)	通过截断将栅格的每个像元值转换为整型
转为浮点型	outFloat = Float(inRaster)	将每个栅格像元的值转换为浮点型表达形式

实例及运行结果

由于Arcpy对一些简单的栅格加、减、乘、除等运算进行了符号重载，因此可以直接用+、-、*、/等运算符来进行栅格和常量、栅格与栅格之间的数学运算，PyCharm中的计算NPP的具体实例如下所示：

// ArcPy计算NPP代码-Python by jjg
#-*- coding: UTF-8 -*-
import arcpy
from arcpy import env
from arcpy.sa import *

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")
#读取NPP模型的输入数据
NDVI_Raster = arcpy.Raster("H:\\CASA_Model_200101\\200101输入数据\\NDVI\\MOD13Q1.A200101.hdf.250m_16_days_NDVI_projected1_mask(0-1).tif")#读取tif文件
SOL_Raster = arcpy.Raster("H:\\CASA_Model_200101\\200101输入数据\\sun_radiation\\Y2001M1_sun_radi.img")#读取tif文件
ET_Raster = arcpy.Raster("H:\\CASA_Model_200101\\200101输入数据\\ET\\MOD16A2.200101ET_500m_MVC_Projected_mask.tif")#读取tif文件
PET_Raster = arcpy.Raster("H:\\CASA_Model_200101\\200101输入数据\\PET\\MOD16A2.A200101.hdf.PET_500m_MVC_projected_mask.tif")#读取tif文件
T_Raster = arcpy.Raster("H:\\CASA_Model_200101\\200101输入数据\\TAVG_CORR\\Y2001M1_TAVG_CORR.img")#读取tif文件
SRmin = 0.5
SRmax = 1.5
FPARmin = 0.001
FPARmax = 0.95
T_opt = 26
e_max = 0.389

SR = (1 + NDVI_Raster) / (1 - NDVI_Raster)
FPAR = FPARmin + (FPARmax - FPARmin) * ((SR - SRmin) / (SRmax - SRmin))
APAR = SOL_Raster * FPAR * 0.5
def func(a):
    if a > -10:
        return 0.8 + 0.02 * a - 0.0005 * a * a
    else:
        return 0

T1 = func(T_opt)
T2 = (1.184/(1+Exp(0.2*(T_opt-10-T_Raster))))*(1/(1+Exp(0.3*(T_Raster -10-T_opt))))
W = 0.5 +0.5* ET_Raster / PET_Raster
e = T1 *T2 * W *e_max
NPP = APAR * e
NPP.save("H:\\CASA_Model_200101\\NPP\\NPP.tif")