Arcpy | 裁剪栅格和并行计算

import arcpy
import numpy

def RastertoNumpy(path,nodata_to_value): #（txt的路径，nodata值）
    txt = open(path,"rb")            #‘rb’：表示以二进制方式读取文件。该文件必须已存在。
    row = ""
    col = ""
    for i in range(0,6):            #头文件是6行
        if i == 0:
            txt.seek(14)            #光标向前移动14个位置，获取ncols,: f.seek(offset, whence)whence: 参考位置，一般参数为0,1,2 0 ：将开头作为参考位置 1 ：将当前作为参考位置 2 ：将末尾作为参考位置
            while(1):
                character = txt.read(1).decode(encoding = "utf-8")
                if character == "\r":
                    break
                col = col + character
        if i == 1:
            txt.seek(14,1)
            while(1):
                character = txt.read(1).decode(encoding = "utf-8")
                if character == "\r":
                    break
                row = row + character
        txt.readline()
    row = int(row)
    col = int(col)
    result = numpy.zeros((row,col),dtype="uint16")
    flowdirectionvalue = ""
    i = 0
    j = 0
    while(1):
        character = txt.read(1).decode(encoding = "utf-8")
        if character not in {" ", "\r", "\n", ""}:
            flowdirectionvalue = flowdirectionvalue + character
        elif character != "":
            if flowdirectionvalue == "-9999":
                flowdirectionvalue = int(nodata_to_value)
            if flowdirectionvalue != "":
                flowdirectionvalue = int(flowdirectionvalue)
                result[i][j] = flowdirectionvalue
                j += 1
            if j == col:                              #换行
                j = 0
                i = i + 1
            flowdirectionvalue = ""
        elif character == "":
            break
    txt.close()
    return result                                     #返回一个np数组


def clip(raster,txt,row,col): 
    np1 = numpy.zeros((1000,1000),dtype = "uint16")
    numpyraster = RastertoNumpy(txt, 247)
    temp_cellsizex = arcpy.GetRasterProperties_management(raster, "CELLSIZEX")
    cellsizex = float(temp_cellsizex.getOutput(0))
    cellsizey = cellsizex
    inraster = arcpy.Raster(raster, True)
    dy = (inraster.extent.YMax - inraster.extent.YMin) / row      #纬度
    dx = (inraster.extent.XMax - inraster.extent.XMin) / col      #经度
    number = 0
    for i in range(0,row):
        for j in range(0,col):
            number += 1
            start_x = j * 1000
            end_x = start_x + 1000
            start_y = i * 1000
            end_y = start_y + 1000
            for m in range(start_y,end_y):
                for n in range(start_x,end_x):
                    np1[m-start_y][n-start_x] = numpyraster[m][n]
            lowerleft_point = arcpy.Point(inraster.extent.XMin+j*dx, inraster.extent.YMax-(i+1)*dy)
            result = arcpy.NumPyArrayToRaster(np1, lowerleft_point, cellsizex, cellsizey, 247)
            path = r"D:\河网\dir2\No" + str(number-1) + "(dir).tif"
            result.save(path)

场景：将大型栅格裁剪成1000*1000的瓦片

Raster：dem.tif

首先将栅格转成txt，得到dem.txt，RastertoNumpy(path,nodata_to_value)可以函数读取txt中的栅格。

clip函数进行裁剪

2.对流向栅格并行计算盆域

import arcpy
import os
from concurrent.futures import  ProcessPoolExecutor

def delineate_watershed(PARAMS):
    inputcontent = PARAMS[0]
    outputcontent = PARAMS[1]
    outbasin = arcpy.sa.Basin(inputcontent)
    outbasin.save(outputcontent)

def multi_process(PARAMS):
    with ProcessPoolExecutor(max_workers = 8) as pool:
        pool.map(delineate_watershed,PARAMS)

path = r"D:\河网\dir2"
arcpy.env.workspace = path
result_path = r"D:\河网\result"
PARAMS = []
for tifFiles in arcpy.ListFiles("*.tif"):                #先将输入数据与输出数据的name放入数组中
    info = os.path.basename(tifFiles).split('(')[0]      #返回path最后的文件名
    inputcontent = path + "/" + tifFiles
    outputcontent = result_path + "/" + info + "(basin).tif"
    PARAMS.append([inputcontent,outputcontent])
multi_process(PARAMS)                                    #在调用盆域计算函数