忍者ブログ

崗上虜囚の備忘録

日本よ!。私の日本への思いです。 コメントに返事を書かないこともあります。悪しからず。 コメントの投稿は日本人だけにしてください。 日本人でない場合は、国籍を書いてください。 注、google chromeで閲覧出来ませんので、filefoxかinternet explorerで閲覧してください

[PR]

×

[PR]上記の広告は3ヶ月以上新規記事投稿のないブログに表示されています。新しい記事を書く事で広告が消えます。

コメント

現在、新しいコメントを受け付けない設定になっています。

地震波形表示プログラム1

初めに
地震波形表示プログラムを作成したのでソースを掲載する。
動作環境はgccが使えることと、GMT(注1)と防災科学技術研究所(以下防災科研とする)のWIN32(注2)ツールが使えてパスが通っていることである。従ってlinuxならターミナル(DOS窓)からコマンドを入力しなければならない。
Windowsの場合はcygwinをインストールしておく必要がある。
 
注1GMTとは地図を描くためのオープンソースのソフトである。地図だけなくグラフも描けるので地震波形表示プログラムで使用した。波形表示出力はポストスクリプクト・ファイルなのでポストスクリプト・ファイルを表示出来るプログラム(winidowsの場合は、 GhostScript +GSViewGSView 等)も必要である。使用したGMTのバージョンは4.5.18であるが、他のバージョンで動くか不明である。
 
GMT 4.5.18のダウンロード先は、
日本語マニュアルは、
等がある。
注2WIN32ツールは防災科研が公開している地震波形データを取り扱うための処理プログラムである。
入手先は にあるが、ダウンロードするにはユーザー登録をしておかなければならない(波形データのダウンロードも同じく)。
WIN32ツールは適当な場所(例えば\cygwin\usr\local\w32 とか)にダウンロードして解凍したら、そこのディレクトリのターミナル窓で
------------------------------
$ make
------------------------------
と入力すれば使えるようになる筈である。
ついでに東大のwinシスーテム・ツールもインストールした方が良いだろう。ダウンロード先は、
↓ を参照。
 
地震波形表示プログラムの説明
今回作成した地震波形表示プログラムは、W32BPFtoPS.cW32IIRtoPS.c の2つである。
どちらもCで書かれており、実効ファイルにするには、gccでコンパイルする必要がある。
例えば
-----------------------------------------------------------
$ gcc -o W32BPFtoPS W32BPFtoPS.c -lm
$ gcc -o W32IIRtoPS W32IIRtoPS.c -lm
------------------------------------------------------------
これらのプログラムはCで書かれているが、動けば良いとの考え方で作ったので、余計な部分も有りあまり教育的ではないが、別プログラムを作成するときには参考にはなるだろう。
 
W32BPFtoPSについて
W32BPFtoPS は、防災科研からダウンロードしたwin32のデータファイル(****.cntの名称になっている)から、指定したID(4バイトのHEX)の波形データを抜き出してバンドパス・フィルタを掛けてポストスクリプト・ファイルに波形を描く。
コマンドの入力の仕方は、まず、W32BPFtoPS と入力すると、ヘルプ機能として以下の画面が出てくる。
 
W32BPFtoPS tbl w32 ps fl fh [t1] [tw] [s]
  tbl : chanel tbl file(in)
  w32 : win32 file     (in)
  ps  : ps file       (out)
  fl    : hi  pass filter file[Hz]: fl=0 fh=0 then no filter
  fh   : low pass filter file[Hz]: fl=x fh=y then Banded pass filter
  t1  : first time for display[sec]
  tw  : time width for display[sec]
  s   : if S then do not erase the work file
 
引数の意味:
tbl は目的のIDが記述されたチャネル・テーブル・ファイル名で、防災科研から波形データをダウンロードした時、"01_**_******.sjis.ch"と"01_**_******.euc.ch"の2つのファイルが添付されてくるが、これが全地点のチャネル・テーブル・ファイルで、そこから目的のIDが記述された箇所を新たなテキストファイルにコピーして(注3)、適当なファイル名にしたものである。
W32 はwin32の波形データファイルの名。
ps  は作成するポストスクリプト・ファイルの名。
fl   はバンドパス・フィルタのローカット周波数[Hz]。
fh  はバンドパス・フィルタのハイカット周波数[Hz]。
     fl と fh を0にした場合、バンドパス・フィルタ無しの生波形描画となる。
t1  波形データファイルの先頭時刻のt1秒後から波形を描画する。(省略した場合は先頭から)
tw  tw が指定された場合は、t1からtw秒間の波形を描画する。(省略した場合は最後まで)
s 
W32BPFtoPSは使用したテンポラリ・ファイルを削除して終了するが、Sと記入した場合は
  
テンポラリ・ファイルを削除しない。
  そこで "XXX.sh.eu" ファイルを編集して、sh 
XXX.sh.eu と入力すれば、編集した
  とおりのポストスクリプトファイルが出来る。
 
注3,チャネル・テーブルの内容の例(
北川地点と日向地点)は以下のようになっている。
-------------------------------------------------------------------------------------------------------
# n.kigh Kitagawa myzh02
6c83  1  0  N.KIGH  U     6   27   170.80  m/s     1.00  0.70  0  1.023e-07  32.6972 131.6829  -192  0  0  Kitagawa
6c84  1  0  N.KIGH  N     6   27   170.80  m/s     1.00  0.70  0  1.023e-07  32.6972 131.6829  -192  0  0  Kitagawa
6c85  1  0  N.KIGH  E     6   27   169.00  m/s     0.99  0.70  0  1.023e-07  32.6972 131.6829  -192  0  0  Kitagawa
# n.hygh Hyuga myzh15
6e43  1  0  N.HYGH  U     6   27   171.20  m/s     1.02  0.70  0  1.023e-07  32.3654 131.5893   -28  0  0  Hyuga
6e44  1  0  N.HYGH  N     6   27   167.00  m/s     1.02  0.70  0  1.023e-07  32.3654 131.5893   -28  0  0  Hyuga
6e45  1  0  N.HYGH  E     6   27   177.00  m/s     1.00  0.70  0  1.023e-07  32.3654 131.5893   -28  0  0  Hyuga
------------------------------------------------------------------------------------------------------------------
  
例1,
---------------------------------------------------------------------------------------------------------------
$ W32BPFtoPS KIG_HYGch.tbl 202201220108m15.w32 KIG_HYG.ps 0 0 40 60
----------------------------------------------------------------------------------------------------------------
意味、202201220108m15.w32の名称のwin32ファイルから、宮崎2地点が記述されたチャネルテーブルのチャネルのデータを抜き出し。先頭40秒から60秒間の生波形をKIG_HYG.psに描画する。
 
注、***.cntファイルは1分ファイルなので、202201220108m15.w32は15個の***.cntファイルを結合したものである。結合するには"catwin_32"か"CatCnts2w32”コマンドを使用する必要がある。CatCnts2w32コマンドについては別に説明する。
 
結果、
例2,
--------------------------------------------------------------------------------------------------------------------
$ W32BPFtoPS KIG_HYGch.tbl 202201220108m15.w32 KIG_HYG2.ps 0.01 1.0 40 60
---------------------------------------------------------------------------------------------------------------------
上記と同様にから先頭40秒から60秒間の波形に0.01〜1.0Hzのバンドパスフィルタを掛けKIG_HYG2.psに描画する。
 
結果、

W32BPFtoPS.c Cソース
ダウンロード


CatCnts2w32コマンドについて:
CatCnts2w32コマンドは複数のwin32ファイルを繋げて一つのwin32ファイルにするものであるこれもソースファイルを以下に示すので、テキストファイルにコピーしてgccでコンパイルし、使えるようにしたら良い。
 
使用例:
-------------------------------------------------------------------
CatCnts2w32 2022012201080103BM.cnt 6  202201220108m6.w32
-------------------------------------------------------------------
意味、2022012201080103BM.cnに続く6つのcntファイルを繋げた202201220108m6.w32のファイルを作成する。
 

CatCnts2w32のCソース:
/*********************************************************************/
//  << CatCnts2w32.c >> : 複数のcntファイル1を一つのw32ファイルにする
/*********************************************************************/
#include    <stdio.h>
#include    <time.h>
#include    <math.h>
#include    <string.h>
#include    <stdlib.h>
 
int     AtoToI(char *dec, int byt);     // 
time_t  IncDayTim(char *dDayTim, char *sDayTim);    // 
int     SplitLin(int max_n,char *bf,char *argv[]);
long    XtoL(char *dt);
void    CharCpy(char *dbf, char *sbf, int n);   // キャラクタ コピ-
 
/*********************************************************************/
int main(int argc ,char *argv[])
{
    FILE    *fp0,*fp2;
    char    cntFnm[360][32],winFnm[64],daytim[16];
    char    bf0[256],bf1[256];
    int     i,mint;
 
    if((argc<3)||(0==strcmp("--h",argv[1]))){
        printf("CatCnts2w32 cnt_file minute [win32_file]\n");
        printf("  cnt_file : first win32file (*****.cnt)\n");
        printf("  minute   : 1..60 (default: minute=1)\n");
        printf("  win32_file: win32 file name\n");
        printf("  .. 2008/10/14 ..\n");
        exit(0);
    }
 
    if(0==(mint=atoi(argv[2]))){ printf("syntax error\n"); exit(0); }
 
    for(i=0;i<mint;i++){                            // get 先頭 cnt ファイル名
        strncpy(cntFnm[i],argv[1],31); cntFnm[i][31]=16;
    }
 
    if(argc>=4){    strcpy(winFnm,argv[3]); }
    else{           strcpy(winFnm,argv[1]); }
 
    CharCpy(daytim,cntFnm[0],12); daytim[12]=0;     // get 先頭時刻
 
    for(i=1;i<mint;i++){                            // 次の1分時刻
        IncDayTim(daytim, daytim);                  // set 次のcnt ファイル名
        CharCpy(cntFnm[i],daytim,12);
    }
    
    sprintf(bf0,"cp %s tmp0_w32.tmp\n",cntFnm[0]);
    if(0==(fp0=popen(bf0,"r"))){    printf("error : %s\n",bf0); exit(0); }
    do{
        if(0==fgets(bf1,255,fp0)){  break; }
        if(0==strncmp("***** ERR",bf1,9)){ printf("error %s\n",bf1); exit(0); }
    }while(1);
    pclose(fp0);
    
    for(i=1; i<mint; i++){
        // -- cat_win32
        sprintf(bf0,"catwin32 tmp0_w32.tmp %s > tmp1_w32.tmp\n",cntFnm[i]);
        printf("%s",bf0);
        if(0==(fp0=popen(bf0,"r"))){    printf("error : %s\n",bf0); exit(0); }
        do{
            if(0==fgets(bf1,255,fp0)){  break; }
            if(0==strncmp("***** ERR",bf1,9)){ printf("error %s\n",bf1); exit(0); }
        }while(1);
        pclose(fp0);
        sprintf(bf0,"cp tmp1_w32.tmp tmp0_w32.tmp\n");
        if(0==(fp0=popen(bf0,"r"))){    printf("error : %s\n",bf0); exit(0); }
        do{
            if(0==fgets(bf1,255,fp0)){  break; }
            if(0==strncmp("***** ERR",bf1,9)){ printf("error %s\n",bf1); exit(0); }
        }while(1);
        pclose(fp0);
        
    }
    sprintf(bf0,"cp tmp1_w32.tmp %s\n",winFnm);
    if(0==(fp0=popen(bf0,"r"))){    printf("error : %s\n",bf0); exit(0); }
    do{
        if(0==fgets(bf1,255,fp0)){  break; }
        if(0==strncmp("***** ERR",bf1,9)){ printf("error %s\n",bf1); exit(0); }
    }while(1);
    pclose(fp0);
    
    if(0==(fp0=popen("rm *.tmp\n","r"))){ printf("error : rm tmp1.w32\n"); exit(0); }
    do{     if(0==fgets(bf1,255,fp0)) break;    }while(1);
    pclose(fp0);
    
    printf("seiko\n");
    exit(0);
}
 
//=====================================================================
// 
int AtoToI(char *dec, int byt)      // 
{
    char bf[6];
    int d;
    
    strncpy(bf,dec,byt); *(bf+byt)=0;
    d=atol(bf);
    return(d);
}
 
//=====================================================================
// yyyymmddhhmm+1 
time_t IncDayTim(char *dTim, char *sTim)
{
    time_t          st,dt;
    struct tm       sT,*dT;
    int             ye,mo,dy,ho,mi,cy;
 
    sT.tm_isdst = -1;   
    sT.tm_sec = 0;
    sT.tm_min = (*(sTim+10) & 0x0F)*10+(*(sTim+11) & 0x0F);
    sT.tm_hour= (*(sTim+8)  & 0x0F)*10+(*(sTim+9) & 0x0F);
    sT.tm_mday= (*(sTim+6)  & 0x0F)*10+(*(sTim+7) & 0x0F);
    sT.tm_mon = (*(sTim+4)  & 0x0F)*10+(*(sTim+5) & 0x0F) - 1;
    sT.tm_year= (*(sTim+2) & 0x0F)*10+(*(sTim+3) & 0x0F)
            + (*sTim & 0x0F)*1000  +(*(sTim+1) & 0x0F)*100 -1900;
    
    st = mktime(&sT);
    dt = st + (time_t)60;
    dT = localtime(&dt);
 
    sprintf(dTim,"%04d%02d%02d%02d%02d",
        dT->tm_year+1900,dT->tm_mon+1,dT->tm_mday,dT->tm_hour,dT->tm_min);
    
    return(0);
}
 
//=====================================================================
// 1ラインの文字列をデリミタ毎分割する                03/06/19
int SplitLin(int max_n,char *bf,char *argv[])
{
    int     n;
    char    *cp;
 
    n=0;
    cp=bf;
    argv[0]=cp;
    while(*cp!=0){
        switch(*cp){
         case 0:        return(n);
         case 0x0d:
            *cp=0;      return(n);
         case 0x0a:
            *cp=0;      return(n);
 
         case ' ':                          // 先頭がスペースなら
            while(*cp==' '){                //   スペースが続くかぎり
                *cp=0;                      //   0を書込
                cp++;                       // ポインタを++
                if(*cp==0) return(n);
            }
            break;
 
         case ',':                          // カンマなら
            while(*cp==','){
                *cp=0; cp++;                //  0を書く、ポインタを++
                if(*cp!=',') break;
                argv[n]=cp;  n++;           //  アドレスを保存、ポインタを++
                if(n>=max_n) return(n);
                if(*cp==0) return(n);
            }
            break;
 
         default:                           // それ以外なら
            argv[n]=cp; n++;                //   先頭アドレスを保存、ポインタを++
            if(n>=max_n) return(n);
            while((*cp!=' ')&&(*cp!=',')){
                if((*cp==0x0a)||(*cp==0x0d)){
                    *cp=0; return(n);
                }
                if(*cp==0) return(n);
                cp++;
            };
            break;
        }
    }
    return(n);
}
 
//-------------------------------------------------------------------------------
long XtoL(char *dt)
{
    int     i,j,l,n;
    long    a,c;
 
    l=strlen(dt);
    a=0; n=0;
    for(i=l-1;i>=0;i--){
        c=*(dt+i);
        if(c<'0') return(-1);
        if(c<='9'){
             c =c & 0x0f;
        }else{
            c =c & 0x0f;
            c =c + 9;
        }
        for(j=0;j<n;j++) c = c << 4;
        a = a + c;
        n++;
    }
    return(a);
}
//-------------------------------------------------------------------------------
void CharCpy(char *dbf, char *sbf, int n)   // キャラクタ コピ-
{
    int     i;
    
    for(i=0;i<n;i++)    *dbf++=*sbf++;
    return;
}
/*********************************************************************/
ここまで。
PR

コメント

お名前
タイトル
文字色
メールアドレス
URL
コメント
パスワード Vodafone絵文字 i-mode絵文字 Ezweb絵文字

カレンダー

03 2024/04 05
S M T W T F S
1 2 3 4 5 6
7 8 9 10 11 12 13
14 15 16 17 18 19 20
21 22 23 24 25 26 27
28 29 30

フリーエリア

最新トラックバック

プロフィール

HN:
崗上虜囚
性別:
非公開

アナライズ

カウンター