ArcPy学习心得（4）(四)大尺度空间范围可燃物载量计算方法

可燃物载量表示在单位面积范围内可燃物的负荷载量总和，一般来说，可燃物载量越高，火灾隐患越大，而不同地物的可燃物载量是不同的，本文结合了可燃物载量的主流计算模型，利用NDVI、坡度、土地利用分类数据相结合的方式来对全国各地的可燃物载量进行计算(部分算法只公开思路，不公开细节)，而这篇文章，主要阐述了如何使用Arcpy在大尺度空间范围下对数据进行更高效的处理及计算。

在进行代码设计之前，让我们先分析一下此次计算所需要的数据都有哪些吧！

根据计算可燃物载量的步骤，我们需要先对NDVI数据、坡度数据、土地利用数据进行读取，只有在获取到数据之后，才能够进行下一步的计算，以下为数据读取部分的代码，其中使用了arcpy的Raster方法对三类数据进行读取。

NDVITif = arcpy.Raster(NDVIPath)
SlopeTif = arcpy.Raster(SlopePath)
ReclassTif = arcpy.Raster(ReclassPath)

 

利用arcpy的Raster方法将数据读取为Raster后，需要使用arcpy的RasterToNumPyArray方法，将数据从Raster转换为Array格式的数据，其中Array为二维数组,以下为读取数据数组的代码。

NoData = Reclass.noDataValue#读取Nodata
#读取Reclass
dataReclass=arcpy.RasterToNumPyArray(Reclass,nodata_to_value=NoData)
rows,cols = dataReclass.shape#获取行数和列数
#读取NDVI
NoData = NDVI.noDataValue#读取Nodata
dataNDVI=arcpy.RasterToNumPyArray(NDVI, nodata_to_value=NoData)
#读取Slope
NoData = Slope.noDataValue#读取Nodata
dataSlope=arcpy.RasterToNumPyArray(Slope, nodata_to_value=NoData)

在获取到数据的数组后，即可对数据进行计算，然而，在大尺度空间范围下，数据计算过程中所占用的运行内存是非常多的，在数据量较大的情况下，会发生内存溢出的情况，针对于此类情况，Numpy的Array数组有一个特性，即可对Array数组进行分块处理，也可理解为切片处理，而在对Array数组进行分块处理时，可选择规则型分块，即根据数组的行数和列数，将数组分成大小相同的矩形，对其进行处理，也可利用条件表达式，对Array数组进行过滤，对过滤后的数组进行处理。

而在此次计算过程中，因为不同土地利用类型下的可燃物载量也不同，因此我们选择将土地利用类型作为筛选条件，来对数组进行过滤，并通过NDVI以及坡度分块计算各地类的可燃物载量。

以下为实现过程的代码(部分算法代码不予展示)。

S = np.zeros([rows,cols],dtype=float)
S[dataReclass==10]=FuelloadQ(dataNDVI[dataReclass==10],dataSlope[dataReclass==10],10)
S[dataReclass==20]=FuelloadQ(dataNDVI[dataReclass==20],dataSlope[dataReclass==20],20)
S[dataReclass==30]=FuelloadQ(dataNDVI[dataReclass==30],dataSlope[dataReclass==30],30)

在上段代码中，数组S为最终的计算结果，在计算过程中，首先利用了条件表达式“dataReclass==n,n=(10,20,30)”，对数组S、dataNDVI、dataSlope进行了分块，分别提取出了其中以农田、林地、草地为下垫面的数组，随后将分块后的数据传入至FuelloadQ方法中对其进行可燃物载量的计算，即可在最大程度上节省脚本运行时所花费的时间成本与空间成本。

最后，将计算好的数组S通过arcpy的NumPyArrayToRaster方法转换为Raster，随后对这个Raster进行保存，即可完成所有计算步骤。以下为保存数组S的代码。

NoData = 0#读取Nodata
ExtentXmin = tif.extent.XMin#取x坐标最小值
ExtentYmin = tif.extent.YMin#取y坐标最小值
lowerLeft = arcpy.Point(ExtentXmin,ExtentYmin)#取得数据起始点坐标(左下角坐标)
Widthcellsize=tif.meanCellWidth#像元宽度
Heightcellsize=tif.meanCellHeight#像元高度
new_raster=arcpy.NumPyArrayToRaster(S,lowerLeft,Widthcellsize,Heightcellsize,value_to_nodata=NoData)
new_raster.save(SavePath)
arcpy.BuildPyramids_management(SavePath)

其中，NumPyArrayToRaster方法所需要的参数有：数组S、数据的起始点坐标(左下角坐标)、像元宽度、像元高度、NoData等参数。

Arcpy中的BuildPyramids_management方法的作用是为保存后的Raster建立金字塔。

以下内容为此次代码的部分内容(不包含可燃物载量算法)。

import numpy as np
import arcpy

def Fuelload(NDVI,Reclass,Slope):
    #获取数据矩阵
    NoData = Reclass.noDataValue#读取Nodata
    dataReclass = arcpy.RasterToNumPyArray(Reclass, nodata_to_value=NoData)
    rows,cols = dataReclass.shape#获取行数和列数
    #读取NDVI
    NoData = NDVI.noDataValue#读取Nodata
    dataNDVI = arcpy.RasterToNumPyArray(NDVI, nodata_to_value=NoData)
    
    #读取Slope
    NoData = Slope.noDataValue#读取Nodata
    dataSlope = arcpy.RasterToNumPyArray(Slope, nodata_to_value=NoData)

    S = np.zeros([rows,cols],dtype=float)
    
    S[dataReclass == 10] = FuelloadQ(dataNDVI[dataReclass == 10],dataSlope[dataReclass == 10],10)#农田
    S[dataReclass == 20] = FuelloadQ(dataNDVI[dataReclass == 20],dataSlope[dataReclass == 20],20)#林地
    S[dataReclass == 30] = FuelloadQ(dataNDVI[dataReclass == 30],dataSlope[dataReclass == 30],30)#草地
    del dataNDVI,dataReclass,dataSlope
    gc.collect()
    return S

def Run(NDVIPath,ReclassPath,SlopePath,SavePath):
    NDVITif = arcpy.Raster(NDVIPath)
    ReclassTif = arcpy.Raster(ReclassPath)
    SlopeTif = arcpy.Raster(SlopePath)

    S = Fuelload(NDVITif,ReclassTif,SlopeTif)

    #保存
    NoData = 0#读取Nodata
    ExtentXmin = tif.extent.XMin#取x坐标最小值
    ExtentYmin = tif.extent.YMin#取y坐标最小值
    lowerLeft = arcpy.Point(ExtentXmin,ExtentYmin)#取得数据起始点坐标(左下角坐标)
    Widthcellsize=tif.meanCellWidth#像元宽度
    Heightcellsize=tif.meanCellHeight#像元高度

    new_raster = arcpy.NumPyArrayToRaster(S,lowerLeft,Widthcellsize,Heightcellsize,value_to_nodata=NoData)
    new_raster.save(SavePath)
    arcpy.BuildPyramids_management(SavePath)