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;

Friday, December 6, 2013

TEXT2SEGY Header

Kode Perl di bawah ini berfungsi untuk meng-injeksikan nilai header yang terdapat dalam file ASCII (.txt) ke dalam header SEGY. 

Berikut adalah feature dari Kode Perl tersebut:
  • Nilai header yang akan di-injeksikan dimulai dari kolom pertama sampai kolom ke-n.
  • User diberikan kebebasan untuk mendefinisikan banyaknya header yang akan di-injeksikan.
  • Urutan header harus sama dengan urutan trace i.e. baris pertama pada file ASCII untuk trace pertama, dst.
  • File ASCII harus integer e.g. jika anda memiliki nilai floating point dengan dua digit dibelakang koma anda harus mengkalikannya dengan 100. Gunakan AWK untuk komputasinya.

Gambar di bawah ini menunjukkan dump dari header sebelum injeksi text.


Gambar di bawah ini menunjukkan 5 baris pertama dari header.txt (pada contoh ini saya memiliki 4 kolom). Juga, perintah text2segyhdr.pl dengan meng-injeksikan kolom pertama dari text file untuk byte 1, kolom 2 untuk byte 9, kolom 3 untuk byte 17, dan kolom 4 untuk byte 25. Perhatikan bagaimana saya menjalankan perintah Perl!


Gambar di bawah  ini adalah dump setelah injeksi. Perhatikan byte 001, 009, 017, 025 sudah berubah.


Copy kode berikut lalu berinama text2segyhdr.pl, jalankan spt contoh saya di atas.
#!/usr/bin/perl
#SEGY READ
$fname=$ARGV[0];
my $num_args = $#ARGV -2;
print ">>>>> Read SEGY...<<<<< \n";
sysopen(IN,$fname,O_RDONLY) or die "Can't Open $fname: $!";
$size= -s IN;
sysread (IN,$binary,$size);
close(IN);
$noofsamples = unpack("n",substr($binary,3220,4));
$bytewidth=4;  #32 bits
$tracecount=($size-3600)/(240+$bytewidth*$noofsamples);
$tracecount = sprintf("%.0f", $tracecount);
print ">>>>> Read TEXT HEADER... <<<<< \n";

#####ASCII TEXT HEADER READ
my $file = $ARGV[1];
my @data;
open( FILE, $file ) or die "Can't open file '$file': $!";
while( <FILE> ) {
chomp;
my @row = split;
push @data, \@row;
}
close( FILE );

###INSERT TEXT TO HEADER
for ($j=0; $j<=$num_args; $j++)
{
for ($i=0; $i<=$tracecount; $i++)
        {$startbyte=3600+($i)*240+($i)*$bytewidth*$noofsamples+$ARGV[$j+2]-1;
         sysopen(IN,$fname,2) or die "Can't sysopen $fname";
         sysseek(IN,$startbyte,0);
         $databin = pack 'N', "$data[$i][$j]";
         syswrite(IN,$databin,4) or die "Can't write $fname";
         close(IN);
print "\e[?25l";
printf ">>>>> Writing Progress...be patient... (%.0f%%)\r",$i/$tracecount*100 ;
}
}
print "\e[?25h";
print "\n";
print ">>>>> Done ! <<<<< \n";

Monday, December 2, 2013

Perl (Part-1)

Perl adalah bahasa program yang sangat populer digunakan baik untuk numerik, graphics, games, dll.

Pada bagian ini saya akan menjelaskan penggunaan Perl yang sering saya gunakan dalam perkerjan sebagai seorang Geoscientist.
Perl dapat diinstall dalam Plaform Linux, Mac maupun Windows. Anda dapat mengikuti petunjuk instalasi di sini.

Untuk Windows misalnya, kita bisa mengunakan Strawberry Perl dan Padre untuk Text Scripting.

Berikut adalah tampilan scripting pada Padre:

Untuk menjalankan program Perl:
Windows:
Jika anda menggunakan Padre cukup menggunakan F5
Command Windows perl namaprogram.pl

Linux:
./namaprogram.pl
pp -o namaprogram namaprogram.pl lalu jalankan ./namaprogram
perlcc namaprogram.pl

1. Hello World
#!/usr/bin/perl
use strict;
use warnings;
use 5.010;

say "Hello World";
print "Hello World\n"


2. Membaca matrix txt dan print pada console
#!/usr/bin/perl
use strict;
use warnings;
use 5.010;

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

#Print output
for my $i (0..$#data)
{
for my $j (0..$#{$data[0]})
{
   print "$data[$i][$j]","\t";
}
   print "\n";
}


3. Membaca dan Menulis 2D Matrix
#!/usr/bin/perl
use strict;
use warnings;
use 5.010;

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

open OUT," > matrixout.txt" or die "$!\n";
#Print output
for my $i (0..$#data)
{
for my $j (0..$#{$data[0]})
{
   print OUT "$data[$i][$j]","\t";
}
   print OUT "\n";
}


4. Matrix Transpose
#!/usr/bin/perl
use strict;
use warnings;
use 5.010;

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

#Print output
for my $j (0..$#{$data[0]})
{
for my $i (0..$#data)
{
print "$data[$i][$j]","\t";
}
print "\n";
}


5. Output setiap baris beberapa kali sesuai dengan integer yang terdapat pada kolom pertama:
#!/usr/bin/perl
use strict;
use warnings;
use 5.010;

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


#Print output
for my $i (0..$#data)
{
for (0..$data[$i][0]-1)      #[0] memilih integer pada kolom pertama
{
for my $j (0..$#{$data[0]})
{
   print "$data[$i][$j]","\t";
}
   print "\n";
}
}


6. Sort kolom 2
#!/usr/bin/perl
use strict;
use warnings;
use 5.010;

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


my @sorted = sort { $a->[1] <=> $b->[1] } @data;  ##sort by col 2
#Print output
for my $i (0..$#sorted)
{
for my $j (0..$#{$sorted[0]})
{
   print "$sorted[$i][$j]","\t";
}
   print "\n";
}


7. Sort berdasarkan kolom 2 dan 3
#!/usr/bin/perl
use strict;
use warnings;
use 5.010;

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

my @sorted = sort { $a->[1] cmp $b->[1] || $a->[2] cmp $b->[2] } @data;  ##sort by col  2 then 3
#Print output
for my $i (0..$#sorted)
{
for my $j (0..$#{$sorted[0]})
{
   print "$sorted[$i][$j]","\t";
}
   print "\n";
}


8. Memilih baris
#!/usr/bin/perl
use strict;
use warnings;
use 5.010;

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


#Print output
for my $i (1..2)  #memilih baris ke 2 s/d 3
{
for my $j (0..$#{$data[0]})
{
   print "$data[$i][$j]","\t";
}
   print "\n";
}


9. Memilih kolom
#!/usr/bin/perl
use strict;
use warnings;
use 5.010;

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


#Print output
for my $i (0..$#data)
{
for my $j (2-1,4-1) #memilih kolom 2 dan kolom 4 untuk di-output
{
   print "$data[$i][$j]","\t";
}
   print "\n";
}


10. Perkalian Matrix dengan subroutine
#!/usr/bin/perl -w
use strict;
use warnings;
my (@firstfactor, @secondfactor, @product);

my $file1 = 'matrixin1.txt';
open( FILE1, $file1 ) or die "Can't open file '$file1': $!";
while( <FILE1> ) {
chomp;
my @row1 = split;
push @firstfactor, \@row1;
}
close( FILE1 );

my $file2 = 'matrixin2.txt';
open( FILE2, $file2 ) or die "Can't open file '$file2': $!";
while( <FILE2> ) {
chomp;
my @row2 = split;
push  @secondfactor, \@row2;
}
close( FILE2 );


@product      = getproduct  (\@firstfactor, \@secondfactor);
for my $i (0..$#product)
{
for my $j (0..$#{$product[0]})
{
   print "$product[$i][$j]","\t";
}
   print "\n";
}

sub getproduct {
my ($a, $b) = @_;
my @product;
my $arow=$#$a;
my $acol=$#{$$a[0]};
my $brow=$#$b;
my $bcol=$#{$$b[0]};
  $acol == $brow or die "The matrices can't be multiplied!\n";
  for my $i (0..$arow) {
    for my $j (0..$bcol) {
      for my $k (0..$acol) {
        $product[$i][$j] += $$a[$i][$k] * $$b[$k][$j];
      }
    }
  }
  return @product;
}


11. Compare files (1)
more file1.txt
1 12345
2 26789
1 4567
5 567
3 987

more file2.txt
5 525
7 728
1 670
1 4567
1 12345

Output
1    12345    1    12345
1    4567    1    4567


#!/usr/bin/perl
$file1 = 'file1.txt';
$file2 = 'file2.txt';

open (A, $file1) || die ("Could not open $file!");
open (B, $file2) || die ("Could not open $file!");
while ($Aline = <A>)
{
        chomp $_;
        my @A = split(/\s+/,$Aline);
        seek B,0,0;
        while ($Bline = <B>)
        {
                chomp $_;
                my @B = split(/\s+/,$Bline);
                        if($A[0] eq $B[0] and $A[1] eq $B[1] ) {
                            print "$A[0]\t$A[1]\t$B[0]\t$B[1]\n";
                        }
        }
}
close (B);
close (A);