0
# Data Access
1

2
Multiple specialized interfaces for accessing seismic data organized by traces, inlines, crosslines, depth slices, and gathers. Each access mode provides a different view of the same underlying data optimized for specific workflows.
3

4
## Capabilities
5

6
### Individual Trace Access
7

8
Direct access to individual seismic traces by trace index. Most fundamental access pattern.
9

10
```python { .api }
11
# Trace access mode (available as f.trace)
12
f.trace[index]           # Read trace by index (returns numpy.ndarray)
13
f.trace[start:end]       # Read multiple traces (returns numpy.ndarray) 
14
f.trace[start:end:step]  # Read traces with step (returns numpy.ndarray)
15
f.trace[index] = data    # Write trace data (if file writable)
16
```
17

18
**Usage Example:**
19

20
```python
21
with segyio.open('data.sgy') as f:
22
    # Read single trace
23
    first_trace = f.trace[0]
24
    print(f"Trace shape: {first_trace.shape}")
25
    print(f"Trace dtype: {first_trace.dtype}")
26
    
27
    # Read multiple traces
28
    traces_10_to_19 = f.trace[10:20]
29
    print(f"Multiple traces shape: {traces_10_to_19.shape}")
30
    
31
    # Read every 10th trace
32
    sampled_traces = f.trace[::10]
33
    
34
    # Iterate over all traces
35
    for i, trace in enumerate(f.trace):
36
        # Process each trace
37
        amplitude = trace.max()
38
        print(f"Trace {i} max amplitude: {amplitude}")
39
        
40
        if i >= 100:  # Limit for example
41
            break
42

43
# Writing traces
44
with segyio.open('data.sgy', 'r+') as f:
45
    # Modify a trace
46
    modified_trace = f.trace[0] * 2.0  # Double amplitudes
47
    f.trace[0] = modified_trace
48
    
49
    # Write synthetic data
50
    synthetic = np.sin(np.linspace(0, 10*np.pi, len(f.samples)))
51
    f.trace[0] = synthetic
52
```
53

54
### Inline Access
55

56
Access seismic data organized by inline number (3D structured files only). Each inline represents a 2D vertical slice parallel to the acquisition direction.
57

58
```python { .api }
59
# Inline access mode (available as f.iline)
60
f.iline[inline_number]              # Read inline slice (returns 2D numpy.ndarray)
61
f.iline[inline_number, crosslines]  # Read subset of crosslines
62
f.iline[inline_number] = data       # Write inline data (if file writable)
63
```
64

65
**Usage Example:**
66

67
```python
68
with segyio.open('3d_volume.sgy') as f:
69
    print(f"Available inlines: {f.ilines}")
70
    
71
    # Read complete inline
72
    inline_100 = f.iline[100]
73
    print(f"Inline shape: {inline_100.shape}")  # (n_crosslines, n_samples)
74
    
75
    # Read subset of crosslines for an inline
76
    partial_inline = f.iline[100, 200:250]  # Crosslines 200-249
77
    
78
    # Iterate over all inlines
79
    for il in f.ilines:
80
        inline_data = f.iline[il]
81
        avg_amplitude = inline_data.mean()
82
        print(f"Inline {il} average amplitude: {avg_amplitude}")
83
    
84
    # Process inline with numpy operations
85
    inline_data = f.iline[150]
86
    
87
    # Apply gain correction
88
    gained = inline_data * 1.5
89
    
90
    # Apply AGC (Automatic Gain Control) 
91
    window_size = 50
92
    agc_data = np.zeros_like(inline_data)
93
    for i, trace in enumerate(inline_data.T):
94
        for j in range(len(trace)):
95
            start = max(0, j - window_size//2)
96
            end = min(len(trace), j + window_size//2)
97
            window_rms = np.sqrt(np.mean(trace[start:end]**2))
98
            if window_rms > 0:
99
                agc_data[i, j] = trace[j] / window_rms
100
```
101

102
### Crossline Access
103

104
Access seismic data organized by crossline number (3D structured files only). Each crossline represents a 2D vertical slice perpendicular to the acquisition direction.
105

106
```python { .api }
107
# Crossline access mode (available as f.xline)
108
f.xline[crossline_number]            # Read crossline slice (returns 2D numpy.ndarray)
109
f.xline[crossline_number, inlines]   # Read subset of inlines
110
f.xline[crossline_number] = data     # Write crossline data (if file writable)
111
```
112

113
**Usage Example:**
114

115
```python
116
with segyio.open('3d_volume.sgy') as f:
117
    print(f"Available crosslines: {f.xlines}")
118
    
119
    # Read complete crossline
120
    crossline_300 = f.xline[300]
121
    print(f"Crossline shape: {crossline_300.shape}")  # (n_inlines, n_samples)
122
    
123
    # Compare inline vs crossline
124
    il_data = f.iline[100]
125
    xl_data = f.xline[300]
126
    
127
    print(f"Inline 100 dimensions: {il_data.shape}")
128
    print(f"Crossline 300 dimensions: {xl_data.shape}")
129
    
130
    # Find intersection point
131
    intersection = f.trace[f.index(inline=100, crossline=300)]
132
    
133
    # Extract same trace from both slices
134
    # Find position of crossline 300 in inline 100
135
    xl_pos = list(f.xlines).index(300)
136
    trace_from_inline = il_data[xl_pos, :]
137
    
138
    # Find position of inline 100 in crossline 300  
139
    il_pos = list(f.ilines).index(100)
140
    trace_from_crossline = xl_data[il_pos, :]
141
    
142
    # These should be identical
143
    assert np.allclose(trace_from_inline, trace_from_crossline)
144
```
145

146
### Fast and Slow Dimension Access
147

148
Access data along the fast and slow dimensions based on file sorting. Provides consistent access regardless of whether file is inline-sorted or crossline-sorted.
149

150
```python { .api }
151
# Fast/slow dimension access (available as f.fast and f.slow)
152
f.fast[index]       # Read along fast dimension
153
f.slow[index]       # Read along slow dimension
154
```
155

156
**Usage Example:**
157

158
```python
159
with segyio.open('data.sgy') as f:
160
    print(f"File sorting: {f.sorting}")
161
    
162
    if f.sorting == segyio.TraceSortingFormat.INLINE_SORTING:
163
        print("Fast=inline, Slow=crossline")
164
        assert f.fast == f.iline
165
        assert f.slow == f.xline
166
    elif f.sorting == segyio.TraceSortingFormat.CROSSLINE_SORTING:
167
        print("Fast=crossline, Slow=inline")
168
        assert f.fast == f.xline
169
        assert f.slow == f.iline
170
    
171
    # Access first line along fast dimension
172
    fast_line = f.fast[f.fast.indices[0]]
173
    
174
    # Access first line along slow dimension  
175
    slow_line = f.slow[f.slow.indices[0]]
176
```
177

178
### Depth Slice Access
179

180
Access horizontal slices at specific sample indices (time/depth slices). Useful for horizon mapping and amplitude analysis.
181

182
```python { .api }
183
# Depth slice access mode (available as f.depth_slice)
184
f.depth_slice[sample_index]          # Read horizontal slice at sample index
185
f.depth_slice[sample_index] = data   # Write depth slice (if file writable)
186
```
187

188
**Usage Example:**
189

190
```python
191
with segyio.open('3d_volume.sgy') as f:
192
    print(f"Sample range: {f.samples}")
193
    
194
    # Read depth slice at 500ms (sample index 125 for 4ms sampling)
195
    sample_index = 125
196
    depth_slice = f.depth_slice[sample_index]
197
    print(f"Depth slice shape: {depth_slice.shape}")  # (n_inlines, n_crosslines)
198
    
199
    # Create amplitude map for interpretation
200
    import matplotlib.pyplot as plt
201
    
202
    plt.figure(figsize=(10, 8))
203
    plt.imshow(depth_slice, cmap='seismic', aspect='auto')
204
    plt.title(f'Amplitude at {f.samples[sample_index]}ms')
205
    plt.xlabel('Crossline')
206
    plt.ylabel('Inline') 
207
    plt.colorbar(label='Amplitude')
208
    
209
    # Find maximum amplitude location
210
    max_pos = np.unravel_index(np.argmax(np.abs(depth_slice)), depth_slice.shape)
211
    max_inline = f.ilines[max_pos[0]]
212
    max_crossline = f.xlines[max_pos[1]]
213
    print(f"Maximum amplitude at IL={max_inline}, XL={max_crossline}")
214
    
215
    # Extract amplitude attributes along a horizon
216
    horizon_samples = np.full(depth_slice.shape, sample_index)  # Flat horizon
217
    
218
    # For variable horizon, load from interpretation
219
    # horizon_samples = load_horizon_picks('horizon.dat')
220
    
221
    # Extract attributes along horizon
222
    horizon_amplitudes = np.zeros(depth_slice.shape)
223
    for i, il in enumerate(f.ilines):
224
        for j, xl in enumerate(f.xlines):
225
            sample_idx = int(horizon_samples[i, j])
226
            if 0 <= sample_idx < len(f.samples):
227
                horizon_amplitudes[i, j] = f.depth_slice[sample_idx][i, j]
228
```
229

230
### Gather Access
231

232
Access pre-stack gathers organized by Common Depth Point (CDP) or other grouping criteria. Essential for velocity analysis and pre-stack processing.
233

234
```python { .api }
235
# Gather access mode (available as f.gather)
236
f.gather[cdp_number]              # Read gather by CDP number
237
f.gather[cdp_number, offsets]     # Read subset of offsets
238
```
239

240
**Usage Example:**
241

242
```python
243
# Pre-stack data with multiple offsets
244
with segyio.open('prestack.sgy') as f:
245
    print(f"Available offsets: {f.offsets}")
246
    
247
    # Get first CDP number
248
    first_cdp = f.header[0][segyio.TraceField.CDP]
249
    
250
    # Read complete gather  
251
    gather = f.gather[first_cdp]
252
    print(f"Gather shape: {gather.shape}")  # (n_offsets, n_samples)
253
    
254
    # Read subset of offsets
255
    near_offsets = f.gather[first_cdp, f.offsets[:3]]  # First 3 offsets
256
    
257
    # Velocity analysis - compute semblance
258
    def compute_semblance(gather, times, velocities):
259
        """Compute semblance for velocity analysis."""
260
        semblance = np.zeros((len(times), len(velocities)))
261
        
262
        for i, t in enumerate(times):
263
            for j, v in enumerate(velocities):
264
                # NMO correction
265
                corrected = np.zeros_like(gather)
266
                for k, offset in enumerate(f.offsets):
267
                    t_nmo = np.sqrt(t**2 + (offset/v)**2)
268
                    # Interpolate to corrected time
269
                    # (simplified - would use proper interpolation)
270
                    sample_idx = int(t_nmo * 1000 / 4)  # 4ms sampling
271
                    if sample_idx < gather.shape[1]:
272
                        corrected[k, i] = gather[k, sample_idx]
273
                
274
                # Compute semblance
275
                stack = corrected.sum(axis=0)
276
                energy_stack = (stack**2).sum()
277
                energy_individual = (corrected**2).sum()
278
                
279
                if energy_individual > 0:
280
                    semblance[i, j] = energy_stack / energy_individual
281
        
282
        return semblance
283
    
284
    # Example velocity analysis
285
    times = np.arange(0, 2.0, 0.004)  # 0-2 seconds
286
    velocities = np.arange(1500, 4000, 50)  # 1500-4000 m/s
287
    
288
    semblance = compute_semblance(gather, times, velocities)
289
    
290
    # Find maximum velocity at each time
291
    optimal_velocities = velocities[np.argmax(semblance, axis=1)]
292
```
293

294
### Multi-dimensional Indexing
295

296
Combined indexing across multiple dimensions for precise data access.
297

298
```python
299
with segyio.open('3d_volume.sgy') as f:
300
    # Get trace index from inline/crossline coordinates
301
    trace_idx = f.index(inline=100, crossline=200)
302
    trace = f.trace[trace_idx]
303
    
304
    # Alternative: direct access via coordinates
305
    trace = f.trace[f.index(100, 200)]
306
    
307
    # Bulk operations on subsets
308
    inline_subset = f.iline[100:110]  # Multiple inlines
309
    crossline_subset = f.xline[200:210]  # Multiple crosslines
310
    
311
    # Advanced slicing
312
    subvolume = f.iline[100:110, 200:210]  # Sub-cube extraction
313
```
314

315
## Data Processing Examples
316

317
### Amplitude Statistics
318

319
```python
320
def compute_amplitude_stats(filename):
321
    """Compute amplitude statistics for SEG-Y file."""
322
    with segyio.open(filename) as f:
323
        all_amplitudes = []
324
        
325
        # Sample every 100th trace for efficiency
326
        for i in range(0, f.tracecount, 100):
327
            trace = f.trace[i]
328
            all_amplitudes.extend(trace)
329
        
330
        amplitudes = np.array(all_amplitudes)
331
        
332
        return {
333
            'min': amplitudes.min(),
334
            'max': amplitudes.max(), 
335
            'mean': amplitudes.mean(),
336
            'std': amplitudes.std(),
337
            'rms': np.sqrt(np.mean(amplitudes**2))
338
        }
339
```
340

341
### Data QC and Validation
342

343
```python
344
def qc_seismic_data(filename):
345
    """Perform basic QC on seismic data."""
346
    issues = []
347
    
348
    with segyio.open(filename) as f:
349
        # Check for dead traces
350
        dead_traces = []
351
        for i in range(min(1000, f.tracecount)):
352
            trace = f.trace[i]
353
            if np.all(trace == 0) or np.all(np.isnan(trace)):
354
                dead_traces.append(i)
355
        
356
        if dead_traces:
357
            issues.append(f"Dead traces found: {len(dead_traces)}")
358
        
359
        # Check for clipped data
360
        clipped_traces = []
361
        for i in range(min(1000, f.tracecount)):
362
            trace = f.trace[i]
363
            if np.any(np.abs(trace) >= 32767):  # Common clipping level
364
                clipped_traces.append(i)
365
        
366
        if clipped_traces:
367
            issues.append(f"Clipped traces found: {len(clipped_traces)}")
368
        
369
        # Check geometry consistency
370
        if not f.unstructured:
371
            expected_traces = len(f.ilines) * len(f.xlines) * len(f.offsets)
372
            if f.tracecount != expected_traces:
373
                issues.append(f"Trace count mismatch: {f.tracecount} vs {expected_traces}")
374
    
375
    return issues
376
```