【无标题】

#批量处理降雨数据,涉及批量的矢量文件处理,批量文件空间插值,批量分区统计,批量写入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("写入完毕")

  • 1
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值