Fortran Coder

查看: 15192|回复: 2
打印 上一主题 下一主题

[文件读写] Fortran批量读取SAC地震数据文件,写SAC文件,并转换成SEGY

[复制链接]

13

帖子

9

主题

0

精华

入门

F 币
75 元
贡献
52 点
跳转到指定楼层
楼主
发表于 2016-2-28 20:04:29 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
[Fortran] 纯文本查看 复制代码
Module SACParameter
    Integer, parameter :: FileUnit = 100
    Character( Len = * ), Parameter :: Directer = '"'//'..\WAV\WAVFormYiliang\yiliang\*.*'//'"'
    Type :: EventPara
        Integer( kind = 4 ) :: StationNumber
        Integer( kind = 4 ) :: FileNumber
        Character( Len = 4 ), allocatable :: StationName(:)        ! 台站名数组   
        Character( Len = 4 ), allocatable :: ComponentName(:,:)        ! 三分量文件名(台站,分量)
        Character( Len = 1 ), allocatable :: stringFileName(:,:)    ! 文件名信息数组二维字符数组
        Real(kind = 4), allocatable :: StationWaveformLong(:,:)    ! 重采样后的波形长度(台站,分量)
        Real(kind = 4), allocatable :: StationMaximumAmplitude(:,:)! 归一化后的最大振幅(台站,分量) 
        Real(kind = 4), allocatable :: StationDELTA(:,:)           ! 重采样采样间隔(台站,分量)
    End Type EventPara
    Type :: SACHeader
        Real( Kind = 4 ), Allocatable :: DELTA(:,:), DEPMIN(:,:), DEPMAX(:,:), SCALE(:,:), ODELTA(:,:), B(:,:), E(:,:),                                                                    &
                             O(:,:), A(:,:), INTERNAL1(:,:), T0(:,:), T1(:,:), T2(:,:), T3(:,:), T4(:,:), T5(:,:), T6(:,:), T7(:,:), T8(:,:), T9(:,:), F(:,:),                             &
                             RESP0(:,:),RESP1(:,:),RESP2(:,:),RESP3(:,:),RESP4(:,:),RESP5(:,:),RESP6(:,:),RESP7(:,:),RESP8(:,:),RESP9(:,:),                                                &
                             STLA(:,:),STLO(:,:),STEL(:,:),STDP(:,:),EVLA(:,:),EVLO(:,:),EVEL(:,:),EVDP(:,:),MAG(:,:),                                                                     &
                             USER0(:,:), USER1(:,:), USER2(:,:), USER3(:,:), USER4(:,:), USER5(:,:), USER6(:,:), USER7(:,:), USER8(:,:), USER9(:,:),                                       &
                             DIST(:,:), AZ(:,:), BAZ(:,:), GCARC(:,:), INTERNAL2(:,:), INTERNAL3(:,:),                                                                                     &
                             DEPMEN(:,:), CMPAZ(:,:), CMPINC(:,:), XMINIMUM(:,:), XMAXIMUM(:,:), YMINIMUM(:,:), YMAXIMUM(:,:),                                                             &
                             ADJTM(:,:), UNUSED1(:,:), UNUSED2(:,:), UNUSED3(:,:), UNUSED4(:,:), UNUSED5(:,:), UNUSED6(:,:)
        Integer( Kind = 4 ), Allocatable :: NZYEAR(:,:), NZJDAY(:,:), NZHOUR(:,:), NZMIN(:,:), NZSEC(:,:), NZMSEC(:,:), NVHDR(:,:),                                                        &
                             NORID(:,:), NEVID(:,:), NPTS(:,:), NSPTS(:,:), NWFID(:,:), NXSIZE(:,:), NYSIZE(:,:), UNUSED7(:,:),                                                            &
                             IFTYPE(:,:), IDEP(:,:), IZTYPE(:,:), UNUSED8(:,:), IINST(:,:), ISTREG(:,:), IEVREG(:,:), IEVTYP(:,:), IQUAL(:,:), ISYNTH(:,:), IMAGTYP(:,:), IMAGSRC(:,:),    &
                             UNUSED9(:,:), UNUSED10(:,:), UNUSED11(:,:), UNUSED12(:,:), UNUSED13(:,:), UNUSED14(:,:), UNUSED15(:,:), UNUSED16(:,:),                                        &
                             LEVEN(:,:), LPSPOL(:,:), LOVROK(:,:), LCALDA(:,:), UNUSED17(:,:)
        Character( Len = 16 ), Allocatable :: KEVNM(:,:)
        Character( Len = 8 ), Allocatable :: KSTNM(:,:), KHOLE(:,:), KO(:,:), KA(:,:), KT0(:,:), KT1(:,:), KT2(:,:), KT3(:,:), KT4(:,:), KT5(:,:), KT6(:,:), KT7(:,:), KT8(:,:), KT9(:,:), &
                             KF(:,:), KUSER0(:,:), KUSER1(:,:), KUSER2(:,:), KCMPNM(:,:), KNETWK(:,:),  KWaveformRD(:,:), KINST(:,:)

        Integer(kind = 4), Allocatable :: Header1(:,:)     ! 头端一(台站,分量)
        Real(kind = 4), Allocatable :: Header2(:,:)        ! 头端二(台站,分量)
        Character(Len = 8), Allocatable :: Header3(:,:)    ! 头端三(台站,分量)
        real( kind = 4 ), Allocatable :: Waveform(:,:,:)   ! 波形(台站,分量,记录)
        Integer(kind = 4), Allocatable :: WaveformLong(:,:)             ! 波形长度(台站,分量)
    End Type SACHeader 
End module SACParameter

Program Main
USE SACParameter
  Implicit None
    type(EventPara) :: EVENT1
    type(SACHeader) :: SACFILE1
  External WriteName, SplitName

  !Character( Len = * ), Intent( In ) :: DirecterPath
  Integer :: i, j, k, m, n
  !Character( Len = * ), Parameter :: Path = '"'//'DirecterPath'//'"'
  Call DoWithWildcard( '..\WAV\WAVFormYiliang\yiliang\*.*' , WriteName , n )
  If ( N >= 0 ) then
    Write(*,*) '共' , N , '个文件!'
  End If
  pause 1
  Call DoWithWildcard( '..\WAV\WAVFormYiliang\yiliang\*.*' , SplitName , n )
  Do i = 1, EVENT1%StationNumber
    Write(*,*) EVENT1%StationName(i)
    Do j = 1, EVENT1%FileNumber
        Write(*,*) EVENT1%ComponentName(i,j)
    End do
  End do
End Program Main


Subroutine DoWithWildcard(cWildcard,CallBack,iTotal)
  Use DFLib,only:GetFileInfoQQ,GetLastErrorQQ,FILE$INFO,FILE$LAST,FILE$ERROR,FILE$FIRST,ERR$NOMEM,ERR$NOENT,FILE$DIR
  USE SACParameter
  Implicit None
    type(EventPara) :: EVENT1
    type(SACHeader) :: SACFILE1
  Interface 
    Subroutine CallBack( FileName , loop )
      Character(*),Intent(In) :: FileName
      Integer,Intent(In) :: loop
    End Subroutine CallBack
  End Interface
  Character*(*),Intent(In)::cWildcard
  Integer,Intent(Out)::iTotal
  TYPE (FILE$INFO) info
  INTEGER(4)::Wildhandle,length,retInt
  Wildhandle = FILE$FIRST
  iTotal = 0
  DO WHILE (.TRUE.)
      length = GetFileInfoQQ(cWildCard,info,Wildhandle)
      IF ((Wildhandle .EQ. FILE$LAST) .OR.(Wildhandle .EQ. FILE$ERROR)) THEN
        SELECT CASE (GetLastErrorQQ())
        CASE (ERR$NOMEM)  !//内存不足
          iTotal = - 1
          Return
        CASE (ERR$NOENT)  !//碰到通配符序列尾
          Return
        CASE DEFAULT
          iTotal = 0
          Return
        END SELECT
      END IF
      If ((info%permit.AND.FILE$DIR).Eq.0) then
        Call CallBack( Trim(info.Name) , iTotal + 1 )
        iTotal = iTotal + 1
      End If
  END DO
  EVENT1%FileNumber = iTotal
  !Allocate(EVENT1%stringFileName(EVENT1%FileNumber))
End Subroutine DoWithWildcard

Subroutine WriteName( FileName , loop )
USE SACParameter
  Implicit None
    type(EventPara) :: EVENT1
    type(SACHeader) :: SACFILE1
  Character( Len = *),Intent(In) :: FileName
  Integer,Intent(In) :: loop
  Write(*,*) loop , FileName
End Subroutine WriteName

Subroutine SplitName( FileName , loop )          
USE SACParameter                                 
  Implicit None                                  
    type(EventPara) :: EVENT1                  
    type(SACHeader) :: SACFILE1                  
  Character(*),Intent(In) :: FileName             
  Integer,Intent(In) :: loop                    
  Integer( kind = 4) :: i, j, k, m = 0, n = 0, FileNameLong
  Character( Len = 1 ) :: separator = "."
  Character( Len = 4 ) :: Station
  FileNameLong = len_trim(FileName)
  j = Index( FileName, separator )
  k = Index( FileName, separator, .True. )
  If( Station /= FileName(1:j - 1)) then
      n = n + 1
      m = 1
      Station = FileName(1:j - 1)
      write(*,*) Station
      pause 2
      EVENT1%StationName(n) = FileName(1:j - 1)
      EVENT1%ComponentName(n,m) = FileName( k + 1:FileNameLong)
  Else
      EVENT1%ComponentName(n,m) = FileName( k + 1:FileNameLong)
      m = m + 1
  End if
  EVENT1%StationNumber = n
End Subroutine SplitName

yiliang.zip

718.69 KB, 下载次数: 24

SAC二进制波形文件

分享到:  微信微信
收藏收藏 点赞点赞 点踩点踩

2015

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1572 元
贡献
676 点

美女勋章热心勋章星光勋章新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

沙发
发表于 2016-2-28 21:40:29 | 只看该作者
你的问题是,所有的可分配数组都没有分配。

我对 SEGY 文件有一定了解,并写了一个关于 SEGY 读取的库。可惜是公司的,不方便给你。

我不清楚 SAC 地震数据格式,也不知道如何读取。

490

帖子

4

主题

0

精华

大宗师

F 币
3298 元
贡献
1948 点

水王勋章元老勋章热心勋章

板凳
发表于 2016-2-29 09:02:54 | 只看该作者
fcode 发表于 2016-2-28 21:40
你的问题是,所有的可分配数组都没有分配。

我对 SEGY 文件有一定了解,并写了一个关于 SEGY 读取的库。可 ...

github上面随便找找就有,为啥lz要执迷于Fortran呢?
传送门:
Matlab代码
Python代码
C#代码
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-11-1 09:33

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表