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:
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:
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:
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:
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:
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:
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!