Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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 is IBM binary file. If I read it as 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 trace, and it is correctly read by mathematica, which is {0, 20, 40, 60, ..., 72320} feet. I don't need to read all the trace headers. What I care most is to read the correct velocity values.

Could someone help me? Thanks a lot!

(* The following code do not work right away, you need to download and unzip the file first and replace the file string path, the file name is timodel_vp.segy.gz, it is about 22 Megabytes after unzip and absolutely free to download. The address is: http://software.seg.org/datasets/2D/Hess_VTI/

Sorry for the inconvenience! *)

nByteFile = 
  FileByteCount[
   "e:/SeismicDataset/HessVTI/VelocityModelFiles/timodel_c11.segy"];

nByteTrace = (nByteFile - 3600)/3617;

nTraceSample = (nByteTrace - 240)/4;


str = OpenRead[
   "e:/SeismicDataset/HessVTI/VelocityModelFiles/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; 
share|improve this question

1 Answer 1

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
share|improve this answer
    
That is great help! Thanks a lot! –  yanfyon Jun 24 '13 at 18:20
    
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! –  yanfyon Jun 24 '13 at 18:26
1  
@yanfyon The new vectorized, compiled version should help with speed. –  Michael E2 Jun 24 '13 at 18:53
1  
@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. –  Oleksandr R. Jun 24 '13 at 18:53
    
You are really genius!!! Thank you so much, Michael. –  yanfyon Jun 24 '13 at 19:00

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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