#批量处理降雨数据,涉及批量的矢量文件处理,批量文件空间插值,批量分区统计,批量写入TXT
代码主要完成将降雨的站点数据(包括站点雨量信息,坐标)转化成矢量数据,然后就行空间插值,将空间插值的结果提取到面的几何中心,然后再将,几何中心的坐标,雨量站信息,插值后的值转化为时间序列,按格式写入txt文件实现的整体代码如下(!!!注意文件放的位置!!!):
# coding=utf-8
"""
将降雨的站点数据(包括站点雨量信息,坐标)转化成矢量数据,然后就行空间插值,将空间插值的结果提取到面的几何中心,然后再将,几何中心的坐标,雨量站信息,
插值后的值转化为时间序列,按格式写入txt文件实现步骤如下:
1:将excel的数据转化为矢量文件
2:将矢量文件进行空间插值
3:找出面的几何中心
4:空间插值的结果做分区统计
5:将面的几何中心(坐标信息)写入txt(雨量站)
6:将分区统计的结果转化成时间序列写入TXT
7:将写好的时间序列进行按雨量站名称排序
8:将站点信息加入子汇水区
"""
import sys
import xlrd
import arcpy
from arcpy import env
from arcpy.sa import *
import os
reload(sys)
sys.setdefaultencoding('utf-8')
# 1:将excel的数据转化为矢量文件
def readExcelT_Shp(indir, outPath):
path = outPath
if os.path.exists(path):
pass
else:
os.mkdir(path)
xl1 = xlrd.open_workbook(indir)
xl1 = xlrd.open_workbook(indir)
sheet = xl1.sheet_by_name('Sheet1')
lines = sheet.nrows
dirOrder = 0
spatial_reference = arcpy.SpatialReference("WGS 1984")
for i in range(1, lines):
if sheet.row_values(i)[1].decode('utf-8') == '双流':
dirOrder = dirOrder + 1
print(i)
outShp = str(dirOrder) + 'raining.shp'
print(outShp)
arcpy.CreateFeatureclass_management(outPath, outShp, "POINT", spatial_reference=spatial_reference)
print('************1***********')
arcpy.AddField_management(outPath + '\\' + outShp, "raining", "FLOAT")
print("生成文件:" + outShp)
print(sheet.row_values(i)[4].decode('utf-8'))
rows = arcpy.InsertCursor(outPath + '\\' + outShp)
point = arcpy.Point()
longitude = sheet.row_values(i)[2]
latitude = sheet.row_values(i)[3]
rainValue = sheet.row_values(i)[5]
row = rows.newRow()
point.X = longitude
print(point.X)
point.Y = latitude
print(point.Y)
pointGeometry = arcpy.PointGeometry(point)
row.shape = pointGeometry
print(rainValue)
print('***********%d************' % i)
row.raining = rainValue
rows.insertRow(row)
del row
del rows
# 2:将矢量文件进行空间插值
def interpolation(workspace, outfile):
env.workspace = workspace
path = outfile
if os.path.exists(path):
pass
else:
os.mkdir(path)
for i in range(1, 63):
print("**********%d**********" % i)
inPointFeatures = str(i) + 'raining.shp'
zField = 'raining'
cellSize = 0.00081411112
arcpy.CheckOutExtension("Spatial")
# Execute Spline
outSplineBarriers = NaturalNeighbor(inPointFeatures, zField, cellSize)
outSplineBarriers.save(outfile + "\\" + str(i) + 'result_NaturalNeighbor.tif')
# 3:找出面的几何中心
def resPlygonCenter(workspace, inpolygon, SpatialReference, output_location):
path = output_location
if os.path.exists(path):
pass
else:
os.mkdir(path)
env.workspace = workspace
featuresList = []
in_polygon = inpolygon
polygon_fields = ['SHAPE@TRUECENTROID', 'name']
print("********************************************")
with arcpy.da.SearchCursor(in_polygon, polygon_fields) as cursor:
for row in cursor:
# 将数据添加到集合中去
featuresList.append([row[0], row[1]])
# 结束循环
fc_name = 'Point_Layer'
outputpath = output_location
sr = arcpy.SpatialReference(SpatialReference)
arcpy.CreateFeatureclass_management(outputpath, fc_name, 'POINT', spatial_reference=sr)
point_layer = os.path.join(output_location, fc_name)
arcpy.AddField_management(point_layer + ".shp", 'NAME', 'TEXT')
point_fields = ['SHAPE@XY', 'NAME']
with arcpy.da.InsertCursor(point_layer + ".shp", point_fields) as cursor:
for record in featuresList:
cursor.insertRow(record)
# 4:空间插值的结果做分区统计
def static_byZone(workspace, zone, field, tiff, outpath):
path = outpath
if os.path.exists(path):
pass
else:
os.mkdir(path)
os.mkdir(txtFile)
txtfile = outpath + "\\txt\\" + str(i) + "staticData.txt"
arcpy.CheckOutExtension("spatial")
env.overwriteOutput = 1
workingDir = workspace
env.workspace = workingDir
in_zone_data = zone
zone_field = field
hour = 9
day = 28
for i in range(1, 63):
hourString = hour
if hour < 10:
hourString = '0' + str(hour)
time = "07/" + str(day) + "/2020" + " " + str(hourString) + ":00"
in_value_raster = workspace + "\\" + str(i) + tiff
out_table = outpath + "\\" + str(
i) + "static.dbf" #
outZSaT = ZonalStatisticsAsTable(in_zone_data, zone_field, in_value_raster, out_table, "DATA", "MEAN")
txtFile = outpath + "\\txt\\"
with open(txtfile, "w") as f:
fields = ['name', 'MEAN'] # 指定要读入的字段
with arcpy.da.SearchCursor(out_table, fields) as rows:
for row in rows:
mean = round(row[1], 3)
print('{0},{1}'.format("R" + row[0], mean))
f.write('R{0} {1} {2}'.format(row[0], time, mean))
f.write('\n') # 写完一行后换行
hour = hour + 1
if hour == 24:
hour = 0
day = day + 1
f.close()
# 5:将面的几何中心(坐标信息)写入txt(雨量站)
def write_CenterPoint(workspace, outfile):
output_location = workspace
env.workspace = output_location
fc_name = 'Point_Layer'
point_layer = os.path.join(output_location, fc_name)
file = open(outfile, 'r+')
flag = 1
fc = point_layer + '.shp'
with arcpy.da.SearchCursor(fc, ["SHAPE@XY", "name"]) as cursor:
for row in cursor:
print("这是第" + str(flag) + "个点要素")
flag += 1
# 打印每个点要素的x,y坐标
x, y = row[0]
name = 'R' + row[1]
print("{}, {} , {}".format(name, x, y))
file.write(("{} {} {}".format(name, x, y)))
file.write("\n")
file.close()
print("写入完毕")
# 6:将分区统计的结果转化成时间序列写入TXT
def read_Static_Wserise(inpath, oupfile):
dirpath = inpath
output_location = oupfile
for i in range(1, 63):
readfile = open(dirpath + '\\' + str(i) + 'staticData.txt', 'r')
writeFiles = open(output_location, 'a')
for lines in readfile:
writeFiles.write(lines)
readfile.close()
writeFiles.close()
# 7:将写好的时间序列进行按雨量站名称排序
def sortSerise(file):
files = open(file, 'r')
list = []
for lines in files:
list.append(lines)
files.close()
files = open(file, 'w')
files.write("")
files.close()
files = open(file, 'a')
list.sort()
i = 0
for item in list:
files.write(item)
i = i + 1
if i == 62:
files.write(';\n')
i = 0
files.close()
# 8:将站点信息加入子汇水区
def writeRage_subcachment(inShpPath, outfile):
env.workspace = inShpPath
fc_name = 'subcachment'
output_location = inShpPath
point_layer = os.path.join(output_location, fc_name)
flag = 1
file = open(outfile, 'r+')
fc = point_layer + '.shp'
with arcpy.da.SearchCursor(fc, ["name", "nodeNumber", "Area", "F_Imp", "width", "slope"]) as cursor:
for row in cursor:
print("这是第" + str(flag) + "个点要素")
flag += 1
# 打印每个点要素的x,y坐标
Name = row[0]
RainGage = 'R' + row[0]
Outlet = row[1]
Area = row[2]
imp = row[3]
width = row[4]
Slope = row[5]
cur = 0
print("{0}, {1} , {2} , {3}".format(Name, RainGage, Outlet, Area))
print("{}, {} , {} , {}".format(imp, width, Slope, cur))
file.write(("{} {} {} {} ".format(Name, RainGage, Outlet, Area)))
file.write(("{} {} {} 0".format(imp, width, Slope)))
file.write("\n")
file.close()
print("写入完毕")
def main():
readExcelT_Shp(r'E:\\Daily Projiect\\SL_interpolationZonalStatic\\data\\rain1.xls',
r'E:\\Daily Projiect\\SL_interpolationZonalStatic\\pointShpfile')
interpolation(r'E:\\Daily Projiect\\SL_interpolationZonalStatic\\pointShpfile',
r'E:\\Daily Projiect\\SL_interpolationZonalStatic\\interpolationResult')
resPlygonCenter(r'E:\\Daily Projiect\\SL_interpolationZonalStatic\\data\\polygon',
r'E:\\Daily Projiect\\SL_interpolationZonalStatic\\data\\polygon\\subcachment.shp',
r'E:\\Daily Projiect\\SL_interpolationZonalStatic\\data\\polygon\\subcachment.prj',
r'E:\\Daily Projiect\\SL_interpolationZonalStatic\\result_PlygonCenter')
static_byZone(r'E:\\Daily Projiect\\SL_interpolationZonalStatic\\interpolationResult',
r"E:\\Daily Projiect\\SL_interpolationZonalStatic\\data\polygon\\subcachment.shp",
"name",
"result_NaturalNeighbor.tif",
r'E:\\Daily Projiect\\SL_interpolationZonalStatic\\staticResult')
write_CenterPoint(r'E:\\Daily Projiect\\SL_interpolationZonalStatic\\result_PlygonCenter',
r'E:\\Daily Projiect\\SL_interpolationZonalStatic\\ragePosition.txt')
read_Static_Wserise(r'E:\\Daily Projiect\\SL_interpolationZonalStatic\\staticResult\TXT',
r'E:\\Daily Projiect\\SL_interpolationZonalStatic\\timeSeries.txt')
sortSerise(r'E:\\Daily Projiect\\SL_interpolationZonalStatic\\timeSeries.txt')
writeRage_subcachment(r'E:\\Daily Projiect\\SL_interpolationZonalStatic\\data\\polygon',
r'E:\\Daily Projiect\\SL_interpolationZonalStatic\\subcachment.txt')
if __name__ == '__main__':
main()
步骤如下(每个步骤都是一个小功能):
- 1.将excel的数据转化为矢量文件
站点信息包括经纬度和降雨量就行
def readExcelT_Shp(indir, outPath):
path = outPath
if os.path.exists(path):
pass
else:
os.mkdir(path)
xl1 = xlrd.open_workbook(indir)
xl1 = xlrd.open_workbook(indir)
sheet = xl1.sheet_by_name('Sheet1')
lines = sheet.nrows
dirOrder = 0
spatial_reference = arcpy.SpatialReference("WGS 1984")
for i in range(1, lines):
if sheet.row_values(i)[1].decode('utf-8') == '双流':
dirOrder = dirOrder + 1
print(i)
outShp = str(dirOrder) + 'raining.shp'
print(outShp)
arcpy.CreateFeatureclass_management(outPath, outShp, "POINT", spatial_reference=spatial_reference)
print('************1***********')
arcpy.AddField_management(outPath + '\\' + outShp, "raining", "FLOAT")
print("生成文件:" + outShp)
print(sheet.row_values(i)[4].decode('utf-8'))
rows = arcpy.InsertCursor(outPath + '\\' + outShp)
point = arcpy.Point()
longitude = sheet.row_values(i)[2]
latitude = sheet.row_values(i)[3]
rainValue = sheet.row_values(i)[5]
row = rows.newRow()
point.X = longitude
print(point.X)
point.Y = latitude
print(point.Y)
pointGeometry = arcpy.PointGeometry(point)
row.shape = pointGeometry
print(rainValue)
print('***********%d************' % i)
row.raining = rainValue
rows.insertRow(row)
del row
del rows
- 2.将矢量文件进行空间插值
这一步直接使用上一步得到矢量文件就行
def interpolation(workspace, outfile):
env.workspace = workspace
path = outfile
if os.path.exists(path):
pass
else:
os.mkdir(path)
for i in range(1, 63):
print("**********%d**********" % i)
inPointFeatures = str(i) + 'raining.shp'
zField = 'raining'
cellSize = 0.00081411112
arcpy.CheckOutExtension("Spatial")
outSplineBarriers = NaturalNeighbor(inPointFeatures, zField, cellSize) #插值的方法在这里改,注意看arcpy的帮助文档,不同插值所需的变量不同
outSplineBarriers.save(outfile + "\\" + str(i) + 'result_NaturalNeighbor.tif')
- 3.找出面的几何中心
data文件里面的面矢量文件的几何中心
def resPlygonCenter(workspace, inpolygon, SpatialReference, output_location):
path = output_location
if os.path.exists(path):
pass
else:
os.mkdir(path)
env.workspace = workspace
featuresList = []
in_polygon = inpolygon
polygon_fields = ['SHAPE@TRUECENTROID', 'name']
print("********************************************")
with arcpy.da.SearchCursor(in_polygon, polygon_fields) as cursor:
for row in cursor:
# 将数据添加到集合中去
featuresList.append([row[0], row[1]])
# 结束循环
fc_name = 'Point_Layer'
outputpath = output_location
sr = arcpy.SpatialReference(SpatialReference)
arcpy.CreateFeatureclass_management(outputpath, fc_name, 'POINT', spatial_reference=sr)
point_layer = os.path.join(output_location, fc_name)
arcpy.AddField_management(point_layer + ".shp", 'NAME', 'TEXT')
point_fields = ['SHAPE@XY', 'NAME']
with arcpy.da.InsertCursor(point_layer + ".shp", point_fields) as cursor:
for record in featuresList:
cursor.insertRow(record)
- 4.空间插值的结果做分区统计
将第二步空间插值的结果分区统计给我们的面矢量,并写入TXT
def static_byZone(workspace, zone, field, tiff, outpath):
path = outpath
if os.path.exists(path):
pass
else:
os.mkdir(path)
os.mkdir(txtFile)
txtfile = outpath + "\\txt\\" + str(i) + "staticData.txt"
arcpy.CheckOutExtension("spatial")
env.overwriteOutput = 1
workingDir = workspace
env.workspace = workingDir
in_zone_data = zone
zone_field = field
hour = 9
day = 28
for i in range(1, 63):
hourString = hour
if hour < 10:
hourString = '0' + str(hour)
time = "07/" + str(day) + "/2020" + " " + str(hourString) + ":00"
in_value_raster = workspace + "\\" + str(i) + tiff
out_table = outpath + "\\" + str(
i) + "static.dbf" #
outZSaT = ZonalStatisticsAsTable(in_zone_data, zone_field, in_value_raster, out_table, "DATA", "MEAN")
txtfile = outpath + "\\txt\\" + str(i) + "staticData.txt"
with open(txtfile, "w") as f:
fields = ['name', 'MEAN'] # 指定要读入的字段
with arcpy.da.SearchCursor(out_table, fields) as rows:
for row in rows:
mean = round(row[1], 3)
print('{0},{1}'.format("R" + row[0], mean))
f.write('R{0} {1} {2}'.format(row[0], time, mean))
f.write('\n')
hour = hour + 1
if hour == 24:
hour = 0
day = day + 1
f.close()
- 5.将面的几何中心(坐标信息)写入txt(雨量站)
def write_CenterPoint(workspace, outfile):
output_location = workspace
env.workspace = output_location
fc_name = 'Point_Layer'
point_layer = os.path.join(output_location, fc_name)
file = open(outfile, 'r+')
flag = 1
fc = point_layer + '.shp'
with arcpy.da.SearchCursor(fc, ["SHAPE@XY", "name"]) as cursor:
for row in cursor:
print("这是第" + str(flag) + "个点要素")
flag += 1
# 打印每个点要素的x,y坐标
x, y = row[0]
name = 'R' + row[1]
print("{}, {} , {}".format(name, x, y))
file.write(("{} {} {}".format(name, x, y)))
file.write("\n")
file.close()
print("写入完毕")
- 6.将分区统计的结果转化成时间序列写入TXT
def read_Static_Wserise(inpath, oupfile):
dirpath = inpath
output_location = oupfile
for i in range(1, 63):
readfile = open(dirpath + '\\' + str(i) + 'staticData.txt', 'r')
writeFiles = open(output_location, 'a')
for lines in readfile:
writeFiles.write(lines)
readfile.close()
writeFiles.close()
- 7.将写好的时间序列进行按雨量站名称排序
def sortSerise(file):
files = open(file, 'r')
list = []
for lines in files:
list.append(lines)
files.close()
files = open(file, 'w')
files.write("")
files.close()
files = open(file, 'a')
list.sort()
i = 0
for item in list:
files.write(item)
i = i + 1
if i == 62:
files.write(';\n')
i = 0
files.close()
-8. 将站点信息加入子汇水区
def writeRage_subcachment(inShpPath, outfile):
env.workspace = inShpPath
fc_name = 'subcachment'
output_location = inShpPath
point_layer = os.path.join(output_location, fc_name)
flag = 1
file = open(outfile, 'r+')
fc = point_layer + '.shp'
with arcpy.da.SearchCursor(fc, ["name", "nodeNumber", "Area", "F_Imp", "width", "slope"]) as cursor:
for row in cursor:
print("这是第" + str(flag) + "个点要素")
flag += 1
# 打印每个点要素的x,y坐标
Name = row[0]
RainGage = 'R' + row[0]
Outlet = row[1]
Area = row[2]
imp = row[3]
width = row[4]
Slope = row[5]
cur = 0
print("{0}, {1} , {2} , {3}".format(Name, RainGage, Outlet, Area))
print("{}, {} , {} , {}".format(imp, width, Slope, cur))
file.write(("{} {} {} {} ".format(Name, RainGage, Outlet, Area)))
file.write(("{} {} {} 0".format(imp, width, Slope)))
file.write("\n")
file.close()
print("写入完毕")