12
$\begingroup$

I am trying to read a binary file. This file records seismic data. Basically its first 3600 bytes are for the general file information, following are for seismic traces, each begins with 240 bytes for information about this trace.

For the file I am reading, it has 3617 traces, each trace has 6240 bytes, with the first 240 bytes for trace headers. Following the trace headers are some velocity values for underground formation, it has 1500 data points. I can view this data using some commercial software... the values should be around 4999.99 feet/s(velocity of sea water) to 14800 feet/s. But what I read range from 590.125 to 743.25, which is obviously not correct.

Later I found out that the binary file I am trying to read stores IBM-format floating point values.

If I read it as an IEEE binary file, for integer data, it is all right, but for floating point data, what I read is totally wrong.

The trace head I read is about the geographic location of the traces, and it is correctly read by Mathematica, giving {0, 20, 40, 60, ..., 72320} feet. I don't need to read all the trace headers.

What I care about most is reading the correct velocity values.

Here's the code I have so far (the input files are from timodel_vp.segy.gz (22 MB) from here: http://software.seg.org/datasets/2D/Hess_VTI/)

nByteFile = 
  FileByteCount[
   "timodel_c11.segy"];

nByteTrace = (nByteFile - 3600)/3617;

nTraceSample = (nByteTrace - 240)/4;


str = OpenRead[
   "timodel_vp.segy", 
   BinaryFormat -> True];
pos = Range[3600, nByteFile - 1, nByteTrace];

traces = {SetStreamPosition[str, # + 76]; 
  BinaryRead[str, "Integer32", ByteOrdering -> +1], 
  SetStreamPosition[str, # + 240]; 
  Table[BinaryRead[str, "Real32", 
    ByteOrdering -> +1], {nTraceSample}]} & /@ pos; 
$\endgroup$

1 Answer 1

13
$\begingroup$

I imported the data this way:

data = Import["ftp://software.seg.org/pub/datasets/2D/Hess_VTI/timodel_vp.segy.gz", 
              "Binary"];

nByteFile = Length@data
nByteTrace = (nByteFile - 3600)/3617
nTraceSample = (nByteTrace - 240)/4
22573680
6240
1500

You have to convert bytes to IBM 32-bit floating-point yourself. This takes a list of bytes, divides into 4-byte groups and converts them.

ibmFloat32[bytes_List] := 
  With[{sign = BitGet[#[[1]], 7],
        exponent = BitClear[#[[1]], 7], 
        significand = FromDigits[#[[2 ;; 4]], 2^8]},
     N @ (-1)^sign * significand/2^24 * 16^(exponent - 64)] & /@ 
   Partition[bytes, 4];

Then we can find the traces. It takes about a minute. (The code could be made more efficient).

pos = Range[3600, nByteFile - 1, nByteTrace];

traces = {FromDigits[data[[# + 76 + 1 ;; # + 80]], 2^8],
          ibmFloat32[data[[# + 240 + 1 ;; # + nByteTrace]]]} & /@ pos;

Example

traces[[2]] // Short
{20, {4999.99, <<1498>>, 11878.}}

This is in the range specified, so I hope it's working.


This is almost 100 times faster, but less readable:

newibmFloat32 = Compile[{{bytes, _Integer, 1}},
   (2. UnitStep[127.5 - #[[1]]] - 1) *
   16^(#[[1]] - UnitStep[#[[1]] - 127.5] 128. - 64.) *
   {2.^-8, 2.^-16, 2.^-24}.#[[2 ;; 4]] &@
      Transpose@Partition[N@bytes, 4]
   ];

traces2 = {FromDigits[data[[# + 76 + 1 ;; # + 80]], 2^8],
           newibmFloat32[data[[# + 240 + 1 ;; # + nByteTrace]]]} & /@ pos;

traces2 == traces
True
$\endgroup$
6
  • $\begingroup$ That is great help! Thanks a lot! $\endgroup$
    – yanfyon
    Jun 24, 2013 at 18:20
  • $\begingroup$ Thanks a lot! Just that I am a bit greedy. The original code takes about 19 seconds to read IEEE binary file, and the core code above takes about 84 seconds to run in my PC. Often this kind of file is bigger than 1 Gigabytes, so speed is also an important consideration. Thanks again for Michael's great help! $\endgroup$
    – yanfyon
    Jun 24, 2013 at 18:26
  • 1
    $\begingroup$ @yanfyon The new vectorized, compiled version should help with speed. $\endgroup$
    – Michael E2
    Jun 24, 2013 at 18:53
  • 1
    $\begingroup$ @yanfyon bit manipulation is a lot faster in C. If you need to convert large files, you could write the conversion routine as a LibraryLink library. Alternatively, you could use such a library to read the file and convert it before passing the contents to Mathematica. $\endgroup$ Jun 24, 2013 at 18:53
  • $\begingroup$ You are really genius!!! Thank you so much, Michael. $\endgroup$
    – yanfyon
    Jun 24, 2013 at 19:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.