Message Boards Message Boards

Working with SEGY file format for storing geophysical data

GROUPS:

Preamble

Several days ago I found this question in the Community. I think it interesting theme for discussion. Because I want to show how you can using the package CustomImportExport. Something of this article will repeat comments on the previous question, but I want save this post to be a stand-alone tutorial for the package and SEGY in Wolfram Language.

Downloading

Clone or download this repository CustomImportExport using git (in the terminal or cmd-window):

git clone https://guthub.com/KirillBelovTest/CustomImportExport.git

Also you can download repository using green button on GitHub

download button on GitHub

Installing

After you can install the package in your system. You may open the file Installer.nb and execute this code:

SetDirectory[NotebookDirectory[]]; 
Get["Installer.m"]; 

Installer.m - script that will create a new folder in $UserBaseDirectory. After this you can call the context "CustomImportExport`" from any notebooks.

Reading SEGY data

Now let try read the .segy-file. For demonstration in the repository exist the file "/CustomImportExport/Resources/MarmousiModel.segy"

If you was installed the package, you may execute following in any notebook:

$HistoryLength = 0; 
<<CustomImportExport`
SetDirectory[$CustomImportExportDirectory]; 
marmousiPath =  FileNameJoin[{"CustomImportExport", "Resources", "MarmousiModel.segy"}]; 
marmousiData = CustomImport[marmousiPath, "SEGY"]

(*
  Out[..] := SEGYData[{SEGYElement["TextHeader"], SEGYElement[..], ..}]
*)

This output contains all data from .segy-file. equivalents of the last line of code:

marmousiData = CustomImport[File[marmousiPath], "SGY"]; 

Working with loaded data

Now we have the variable marmousiData. This is a object with type SEGYData. For this type in the package defined special operations. For example you can get same of elements:

marmousiData["TextHeader"]

(* Out[..] := SEGYElement["TextHeader", <3200 whitespaces>] *)

Similarly, you can get other elements:

  1. TextHeader
  2. BinaryHeader
  3. TraceHeaders
  4. Traces
marmousiData["BinaryHeader"]
(* Out[..] := SEGYElement["BinaryHeader", {"JobID" -> 0, "LineNumber" -> 0, ..}] *)

marmousiData["TraceHeaders"]
(* Out[..] := SEGYElement["TraceHeaders", {"tracl", ..} -> {{0, 1, ..}, ..}] *)

marmousiData["Traces"]
(* Out[..] := SEGYElement["Traces", {{0.0, 1.0, ..}, ..}] *)

I created an additional type for convenience. SEGYElement also has additional functionality:

marmousiTraceHeaders = marmousiData["TraceHeaders"]; 

marmousiTraceHeaders[]
(* Out[..] := {"tracl", ..} -> {{0, 1}, ..} 
returns list of rules - second element in the marmousiTraceHeaders*)

And this will work for all SEGYElement-items SEGYElement["Name", data][] - returns internal data of the object. In addition, the elements have other ways of getting data. For example - getting one of the key from the BinaryHeader (several ways):

marmousiBinaryHeader = marmousiData["BinaryHeader"]; 
marmousiData["BinaryHeader"]["JobID"]
marmousiData["BinaryHeader", "NumberDataTraces"]
marmousiBinaryHeader["NumberOfSamplesForReel"]
marmousiBinaryHeader[{"NumberDataTraces", "NumberOfSamplesForReel"}]

(*
  Out[..] := 0
  Out[..] := 298
  Out[..] := 300
  Out[..] := {298, 300}
*)

Getting information from the TraceHeaders-element:

marmousiTraceHeaders = marmousiData["TraceHeaders"]; 
marmousiData["TraceHeaders"][All]
marmousiData["TraceHeaders", 1, "gx"]
marmousiTraceHeaders[1 ;; 10, {"gx", "gy"}]
marmousiData["TraceHeaders"][{1, 3}]

(*
  Out[..] := {{1, 1, 1, ..}, ..} - returns all data without keys. 
  Out[..] := 0 - only one value
  Out[..] := {{0, 0}, {0, 0}, ..} - mixed span - list of keys and only part of the indexes
  Out[..] := {{"tracl" -> 1, ..}, ..} - only 1 and 3 elements
*)

And the last data element:

marmousiTraces = marmousiData["Traces"]; 
marmousiData["Traces"][6]
marmousiData["Traces", {1, -1}]
marmousiTraces[1 ;; 10 ;; 2]

(*
  Out[..] := {1500., 1500., ..} - first trace 
  Out[..] := {{1500., 1500., ..}, {..}} - first and last traces
  Out[..] := {{1500., 1500., ..}, {..}, ..} - span
*)

Helper Functions

In the package there are several auxiliary functions:

?SEGYDescription
(* SEGYDescription[segyData]. represent the text header of the .segy file in a convenient format *)

SEGYDescription[marmousiData]
(* Out[..] := <3200 white spaces. If in the file exist description - function will show it>*)

Redefined function ArrayPlot:

ArrayPlot[marmousiData]
(* equivalent: 
  ArrayPlot[marmousiData["Traces"]]
  ArrayPlot[Transpose[marmousiData["Traces"][]]]
*)

enter image description here

And you can set a new value to key of the TraceHeaders (in the future is will be working for other elements):

marmousiData["TraceHeaders", 1, "gx"]
SEGYSetValue[marmousiData["TraceHeaders", 1, "gx"], 1]; 
marmousiData["TraceHeaders", 1, "gx"]

(* Out[..] := 0 *)
(* Out[..] := 1 *)

Optimization

To convert IBM 32 Float numbers, the compiled function is used. By default, it is compiled into bytecode, but if you have a C-compiler, you can use it. On average, this gives a 2-fold gain in the speed of data import:

CustomImport[marmousiPath, "SEGY", "ConvertOptions" -> {"CompilationTarget" -> "C"}]

Delayed Loading SEGY Data

This is a very interesting thing that I realized specifically for working with large files (more RAM available). You can not store such a file in the repository.Therefore, the demonstration will be conducted on the same Marmousi model. The function CustomImport has an additional option for calling:

(* by default loading has value "Memory" *)
unloadedMarmousiData = CustomImport[marmousiPath, "SGY", "Loading" -> "Delayed"] 

(*
  Out[..] := SEGYData[{
    SEGYElement["TextHeader", ..], 
    SEGYElement["BinaryHeader", ..], 
    SEGYElement["TraceHeadersUnloaded", {keys..} -> 
      {"File" -> <filepath>, "Convert" -> convertfunc, "Bytes" -> 240, "Positions" -> {poslist}}], 
    SEGYElement["TracesUnloaded", {keys..} -> 
      {"File" -> <filepath>, "Convert" -> convertfunc, "Bytes" -> tracelength, "Positions" -> {poslist}}]
  }]
*)

Unloaded elements contain only references to sections of the indexed file, which allows not storing large amounts of data entirely in memory. You can get the data with a special loader:

CustomDataLoad[unloadedMarmousiData["TracesUnloaded", 1], "SEGY"]
CustomDataLoad[unloadedMarmousiData["TracesUnloaded", {1, 3}], "SEGY"]
CustomDataLoad[unloadedMarmousiData["TracesUnloaded", 1 ;; 3], "SEGY"]

(*  
  Out[..] := {1500., 1500., ..} - only first trace
  Out[..] := {{1500., 1500., ..}, {..}} - traces 1, 3
  Out[..] := {{1500., 1500., ..}, {..}, {..}} - traces 1, 2 and 3
*)

Exporting SEGY Data

This was much more difficult than import, but in this package, data export to SEGY is also available. Only exported data must fully correspond to the internal structure of SEGYData. For example you can create the copy of the MarmousiModel:

CustomExport[StringReplace[marmousiPath, ".segy" -> "Copy.segy"], marmousiData, "SEGY"]

(* Out[..] := <path to copy> *)

P.S. Sorry for my English. I'm ready to answer questions and add features to this article and expand the functionality of the package. Offer your use cases.

POSTED BY: Kirill Belov
Answer
3 months ago

enter image description here - Congratulations! This post is now a Staff Pick as distinguished by a badge on your profile! Thank you, keep it coming!

POSTED BY: Moderation Team
Answer
3 months ago

Thank you for this simple and straightforward step by step loading of SEGY data. it will be very useful.

POSTED BY: Oluseun Sanuade
Answer
3 months ago

Group Abstract Group Abstract