In [4]:

    

import numpy as np


    

In [5]:

    

with open('../data/OpendTect_seismic_w_header.txt') as f:
    raw = f.readlines()

for l in raw[:3]:
    print(l)


    

    




500	4	251

500	800	4010	2017	-2906	-3810	-574	-67	111	2478	161	-3715	-23	4081	1502	-839	-1140	-1389	1691	2517	-4329	-6591	137	2067	-1882	17	6058	5556	-1640	-4788	-234	353	-6036	-3770	6272	5767	-297	2328	3853	-2052	-3509	-1098	-2676	-2200	2465	3038	-388	-2323	-1199	2789	4538	-150	-4094	-1897	554	683	1461	1556	-107	-549	639	1093	127	-1069	-1343	-1337	-1765	-1351	148	978	562	-718	-1560	-88	1567	95	-1701	-509	1154	1102	327	-440	-230	770	313	-285	1736	2890	-26	-2190	-443	708	-909	-1602	-430	-7	428	2917	4634	1716	-3368	-4049	542	3138	-436	-3346	-303	1879	-1513	-3469	-759	162	-1730	-1244	558	862	2528	5111	3048	-2938	-5942	-2792	3402	5929	1410	-3193	-60	5355	3403	-2902	-4891	-2914	-1999	-1045	2166	3105	-1211	-4331	-1775	1090	378	-991	-1379	-1105	399	2348	3768	3853	941	-2149	-280	2832	1455	-946	-874	-78	1399	2064	-817	-2932	-1238	-12	894	3102	1488	-2868	-2743	-1172	-2787	-3042	-116	1914	1522	-1225	-3759	-2090	529	394	1867	3709	-832	-5309	-982	4119	3107	2439	4313	2441	-2228	-2061	2481	3273	-962	-3057	-1322	-844	-1079	1384	1946	-2617	-4526	-208	1962	-1356	-2986	-660	385	-649	-323	1125	1145	-35	-261	1164	2318	1321	-206	98	810	5	-618	302	760	-28	299	2226	2806	816	-1161	-1350	-1182	-1925	-2312	-1011	909	1198	-579	-1814	-312	1370	-179	-2701	-2129	512	1777	1202	100	-533	-608	-688	-350	817	721

500	801	2860	1944	-1827	-4126	-2270	-425	1154	3166	502	-3359	233	4488	1180	-2849	-2941	-1271	2255	2532	-3520	-4864	938	1161	-3033	367	7021	5164	-2034	-3318	615	-1258	-6938	-2692	6675	4986	-1077	1217	3253	-1006	-2650	-1246	-1993	-1178	2308	2939	-421	-3701	-1882	4322	4827	-2232	-4405	242	1249	-438	1452	2416	-824	-2294	400	1963	-257	-2169	-890	149	-974	-762	1095	1180	-183	-1245	-1616	-74	1750	463	-1501	-662	675	638	212	-721	-1234	8	863	762	2027	2604	-246	-2463	-966	332	-669	-971	19	313	872	2814	3254	178	-3017	-2036	1601	2275	-1096	-2579	498	2048	-1307	-3759	-2059	-720	-1214	-372	940	611	1507	4218	3374	-2172	-5764	-2829	3233	5270	1155	-2245	664	4429	2258	-2397	-3697	-2834	-2522	-984	2634	3476	-878	-3826	-993	1590	-138	-1761	-1330	-1049	-330	1616	3289	3346	1004	-1513	-126	2635	1533	-1189	-1682	-743	1012	1977	-242	-2217	-1011	12	735	2645	1430	-2440	-2499	-760	-1816	-2016	787	2696	1668	-1466	-3716	-2133	-163	-891	895	4788	1423	-6304	-5389	2568	5548	3087	2885	3884	994	-2585	-1253	1528	-29	-2387	-581	1882	1349	-107	-2015	-4942	-4562	794	4113	1111	-2427	-1709	-287	-1047	-1165	802	1992	1007	-45	640	2114	1781	-617	-1566	425	1666	116	-1045	66	911	393	800	2281	1904	-543	-1709	-658	-524	-2269	-2413	377	2029	33	-1867	-495	1074	-495	-2617	-2062	-351	515	1437	2108	637	-1678	-1644	181	574	-712

In [6]:

    

header = raw.pop(0).split()


    

In [7]:

    

len(header)


    

    
Out[7]:





3

In [8]:

    

import os, gzip

def read_odt_seismic(filename, start_time=None, sample_interval=None):
    """
    Read an ascii seismic volume from OdT.
    
    Returns a timebase, arrays representing inline and xline numbers,
    and an array of data.
    
    """
    # Sniff file compression.
    file_extension = os.path.splitext(filename)[1]
    
    if file_extension[-2:] == 'gz':
        with gzip.open(filename, 'rb') as f:
            raw = f.readlines()
    else:
        with open(filename, 'r') as f:
            raw = f.readlines()

    # Sniff for a header.
    if len(raw[0].split()) == 3:
        s, i, n = raw.pop(0).split()
        start_time, sample_interval, number_samples = float(s), float(i), float(n)

    # Gather the data.
    inlines, xlines, traces = [], [], []

    for line in raw:
        l = line.split()
        inlines.append(int(l.pop(0)))
        xlines.append(int(l.pop(0)))
        trace = np.loadtxt(l)
        traces.append(trace)

    # Prepare everything for return.
    number_inlines = max(inlines) - min(inlines) + 1
    number_xlines = max(xlines) - min(xlines) + 1
    number_samples = len(traces[0])  # Nevermind what the header said
    
    if start_time and sample_interval:
        end_time = start_time + (number_samples-1)*sample_interval
        time_basis = np.linspace(start_time, end_time, number_samples)
    else:
        time_basis = None
    
    data = np.array(traces)
    data = np.reshape(data, (number_inlines, number_xlines, number_samples))

    return time_basis, np.array(inlines), np.array(xlines), data


    

In [10]:

    

t, i, x, data = read_odt_seismic("../data/OpendTect_seismic_w_header.txt.gz", 500, 4)


    

In [11]:

    

import matplotlib.pyplot as plt
%matplotlib inline

plt.imshow(data[..., 100], cmap='gray')
plt.show()


    

In [12]:

    

t


    

    
Out[12]:





array([  500.,   504.,   508.,   512.,   516.,   520.,   524.,   528.,
         532.,   536.,   540.,   544.,   548.,   552.,   556.,   560.,
         564.,   568.,   572.,   576.,   580.,   584.,   588.,   592.,
         596.,   600.,   604.,   608.,   612.,   616.,   620.,   624.,
         628.,   632.,   636.,   640.,   644.,   648.,   652.,   656.,
         660.,   664.,   668.,   672.,   676.,   680.,   684.,   688.,
         692.,   696.,   700.,   704.,   708.,   712.,   716.,   720.,
         724.,   728.,   732.,   736.,   740.,   744.,   748.,   752.,
         756.,   760.,   764.,   768.,   772.,   776.,   780.,   784.,
         788.,   792.,   796.,   800.,   804.,   808.,   812.,   816.,
         820.,   824.,   828.,   832.,   836.,   840.,   844.,   848.,
         852.,   856.,   860.,   864.,   868.,   872.,   876.,   880.,
         884.,   888.,   892.,   896.,   900.,   904.,   908.,   912.,
         916.,   920.,   924.,   928.,   932.,   936.,   940.,   944.,
         948.,   952.,   956.,   960.,   964.,   968.,   972.,   976.,
         980.,   984.,   988.,   992.,   996.,  1000.,  1004.,  1008.,
        1012.,  1016.,  1020.,  1024.,  1028.,  1032.,  1036.,  1040.,
        1044.,  1048.,  1052.,  1056.,  1060.,  1064.,  1068.,  1072.,
        1076.,  1080.,  1084.,  1088.,  1092.,  1096.,  1100.,  1104.,
        1108.,  1112.,  1116.,  1120.,  1124.,  1128.,  1132.,  1136.,
        1140.,  1144.,  1148.,  1152.,  1156.,  1160.,  1164.,  1168.,
        1172.,  1176.,  1180.,  1184.,  1188.,  1192.,  1196.,  1200.,
        1204.,  1208.,  1212.,  1216.,  1220.,  1224.,  1228.,  1232.,
        1236.,  1240.,  1244.,  1248.,  1252.,  1256.,  1260.,  1264.,
        1268.,  1272.,  1276.,  1280.,  1284.,  1288.,  1292.,  1296.,
        1300.,  1304.,  1308.,  1312.,  1316.,  1320.,  1324.,  1328.,
        1332.,  1336.,  1340.,  1344.,  1348.,  1352.,  1356.,  1360.,
        1364.,  1368.,  1372.,  1376.,  1380.,  1384.,  1388.,  1392.,
        1396.,  1400.,  1404.,  1408.,  1412.,  1416.,  1420.,  1424.,
        1428.,  1432.,  1436.,  1440.,  1444.,  1448.,  1452.,  1456.,
        1460.,  1464.,  1468.,  1472.,  1476.,  1480.,  1484.,  1488.,
        1492.,  1496.,  1500.])

In [ ]:

    

# import h5py
# h5f = h5py.File('../data/seismic.h5', 'w')
# h5f.create_dataset('data', data=data, compression='lzf')
# h5f.create_dataset('time', data=t, compression='lzf')
# h5f.close()