Monday, December 23, 2013

ASCII2SEGY

Kode Perl di bawah ini berguna untuk mengkonversi data ascii menjadi segy.
1. Jumlah kolom merupakan jumlah trace
2. Jumlah baris merupakan jumlah sampel setiap trace
3. Sampling interval ditentukan oleh user pada command line
4. Data disimpan dalam format IEEE floating point single precision (4 bytes)

Setiap aplikasi seperti Seismic Unix, Petrel, Tomoplus, dll. memiliki konvensi sendiri dalam format code byte 3225-3226. Pada kode di bawah saya menggunakan format code=5 (IEEE). Berikut adalah format code yang umum digunakan:
1 - IBM floating point (4 bytes)
2 - fixed point (4 bytes)
3 - fixed point (2 bytes)
4 - fixed point w/gain code (4 bytes)
5 - IEEE floating point single precision (4 bytes)

Untuk menghindari kesalahan digit dan konversi, sebaiknya data ASCII memiliki nilai puluhan, ratusan atau ribuan. Contoh nilai amplitudo 0.0023 maka saya kalikan 10000 sehingga diperoleh 23.

Anda dapat melakukan scaling ASCII dengan:
awk '{for(i=1;i <=NF;++i)$i*=10000;print}' input.txt > output.txt

Kode awk di atas melakukan scaling ASCII*10000 untuk seluruh nilai amplitudo.

Copy kode Perl di bawah ini lalu beri nama ascii2segy.pl

Jalankan perintah:
./ascii2segy.pl file 4000

file adalah nama file ASCII, 4000(micro sec) adalah sampling rate.

Konversi ke SU:
segyread tape=file.sgy  verbose=1 endian=0 > file.su

Scale down dengan suhtmath:
suhtmath < file.su key=ns op=div scale=0 const=10000 > file_scale_back.su

Tampilkan:



Berikut adalah kode Perl:
#!/usr/bin/perl
my $num_args = $#ARGV;
die "usage: ascii2segy.pl filein samplingrate \n" if ($num_args+1 eq 0);
my $file = $ARGV[0];
$SR=$ARGV[1];

my @data;
#reading ascii
open( FILE, $file ) or die "Can't open file '$file': $!";
while( <FILE> ) {
chomp;
my @row = split;
push @data, \@row;
}
close( FILE );

$nooftrace=1+$#{$data[0]};
$noofsamples=1+$#data;

print "no of traces : $nooftrace\n";
print "no of samples: $noofsamples\n";
print "Sampling Interval: $SR\n";


open OUT," > $file.sgy" or die "$!\n";
###EBCDIC Header 3200 bytes
$ebcdic = pack("C3200",64,64,227,200,201,226,64,226,199,232,64,201,226,64,199,197,213,197,217,193,227,197,196,64,198,217,
214,212,64,193,226,195,201,201,37,195,64,64,64,64,194,232,64,193,199,228,226,64,193,194,196,228,211,211,193,200,
37,195,64,64,64,64,129,135,164,162,129,130,132,164,147,147,129,136,124,135,148,129,137,147,75,131,150,148);
print OUT "$ebcdic";

###Binary header 400 bytes
$bin3201 = pack("N","1");
$bin3205 = pack("N","1");
$bin3209 = pack("N","1");
$bin3213 = pack("n","$nooftrace");
$bin3215 = pack("n","0");
$bin3217 = pack("n","$SR");
$bin3219 = pack("n","$SR");
$bin3221 = pack("n","$noofsamples");
$bin3223 = pack("n","$noofsamples");
$bin3225 = pack("n","5"); #5(seismic unix) 6(tomoplus)
$bin3227 = pack("n","1");
$binres = pack("N93","1");
print OUT "$bin3201";
print OUT "$bin3205";
print OUT "$bin3209";
print OUT "$bin3213";
print OUT "$bin3215";
print OUT "$bin3217";
print OUT "$bin3219";
print OUT "$bin3221";
print OUT "$bin3223";
print OUT "$bin3225";
print OUT "$bin3227";
print OUT "$binres";

########################
my $databin;
for my $j (0..$#{$data[0]})
{
$k=$j+1;
$head1= pack("N60","$k");
print OUT "$head1";
for my $i (0..$#data)   
{
$nolsatu[$i][$j]=sprintf ("%b",unpack('L',pack('f',$data[$i][$j])));
if($data[$i][$j] > 0)
   {
$nol[$i][$j]=0;
$nol[$i][$j].= $nolsatu[$i][$j];
$databin = pack 'B32', "$nol[$i][$j]";
print OUT $databin;
   } else
   {
$databin = pack 'B32', "$nolsatu[$i][$j]";
print OUT $databin;
   }
}
progress_bar( $k, $nooftrace, 25, '=' );
}
print "\e[?25h";
print "\n";
print "...process done!: $file.sgy has been created\n";
sub progress_bar {
    my ( $got, $total, $width, $char ) = @_;
    $width ||= 25;
    $char  ||= '=';
    my $num_width = length $total;
    local $| = 1;
    print "\e[?25l";
    printf "|%-${width}s| write %${num_width}s traces of %s (%.2f%%)\r",
        $char x (($width-1)*$got/$total). '>', $got, $total, 100*$got/$total;
}

close OUT;

No comments: