6
$\begingroup$

Do you know if exist a library or a package in order to import and work with SEED Seismographic data in Mathematica?

$\endgroup$
2
  • 1
    $\begingroup$ You will have to write your own importers. I wrote my own SAC reader, which I can post, but rdseed is your best option for working with SEED files (use it to convert to SAC). $\endgroup$
    – rm -rf
    Aug 24, 2014 at 0:38
  • 1
    $\begingroup$ Thanks! Can you share your SAC reader? $\endgroup$
    – moaimx
    Aug 24, 2014 at 0:57

1 Answer 1

10
$\begingroup$

The best tool to work with the SEED format would be IRIS' own rdseed program. It might be possible to use Jrdseed with JLink and Mathematica, but I have not tried this. Typically SEED files are converted to other convenient formats for data processing, of which the SAC (Seismic Analysis Code) format is the most common one.

Implementation:

The following is my implementation of a SAC format reader which I used regularly in my work:

(* Helper functions *)
packed = Developer`ToPackedArray;
fromJulianDay = DatePlus[{#, 1, 1}, #2 - 1] &;
monthOfJulianDay = fromJulianDay[#, #2][[2]]&;

(* ReadSAC *)
ReadSAC[file_String, flag: All | "RawHeader" | "Header" | "ShowHeader" | "Data"] :=
        Module[{str = OpenRead[file, BinaryFormat -> True], h1, h2, h3, header, niceHeader, data, fixTrim, tidyUp},
            fixTrim = (ToCharacterCode[#] /. {Longest[x__?Positive], ___} :> FromCharacterCode[{x}])&;
            tidyUp = # /. {
                (HoldPattern["UNUSED" -> _] | HoldPattern["INTERNAL" -> _]) :> Sequence[],
                Rule[_, -12345 | -12345. | "-12345"] -> Sequence[], (* Undefined *)
                Rule[l_?(StringMatchQ[#, "L"~~__]&), r_] :> l -> (r /. {1 -> True, 0 -> False}),
                Rule[l_, r_String] :> l -> fixTrim@r (* strip trailing chars after null char *)
            }&;

            h1 = Thread[
                {
                    "DELTA"  , "DEPMIN"   , "DEPMAX" , "SCALE"    , "ODELTA"   , "B"        , "E"        , "O"        ,
                    "A"      , "INTERNAL" , "T0"     , "T1"       , "T2"       , "T3"       , "T4"       , "T5"       ,
                    "T6"     , "T7"       , "T8"     , "T9"       , "F"        , "RESP0"    , "RESP1"    , "RESP2"    ,
                    "RESP3"  , "RESP4"    , "RESP5"  , "RESP6"    , "RESP7"    , "RESP8"    , "RESP9"    , "STLA"     ,
                    "STLO"   , "STEL"     , "STDP"   , "EVLA"     , "EVLO"     , "EVEL"     , "EVDP"     , "MAG"      ,
                    "USER0"  , "USER1"    , "USER2"  , "USER3"    , "USER4"    , "USER5"    , "USER6"    , "USER7"    ,
                    "USER8"  , "USER9"    , "DIST"   , "AZ"       , "BAZ"      , "GCARC"    , "INTERNAL" , "INTERNAL" ,
                    "DEPMEN" , "CMPAZ"    , "CMPINC" , "XMINIMUM" , "XMAXIMUM" , "YMINIMUM" , "YMAXIMUM" , "UNUSED"   ,
                    "UNUSED" , "UNUSED"   , "UNUSED" , "UNUSED"   , "UNUSED"   , "UNUSED"
                } -> BinaryReadList[str, "Real32", 70]
            ];

            h2 = Thread[
                {
                    "NZYEAR" , "NZJDAY"   , "NZHOUR" , "NZMIN"  , "NZSEC"  , "NZMSEC" , "NVHDR"  , "NORID"   , "NEVID"   ,
                    "NPTS"   , "INTERNAL" , "NWFID"  , "NXSIZE" , "NYSIZE" , "UNUSED" , "IFTYPE" , "IDEP"    , "IZTYPE"  ,
                    "UNUSED" , "IINST"    , "ISTREG" , "IEVREG" , "IEVTYP" , "IQUAL"  , "ISYNTH" , "IMAGTYP" , "IMAGSRC" ,
                    "UNUSED" , "UNUSED"   , "UNUSED" , "UNUSED" , "UNUSED" , "UNUSED" , "UNUSED" , "UNUSED"  , "LEVEN"   ,
                    "LPSPOL" , "LOVROK"   , "LCALDA" , "UNUSED"
                } -> BinaryReadList[str, "Integer32", 40]
            ];

            h3 = Thread[
                {
                    "KSTNM"  , "KEVNM"  , "KHOLE" , "KO"  , "KA"  , "KT0" , "KT1"    , "KT2"    , "KT3"    , "KT4"    ,
                    "KT5"    , "KT6"    , "KT7"   , "KT8" , "KT9" , "KF"  , "KUSER0" , "KUSER1" , "KUSER2" , "KCMPNM" ,
                    "KNETWK" , "KDATRD" , "KINST"
                } -> Map[
                    Composition[StringTrim, StringJoin],
                    Internal`PartitionRagged[BinaryReadList[str, "Character8", 192], {8, 16} ~ Join ~ ConstantArray[8, 21]]
                ]
            ];

            header = tidyUp@Join[h1, h2, h3];
            niceHeader = {
                        "Type" -> "SAC", "StationCode" -> "KSTNM", "Network" -> "KNETWK", "Component" -> "KCMPNM",
                        "Location" -> {"STLA", "STLO"}, "SamplingFrequency" -> 1 / "DELTA", "Samples" -> "NPTS",
                        "Date" -> Hold[fromJulianDay["NZYEAR", "NZJDAY"] ~ Join ~ {"NZHOUR", "NZMIN", "NZSEC"}],
                        "Year" -> "NZYEAR", "Month" -> Hold[monthOfJulianDay["NZYEAR", "NZJDAY"]]
                    } /. header // ReleaseHold // Association;
            If[!StringMatchQ[ToString@flag, "*Header"], data = packed@N@BinaryReadList[str, "Real32"]];

            Close@str;

            Switch[flag,
                "Data", data,
                "Header", niceHeader,
                "RawHeader", header,
                "ShowHeader", header /. Rule -> List // TableForm // Panel,
                All, <|niceHeader, "FileName" -> file, "Data" -> data|>
            ]
        ]

ReadSAC[file_String] := ReadSAC[file, All]

Usage:

You can read the header info for a SAC file as ReadSAC["filename", "Header"] (I'm not including a SAC file for this example):

This is obviously a "cleaned" up version of the header and if you'd like to view the raw header (as described in the SAC format), you can use "RawHeader" as the second argument.

The entire file can be read with ReadSAC["filename"] or ReadSAC["filename", All] and the output is again an association. In addition to the keys found in the header above, the association also has "FileName" and "Data" as keys when reading the entire file.

You can also just ask for the seismic data directly with "Data" as the second argument, and in this case, the output is a list of data points (that have not been corrected for instrument response).

$\endgroup$

Your Answer

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

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