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) #在调用盆域计算函数