Modeling Earthquakes

As promised, here is what I learned about how to model earthquakes.

Disclaimer! The process is extremely tedious (or at least that’s what I think).

Before Gareth left on vacation, he showed me how to model an earthquake that occurred in Huarichancara, Peru back in 2016. Here is the end result:

Astrid blog 2 img 1.jpg

I will explain what you are looking at later on.

There are numerous steps to modeling, but I’ll break it down.

In my previous post I discussed the process involved with obtaining interferograms. These interferograms are layered through what is called BIL (Band Interleaved by Line); this is a method that allows for image encoding of data obtained from satellites.

The layering is:    
    Phase
    Amplitude
    Phase
    Amplitude, and so on.

The phase layer is what is used to actually model the earthquake, whereas the amplitude is the topographic aspect used to make the interferogram.

The first step in the modeling process is to separate the phase and amplitude. Before I get too far ahead of myself, I must first make a modeling directory for each of my earthquakes. Then within that, I must make another directory for data and okinv (I’ll explain what that is later) and within my data directory I make another one for ascending or descending (depending on where the earthquake is).

As I mentioned, I need to split my amplitude and phase of the interferogram. The command to do so in the Linux terminal (which is what I use to organize files, make directories, and run programs on the computer) is rmg2mag_phs. This is what allows me to split the unwrapped interferogram into its component parts.

More specifically, what I type into the terminal is this:

rmg2mag_phs filt_topophase.unw.geo filt_topophase.amp.geo filt_topophase.phs.geo 4747

As you can see, I split the interferogram into the phase and amplitude. The number at the end represents the number of columns in the interferogram which is imperative to add at the end of the command.

Gareth showed me an easy way to obtain that column value. Which is this command:

head filt_topophase.unw.geo.vrt

The “head” command shows the first several lines in the file.

After I have made the necessary directories and have now split my phase and amplitude, the next step is to crop and mask the interferogram.

Gareth ran through the steps with me, and recently I have been working on modeling an earthquake that occurred near Hongtu, China. For that earthquake, here is what the unwrapped interferogram looks like:

Astrid blog 2 img 2.jpg

Now that I have the interferogram split, I can now go to MatLab to edit the script crop_data.m where I input the information for the Hongtu, China earthquake. Fortunately for me, Gareth created all the scripts with me, and I can just copy those into the appropriate directories then edit as necessary.

Here is how the crop_data.m script appears on MatLab:

% crop_data.m -- a script to load in various raster data files, view and crop them

% gjf, 03-feb-2017
close all
infile = 'filt_topophase.phs.geo'; % interferogram unwrapped phase file
incols = 5400;% number of columns in that file
outfile = 'asc_crop_1200.r4'; % name for our cropped interferogram file

% open inteferogram file for reading...

infileid = fopen(infile,'r');

% and read it (and transpose it at the same time!)
% (this assumes your file is a float32 [= 4 byte real] file)

indata = fread(infileid,[incols,Inf],'float32')';

% convert to los displacement in meters!

indata=indata*0.0555/(4*pi); % sentinel-1 wavelength is 5.55 cm

% and plot it!
figure
imagesc(indata)
colorbar
axis image

% crop the interferogram down to the part with the earthquake

cropdata=indata(3000:4200,1100:2300);

figure
imagesc(cropdata)
colorbar
axis image

% write out the cropped file

outfileid=fopen(outfile,'w');
fwrite(outfileid,cropdata','float32');% re-transpose it here

% and close your files

fclose(infileid);
fclose(outfileid);

The %’s are put on side notes or on anything I don’t want MatLab to run.

Here is the cropped and uncropped images:

Astrid blog 2 img 3.jpg

 

After this, I need to then mask out anything that is not relevant to the earthquake such as lakes to exclude points of low correlation.

I must do the same for the correlation file and input the correct information in the crop_data_correlation.m script where information on the new registration point (to know where the whole interferogram is referenced to) for my cropped image is inputted.

Here is that script:

% crop_data.m -- a script to load in various raster data files, view and crop them

% gjf, 03-feb-2017
close all
intfile = 'filt_topophase.phs.geo';% interferogram file to be cropped and masked
corfile = 'topophase.cor.geo.r4'; % correlation file (information for mask)
incols = 5400;% number of columns in that file
outfile = 'asc_crop_mask_1200.r4'; % name for our cropped, masked interferogram file
regx = 101.14166666666667 ; % longtidue of registration point
regy = 38.63527777777778 ;% latitude of registration point
dx =0.0002777777777777778 ; % longitude pixel size
dy = 0.0002777777777777778 ;% latitude pixel size

mincor = 0.25; % minimum correlation for mask

startx = 1100;% x coordinate of the start of your crop
starty = 3000;% y coordinate of the same thing
cropsize = 1200;

% open inteferogram file for reading...

intfileid = fopen(intfile,'r');
corfileid = fopen(corfile,'r');

% and read it (and transpose it at the same time!)
% (this assumes your file is a float32 [= 4 byte real] file)

intdata = fread(intfileid,[incols,Inf],'float32')';
cordata = fread(corfileid,[incols,Inf],'float32')';

% convert to los displacement in meters! (only for interferograms!)

intdata=intdata*0.0555/(4*pi); % sentinel-1 wavelength is 5.55 cm

% and plot it!

%imagesc(indata)
%colorbar
%axis image

% crop the interferogram down to the part with the earthquake

cropint=intdata(starty:starty+cropsize-1,startx:startx+cropsize-1); % values substituted with variables
cropcor=cordata(starty:starty+cropsize-1,startx:startx+cropsize-1); % values substituted with variables

% calculate new registration point

shiftx = (startx-1)*dx;% shift in x direction for cropped image to new registration point
shifty = (-1)*(starty-1)*dy;% shift in y direction for same thing

newregx = regx + shiftx; % to find new registration point in x direction
newregy = regy + shifty; % to find new registration point in y direction

% calculate a mask

mask = cropcor>mincor ;
maskint = cropint.*mask ;

figure
imagesc(cropint)
colorbar
axis image

figure
imagesc(maskint)
colorbar
axis image

% write out the cropped file

outfileid=fopen(outfile,'w');
fwrite(outfileid,maskint','float32');% re-transpose it here

% and close your files

fclose(intfileid);
fclose(corfileid);

fclose(outfileid);

percent_masked = (1-(nnz(mask)/cropsize^2))*100;


disp(' ');
disp(sprintf('%6.3f percent of pixels masked at threshold %f',percent_masked,mincor))
disp(sprintf('registration coord.: %f%f',newregx,newregy))
disp(' ');

cropregdata=[regx regy; dx dy; cropsize cropsize];
save cropregdata.dat cropregdata -ascii

Now when I run this script, this is the image that appears:

Astrid blog 2 img 4.jpg

The top image is the unmasked and the bottom is the masked.

At this point, what has been done is the interferogram and correlation files have been cropped and masked and the new registration point has been found in the cropped image.

Next, the LOS (Line of Sight) needs to be be split into the incidence and azimuth which then need to be cropped.

From then on, I need to start the quadtree decomposition. This process is data downsampling for image compression. At the moment, I am having difficulties with this part of the process.

When I figure this out thoroughly (meaning asking Gareth tons of questions when he returns), then I will post more!