Code covered by the BSD License  

Highlights from
Large Data in MATLAB: A Seismic Data Processing Case Study

image thumbnail

Large Data in MATLAB: A Seismic Data Processing Case Study

by

 

These are the files used in the webinar on Feb. 23, 2011.

setup.m
%% Set Up
%% Data Sources
% Two sources of data are used.  
%
% The fault model is a slice from an SEG/EAGE model which was take from
% <http://utam.gg.utah.edu/Inter.LAB1/CH2.lab/lab.mig.pre/lab.html MATLAB 
% example>.
%
% The salt tooth model is from the BP Benchmark data set from:
% <http://software.seg.org/datasets/2D/2004_BP_Vel_Benchmark>
%
% You will need to download the BP Benchmark files to run the
% saltModelMigrationRTM.m and migrateExample.m files.
%% Create Directories
% There are multiple directories used to store data used/generated from the
% scripts in this example set.  Thise section generates the directories
% referenced in other scripts.
d = {'benchmark';               % BP Benchmark data set location
     'faultModelData';          % SEG/EAGE Model data and simulation data
     'saltToothModelData';      % Simulation data for salt model
     'videos'};                 % location for generated videos
 
for i = 1:length(d)
    mkdir(d{i});
end
%% Download Fault Model Data
% Download and create the velocity model dataset
url = 'http://utam.gg.utah.edu/Inter.LAB1/CH2.lab/lab.mig.pre/';
web(url);
velocityModel = str2num(urlread([url,'/velvector']));
velocityModel = reshape(velocityModel,501,201)';
velocityModel = velocityModel(40:139,250:349);
imagesc(velocityModel)
save('faultModelData\velocityModel.mat','velocityModel');

%% Download Benchmark Data
% Download the files needed from the BP Benchmark data set.
cd benchmark
url = 'http://software.seg.org/datasets/2D/2004_BP_Vel_Benchmark';
web(url)
f = {'2004_Benchmark_READMES.pdf'
    'README_Modification'
    'central_shot_674.gif'
    'eage_abstract.pdf'
    'shots0001_0200.segy.gz'
    'shots0201_0400.segy.gz'
    'shots0401_0600.segy.gz'
    'shots0601_0800.segy.gz'
    'shots0801_1000.segy.gz'
    'shots1001_1200.segy.gz'
    'shots1201_1348.segy.gz'
    'vel_6.25m.gif'
    'vel_z6.25m_x12.5m_exact.segy.gz'};

wb = waitbar(0,'Downloading BP Benchmark Data...');
for i = 1:length(f)
    waitbar(i/length(f)-1,wb,['Downloading...',f{i}]);
    if ismember(i,[1:4,length(f)-1])
        urlwrite([url,'/',f{i}],f{i});
    else
        urlwrite(['ftp://seismic.seg.org/pub/datasets/2D/2004_BP_Vel_Benchmark/',f{i}],f{i});
    end
end
close(wb)
disp('Uncompressing files...')
gunzip(f([5:11,end]));
disp('Done')
cd ..
%% Generate traveltime.dat
% Use 2d ray tracing to generate traveltime array for salt tooth segment
% used in migrateExample.m
addpath benchmark
addpath fileReader
addpath migration

%% Read in velocity model data
% SegYFileReader open the SEG Y file and reads the header information.
% Creates an object used to access trace data.

tic
V = SegYFileReader('vel_z6.25m_x12.5m_exact.segy',true,false);
toc
%% Pull out "tooth" segment for processing
idx = 2200:3400;
tic
V = V(:,idx);
toc
subplot(1,2,1)
imagesc(V), axis equal tight
colormap(seismic)
hshot = plot(1,1,'w*');

[nz,nx] = size(V);
wkrs = 4;
fid = fopen('benchmark\travelTime.dat','w');
travelTime = zeros(nz,nx);
subplot(1,2,2)
for ixs = 1:nx
    travelTime(:,:) = ray2d(V,[1 ixs],dx);
    imagesc(travelTime(:,:))
    title(['Traveltime for shot ',num2str(ixs)])
    set(hshot,'XData',ixs);
    drawnow
    fwrite(fid,travelTime,'double');
end
close(fid);
disp('Done-You can now run the demos! :)')

Contact us