Tuesday, November 19, 2013

SEGY2ASCII (PERL)

Pada bagian ini saya telah menjelaskan penggunaan SegMAT yang berbasis Matlab atau Octave untuk melakukan konversi file segy ke ascii. Namun penggunaan kedua software di atas sangat lambat untuk konversi data yang besar.

Konversi segy ke ascii sebenarnya bisa dilakukan dengan menggunakan  Seismic Unix sebagaimana yang telah saya jelaskan di sini. Tentu saja kita perlu meng-konversi dari segy ke su terlebih dahulu. Outputnya adalah satu kolom untuk seluruh trace seismik, oleh karena itu kita perlu menyusun kembali sehingga diperoleh susunan trace yang benar.

Kode dibawah ini adalah kode PERL untuk konversi segy ke ascii (32-bit IBM dan IEEE Single Precision).
  1. Copy kode tersebut dengan text editor lalu bernama misalnya segy2ascii.pl
  2. Eksekusi: ./segy2ascii.pl file.sgy pada mesin lain segy2ascii.pl file.sgy

#!/usr/bin/perl
my $num_args = $#ARGV;
die "usage: segy2ascii.pl file.sgy \n" if ($num_args+1 eq 0);

$fname=$ARGV[0];
sysopen(IN,$fname,O_RDONLY) or die "Can't open $fname: $!";
$size= -s IN;
open my $fh, '<', $fname or die "Can't open ";
$bytewidth=4;
my $cacing;
read($fh,$cacing,3600,0);
$noofsamples = unpack("n",substr($cacing,3220));
$tracecount=($size-3600)/(240+$bytewidth*$noofsamples);

$tracecount = sprintf("%.0f", $tracecount);
print "no of traces : $tracecount\n";
print "no of samples: $noofsamples\n";

print " Type 1 = IBM 32-Bits floating point or 2 = IEEE 32 bits floating point\n";
$formatcode= <STDIN>;
if($formatcode==1)
   {
   print " Converting $fname to ASCII with IBM 32-Bits floating point format...\n";
   } elsif($formatcode==2)
   {
   print " Converting $fname to ASCII with IEEE 32-Bits floating point format...\n";
   } else
   {
   print " Type 1 = IBM 32-Bits floating point or 2 = IEEE 32 bits floating point\n";
   exit;
   }

seek (IN,3840,1);
my $buffer;
my $i=1;
open OUT," > $fname.txt" or die "$!\n";
while (sysread(IN, $buffer, 240+($noofsamples*$bytewidth),0))
{
for ($j=0; $j<=$noofsamples-1; $j++)
{
$amp[$j]=unpack("B*",substr($buffer,$j*4));

if($formatcode==1)
   {
   $trace_value[$j] = conv_bit_string_2_ibm32float($amp[$j]);
   } elsif($formatcode==2)
   {
   $trace_value[$j] = conv_bit_string_2_ieee32float($amp[$j]);
   }

print OUT $trace_value[$j],"\t";
}
print OUT "\n";
progress_bar( $i, $tracecount, 25, '=' );
$i++;
}
print "\e[?25h";
print "\n";
print "...process done!:   $fname.txt (rows(traces)=$tracecount cols(samples)=$noofsamples) 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| writeread %${num_width}s traces of %s (%.2f%%)\r",
        $char x (($width-1)*$got/$total). '>', $got, $total, 100*$got/$total;
}
sub bin2dec {
return unpack("N", pack("B32", substr("0" x 32 . shift, -32)));
}

sub conv_bit_string_2_ibm32float {
$bit_string_ibm = shift;
$first_digit_ibm = substr($bit_string_ibm, 0, 1);
if ( $first_digit_ibm eq "0" ) {
$sign_bit_ibm = 1;
} elsif ( $first_digit_ibm eq "1" ) {
$sign_bit_ibm = -1;
}
$bin_exponent_ibm = substr($bit_string_ibm, 1, 7);

$exponent_ibm = bin2dec($bin_exponent_ibm);
$bin_fraction_ibm = substr($bit_string_ibm, 8, 31);
@bit_chars_ibm = unpack("A1" x length($bin_fraction_ibm), $bin_fraction_ibm);
$place_holder_ibm = -1;
$fraction_ibm = 0;
foreach $bit_ibm ( @bit_chars_ibm ) {
$fraction_ibm += $bit_ibm * (2 ** $place_holder_ibm);
$place_holder_ibm += -1;
}
$ibm_float_ibm = ($sign_bit_ibm ** 1) * (16 ** ($exponent_ibm - 64)) *($fraction_ibm);
return sprintf("%.4f", $ibm_float_ibm);
}

sub conv_bit_string_2_ieee32float {
$bit_string_ieee = shift;
$first_digit_ieee = substr($bit_string_ieee, 0, 1);
if ( $first_digit_ieee eq "0" ) {
$sign_bit_ieee = 1;
} elsif ( $first_digit_ieee eq "1" ) {
$sign_bit_ieee = -1;
}
$bin_exponent_ieee = substr($bit_string_ieee, 1, 8);
$exponent_ieee = bin2dec($bin_exponent_ieee);
$bin_fraction_ieee = substr($bit_string_ieee, 9, 31);
@bit_chars_ieee = unpack("A1" x length($bin_fraction_ieee), $bin_fraction_ieee);
$place_holder_ieee = -1;
$fraction_ieee = 0;
foreach $bit_ieee ( @bit_chars_ieee ) {
$fraction_ieee += $bit_ieee * (2 ** $place_holder_ieee);
$place_holder_ieee += -1;
}
$ieee_float_ieee = ($sign_bit_ieee ** 1) * (1+$fraction_ieee)*(2**($exponent_ieee - 127));
return sprintf("%.4f", $ieee_float_ieee);
}

exit;

No comments: