Python整理PEER所下载的地震波源数据——提取地震波至txt+生成地震波反应谱

介绍

在进行动力时程分析前,选取合适的地震波是必须的。当前最热门的地震波下载网站当属PEER(太平洋地震工程研究中心)。从此网站下载地震波以后,往往需要将地震波数据进行二次转化,一方面得到用于导入其他软件的特种格式,另一方面生成地震波的反应谱。得到反应谱并不难,但是如何及时得到反应谱而略掉繁琐的‘点、点、点~’呢?本篇文章将解决此问题。

gitee代码托管

全部功能实现代码已托管至我的gitee账号Hongtao_He,如有代码问题,可在gitee中提交,或与我联系。对于地震波数据整理所需的关键函数已全部展示于上,亦可根据提示编辑属于自己的实现程序。
创作不易,欢迎关注、点赞、转发、支持,谢谢!

1. 实现目标

本篇文章的前提是处理从PEER上下载的地震波,目标是将所下载的地震波的地震波数据(存于txt中)、dt(时间步长)、ntps(时间步)、加速度反应谱四类信息进行处理。得到的结果,将地震波数据按不同时间步长整理的文件夹;含有地震波名称、dt、npts、加速度谱信息的表格。

运行结果

1 已下载的文件

Python整理PEER所下载的地震波源数据——提取地震波至txt+生成地震波反应谱
为了展示,此处只下载了两条地震波,并在此目录文件内及下一级子目录文件内将地震波进行解压。可以看到,在这种状态下,当每级目录下的地震波若有100条时,或者,一条一条下载所得的地震波,处理起来的工程量不小。虽然能处理地震波的软件不少,但是地震波导入以后再依次复制出目标信息的方式非常的繁琐。

2 实现的效果

运行后,处理后的信息会汇总在“Processed Output Summary”文件夹中。
点击运行后,会将处理后的信息进行汇总于“处理后输出汇总”文件夹内。
“处理后输出汇总”文件夹,是根据地震波所在的目录,将地震波按目录重新创建文件夹,并将其按时间步长类型进行整理。“地震波信息汇总.xlsx”表格保存有每级目录下的地震波的名称、dt、npts、加速度谱。
Python整理PEER所下载的地震波源数据——提取地震波至txt+生成地震波反应谱
Python整理PEER所下载的地震波源数据——提取地震波至txt+生成地震波反应谱

二、代码流程

1 需要的库

import os
import glob
import openpyxl
from pyexcel.cookbook import merge_all_to_a_book as sss
import math

2 主要函数:共5个

Python整理PEER所下载的地震波源数据——提取地震波至txt+生成地震波反应谱

函数1:得到‘.py’程序所在的文件目录及下一级子目录的路径列表a,返回列表a

def file_name(file_dir): # 将目标路径和它下一级的子目录文件夹路径全部输出至列表
    a1 = 1
    a = []
    for root, dirs, files in os.walk(file_dir):
        if a1 == 1:
            b = dirs
            a.append(file_dir)
            for i in range(0,len(b)):
                h = file_dir + '\\' + b[i]
                a.append(h)
            break
    return a

函数2:将每一路径下的‘.csv’文件转为‘.xlsx’(‘.csv’文件存有地震波名称,此步转化方便将地震波名称从‘.xlsx’中提取出来)

def csv_to_xlsx(path_1):  # path_1为列表
    h1 = file_name(path_1)
    bb1 = 1
    print('共有' + str(len(h1)) + '个含有' + '.scv' + '文件夹目录')
    for i in range(0,len(h1)):
        print('现在处理第' + str(bb1) + '个文件夹内的所有'+ ''''.csv'文件!''')
        b1 = h1[i] + '\\' '*.csv'
        b2 = glob.glob(b1)  # 获取bb文件夹内的所有‘.csv’文件的路径
        print(b2)
        b3 = []
        if b2 == []:
            continue
        else:
            for j in range(0,len(b2)):
                b3.append(b2[j])
                sss(b3, h1[i] + '\\' + b2[j].split('\\')[-1].replace('.csv','.xlsx'))
                b3 = []
        bb1 += 1
    print('’.scv‘文件处理完毕')

函数3:处理地震波源数据,返回地震波的dt、npst、数据列表

def Ssw(inFilename):
    dt = 0.0
    npts = 0
    a = []
    inFileID = open(inFilename, 'r')
    flag = 0
    for line in inFileID:
        if line == '\n':
            continue
        elif flag == 1:  # 此处用flag来进行标记,读取到带有DT的行后,后面的行都标记为1(进而将后面的行进行写出)
            ###################################
            words = line.split()
            lengthLine = len(words)
            if words[0] == '***':
                break
            else:
                for i in range(0, lengthLine):
                    a.append(float(words[i]))  # 将地震波数据提取并存储至列表a中
            #######################################
        else:  # 此处为对数据进行读取,并标记开始进行数据写入的开始行(只不过并没有标记处是从哪一行开始,而是以从这一行开始后面都是需要写入的地震波数据
            words = line.split()
            lengthLine = len(words)
            if lengthLine >= 4:
                if words[0] == 'NPTS=':
                    for word in words:
                        if word != '':
                            if flag == 1:
                                dt = float(word)
                                break  # 此处的break中断的点是单个line的判断
                            if flag == 2:
                                npts = int(word.strip(','))
                                flag = 0
                            if word == 'DT=' or word == 'dt':
                                flag = 1
                            if word == 'NPTS=':
                                flag = 2
                elif words[-1] == 'DT':
                    count = 0
                    for word in words:
                        if word != '':
                            if count == 0:
                                npts = int(word)
                            elif count == 1:
                                dt = float(word)
                            elif word == 'DT':
                                flag = 1
                                break  # 此处的break中断的点是单个line的判断
                            count += 1
    return dt, npts, a

函数4:根据所给路径创建文件夹

def mkdir(path): # 此项为创文件
    import os
    # 去除首位空格
    path = path.strip()  # 删除空白符(包括'\n', '\r',  '\t',  ' ')
    path = path.rstrip("\\")  # 去除尾部 \ 符号
    # 判断路径是否存在
    isExists = os.path.exists(path)
    if not isExists:
        # 如果不存在则创建目录
        # 创建目录操作函数
        os.makedirs(path)
        print(path + ' 创建成功')
        return True
    else:
        # 如果目录存在则不创建,并提示目录已存在
        print(path + ' 目录已存在')
        return False

函数5:由保存有地震波的列表,计算地震波反应谱(此处参考Dino陈博所推荐的Nigam积分法)

def nigam_clac(k,m,dt,damp,nstep):
    tt=w=h=w2=hw=0.0
    hw=wd=e=a11=a12=a21=a22=ss=cc=s1=c1=s2=c2=s3=c3=0.0
    wdt=cwdt=swdt=b11=b12=b21=b22=x=sa=sv=sd=0.0
    ddx=dx=dxf=xf=ddym=ddyf=0.0
    nn=mm=0
    dis=[]
    vel=[]
    acc=[]
    tim=[]
    pi = math.pi
    ##
    tt =2*pi*math.sqrt(m/k)
    w =2*pi/tt   ##cir freq
    h =damp
    w2 =w*w
    hw =h*w
    wd =w*math.sqrt(1-h*h)
    wdt =wd*dt
    e =math.exp(-hw*dt)
    cwdt =math.cos(wdt)
    swdt =math.sin(wdt)
    ##
    a11 =e*(cwdt+hw*swdt/wd)
    a12 =e*swdt/wd
    a21 =-e*w2*swdt/wd
    a22 =e*(cwdt-hw*swdt/wd)
    ss =-hw*swdt-wd*cwdt
    cc =-hw*cwdt+wd*swdt
    s1 =(e*ss+wd)/w2
    c1 =(e*cc+hw)/w2
    s2 =(e*dt*ss+hw*s1+wd*c1)/w2
    c2 =(e*dt*cc+hw*c1-wd*s1)/w2
    s3 =dt*s1-s2
    c3 =dt*c1-c2
    b11 =-s2/wdt
    b12 =-s3/wdt
    b21 =(hw*s2-wd*c2)/wdt
    b22 =(hw*s3-wd*c3)/wdt
    ##
    acc.append(2*hw*gm[0]*dt)
    vel.append(-1*gm[0]*dt)
    dis.append(0.0)
    tim.append(0.0)
    x =0.0
    sa =0.0
    sv =0.0
    sd =0.0
    nn =nstep
    for mm  in range(nn):
        tim.append(mm*dt)
        dxf =dx
        xf =x
        ddym =gm[mm+1]
        ddyf =gm[mm]
        x = a12*dxf+a11*xf+b12*ddym+b11*ddyf
        dx =a22*dxf+a21*xf+b22*ddym+b21*ddyf
        ddx =-2*hw*dx-w2*x
        acc.append(ddx)
        vel.append(dx)
        dis.append(x)
        sa =max(sa,abs(ddx))
        sv =max(sv,abs(dx))
        sd =max(sd,abs(x))
    ##
    alpha =sa/9800.0
    kdisp =sd
    ##return tim,dis,vel,acc
    return alpha,kdisp

3 实现思路

1 从‘.py’文件所在目录,得到当前目录及下一级子目录的路径列表a1
2 遍历路径列表a1,并将其下的‘.csv’文件转化为‘.xlsx’,并保存在当前目录下
3 创建输出目录并得到其路径,遍历路径列表a1,先将先前生成的‘.xlsx’文件打开并得到指定方向的地震波名称列表a2,根据地震波名称列表a2,逐个处理地震波源数据文件,并将所得信息进行输出至所创建的输出目录中。
4 完成地震波文件的处理

文章出处登录后可见!

已经登录?立即刷新

共计人评分,平均

到目前为止还没有投票!成为第一位评论此文章。

(0)
心中带点小风骚的头像心中带点小风骚普通用户
上一篇 2022年3月16日 上午10:37
下一篇 2022年3月16日 上午10:53

相关推荐