ArcPy学习心得系列（1）自定义重采样方法

题目：

将原图像的4个像元合并为一个像元，即取2*2窗口的四个像元，取其平均值作为输出图像的像素，输出图像坐标应与原图像相同，输出像元大小应为原图像两倍。

实现思路：

在实现以上需求之前，让我们先来构思一下做出这个结果需要哪些数据吧！

1.根据需求，我们首先要获取的是栅格数据的数组，其次才能对栅格数据进行处理，因此，我们要对栅格数据的数组信息进行提取。

2.在获取栅格数据的数组之前，数组中难免会混进去一些我们不需要的数据，一般而言我们称其为NoDATA，因此在获取整个栅格数组之前，我们还要利用NoDATA对获取到的数组进行过滤。

3.取每四个像元的平均值的时候，因为我们对整个数组并不了解，如果输入进来的数据的行数和列数大于4的时候，我们就需要使用迭代的方式对整个数组进行遍历，而这种迭代的方式需要先获取到栅格数据的行数和列数才能进行(栅格数据的行数和列数也代表了数组的行数和列数)，因此我们需要在迭代开始之前取得数据的行数和列数。

4.输出图像坐标系应与原图像相同，在这一需求上，我们除了要定义数据的输出坐标系外，还要对我们输出数据的相对位置做出定义，也可理解为控制点，在arcpy中，新生成的数组在转化为栅格数据时，需指定初始坐标位置(左下角位置)，因此在在新生成的数据转化为栅格数据之前，我们需要先取得原图像的左下角位置坐标。

5.在将四个像元合为一个像元后，原图像的像元大小也应该做出相应改变，理应为原先图像像元大小的二倍，因此我们需要取得原始图像的像元大小。

梳理一下：

以下为我们需要在原图像中提取到的信息

1.数组信息

2.NoDATA

3.行数和列数

4.左下角坐标位置

5.栅格像元大小

以下为提取原始图像信息的方法

首先，利用ArcPy中的Raster方法读取栅格

tif = arcpy.Raster(tiffFileName)#读取栅格，从读取到的tif中取NoData、行数、列数、像元大小、数据范围这些属性，提取方法为
NoData = tif.noDataValue#读取Nodata
rows,cols = data.shape#获取行数和列数
cellsize=tif.meanCellWidth#像元大小
ExtentXmin = tif.extent.XMin#取x坐标最小值
ExtentYmin = tif.extent.YMin#取y坐标最小值
lowerLeft = arcpy.Point(ExtentXmin,ExtentYmin )#取得数据起始点坐标(左下角坐标)

 利用RasterToNumPyArray方法将栅格数据转换为array数组

 

data = arcpy.RasterToNumPyArray(tif, nodata_to_value=NoData)#数据转array

 

接下来，我们需要对新生成数据的组织结构进行分析：

1）新数据每个单独像元的值由原图像四个像元的均值组成；

2）像元大小为原先图像像元大小的二倍；

3）行数和列数为原先图像行列数的二分之一+n(当原图像的行(列)为偶数时n=0，否则n=1)

我们在进行计算之前，先创建一个行列数为原图像1/2大小的数组arrdata，其中%表示取余，dtype设置为原先图像的dtype

arrdata = np.zeros([int(rows/2)+rows%2,int(cols/2)+cols%2],dtype=data.dtype)#创建一个新数组，行列数等于原来的一半 

 随后利用迭代法结合所给出的需求对新的栅格数组进行赋值，每次前进两格连一下

for x in range(rows%2,rows,2):#进行迭代，从rows%2开始，到rows结束，每次前进2    for y in range(cols%2,cols,2):#进行迭代，从cols%2开始，到cols结束，每次前进2        newdata = data[x:x+2,y:y+2]#取得2*2像元大小        if NoData in newdata:#如果取到的四个像元中有NoDATA则跳过本次循环            continue        arrdata[x/2-x%2,y/2-y%2] = np.average(newdata)#给新数组赋值，新值为2*2像元的均值

 

 

迭代完成后，根据原图像像元大小计算新像元大小newcellsize = cellsize*2，利用NumPyArrayToRaster方法将数组转化为栅格数据

 

 

#将numpy转栅格，输入数据、左下角坐标、像元大小、nodatanew_raster = arcpy.NumPyArrayToRaster(arrdata,lowerLeft,newcellsize,newcellsize, value_to_nodata=NoData)arcpy.env.overwriteOutput = True#设置保存时可以覆盖已有数据new_raster.save(outputPath)#保存数据arcpy.BuildPyramids_management(outputPath)#建立金字塔

 至此，所有需求解决完成，以下为此次分享代码的全部内容

import arcpy
import numpy as np

def ResampleF(tiffFileName,outputPath):
    arcpy.env.outputCoordinateSystem = tiffFileName  # 输出坐标系与输入相同
    tif = arcpy.Raster(tiffFileName)#读取栅格
    NoData = tif.noDataValue#读取Nodata
    rows,cols = data.shape#获取行数和列数
    cellsize=tif.meanCellWidth#像元大小
    ExtentXmin = tif.extent.XMin#取x坐标最小值
    ExtentYmin = tif.extent.YMin#取y坐标最小值
    lowerLeft = arcpy.Point(ExtentXmin,ExtentYmin )#取得数据起始点范围
    
    data = arcpy.RasterToNumPyArray(tif, nodata_to_value=NoData)#数据转array
    
    arrdata = np.zeros([int(rows/2)+rows%2,int(cols/2)+cols%2],dtype=data.dtype)#创建一个新数组，行列数等于原来的一半
    
    for x in range(rows%2,rows,2):#进行循环，从rows%2开始，到rows结束，每次前进2
        for y in range(cols%2,cols,2):#进行循环，从cols%2开始，到cols结束，每次前进2
            newdata = data[x:x+2,y:y+2]#取得2*2像元大小
            if NoData in newdata:#如果有Nodata则跳过循环
                continue
            arrdata[x/2-x%2,y/2-y%2] = np.average(newdata)#给新数组赋值，赋值为2*2像元均值

            
    newcellsize = cellsize*2#计算新数据的像元大小
    
    #将numpy转栅格，输入数据、范围、像元大小、nodata
    new_raster = arcpy.NumPyArrayToRaster(arrdata,lowerLeft,cellsize*2,cellsize*2, value_to_nodata=NoData)

    arcpy.env.overwriteOutput = True#可以覆盖数据
    new_raster.save(outputPath)#保存数据
    arcpy.BuildPyramids_management(outputPath)#建立金字塔