-
Notifications
You must be signed in to change notification settings - Fork 64
Pisco earthquake 3: Preparing the InSAR LOS data
The InSAR scene for this event has about 12 million points. It's not feasible to invert this many points so we must subsample the data. I prefer using the QuadTree technique, this is a gradient based technique that will produce more points in areas of high gradient, where the signal is changing rapidly and thus is most important to model correctly. A good summary of InSAR downsampling methods can be found in Lohman and Simmons (2005). Sadly, there is no QuadTree sampler function in Python so we must do this in Matlab by using the function qtdecomp
. I have put a Matlab script titled quadtree_los_pisco.m
in the examples folder. Mudpy uses geographic coordinates so make sure you perform the down sampling on a lat/lon file and not one in UTM coordinates.
Before doing the QuadTree decomposition we must make a file with the coordiantes of a polygon outlining the area covered by the data. Again we can use Matlab functions inpoly
and impoly
. I have created a Matlab script insar_paths.m
in the examples folder to interactively select the regions. You can run it as many times as you need, changing the fout
variable to change the file name for each region, for example for the mainland I set fout='/Users/dmelgar/Slip_inv/Pisco_insar/data/region_path.txt'
and the run the script, a figure will come up where you can click to select the path, here's a screen shot of what that looks like when you've made the path:
If there are more regions (for example the islands) run the script again by changing the output filename, for example set fout='/Users/dmelgar/Slip_inv/Pisco_insar/data/region_path2.txt'
and now click on the path to go around the island, this is what you should end up with:
You can do this for the remaining islands, in the interest of time I'll use only these two regions. Now let's run the QuadTree decomposition using quadtree_los_pisco.m
script. Notice there are three parameters you can set
%QuadTree parameters, these are set by trial and error
mindim=16
maxdim=256
thresh=0.02
mindim
and maxdim
specify the minimum and maximum size (in pixels) that each QuadTree can be. You don;t need to change these by much. The important parameter is thresh
this specifies the varaince change before refinement takes place. If thresh is low, there will be lots of points produced. If thresh is high fewer points will be produced. For example if thresh=0.05
we get the following subsampled scene:
This produces 416 subsampled grids or points. This is a bit small, a good number for MudPy is somewhere between 1000-2000 InSAR points, so we set thresh=0.025
and we get the following scene
This has 1431 points, much better. Notice that in the "flat" areas where the LOS measurements aren't changing much there are fewer points, but int he regions of abrupt changes there are more points. This si why QuadTree is preferable over just regular subsampling. This will have also saved a file with this information in $MUD/Pisco_insar/data/los_qtree.txt
with the subsampled data.
Each InSAR subsampled point is treated like an individual "station" as if it were a GPS measurement of the static offset. The final step in preparing the InSAR data is then to make individual .los
files for every subsampled point and place them in the correct folder. For this we return to python, MudPy has an insartools
module which can be used here simply do:
>>> from mudpy import insartools
>>> home='/Users/dmelgar/Slip_inv/'
>>> project_name='Pisco_insar'
>>> quadtree_file='/Users/dmelgar/Slip_inv/Pisco_insar/data/los_qtree.txt'
>>> gflist_file='/Users/dmelgar/Slip_inv/Pisco_insar/data/station_info/insar_los.gflist'
>>> prefix='PI'
>>> insartools.quadtree2mudpy(home,project_name,quadtree_file,gflist_file,prefix)
This should create the following plot. This is just to verify that you are using the correct file and everything looks ok.
The variable prefix
should be a 2 letter code that will be used for the "station" names of each InSAR point. In this case I have used "PI". Now you should look in $MUD/data/statics
you should have 1413 files labeled PI000.los
, PI0001.los
, PI0002.los
and so on. These are the individual station files which will be used by the inversion. If you open one of these files its contents should be:
# LOS(m),los unit vector (positive towards satellite) n,e,u
7.095550000000000468e-02 2.756000000000000116e-01 9.613000000000000433e-01 4.527999999999999803e-01
The first number is the LOS change, the data, and the next three numbers are the north, east and up components of the LOS unit vector that describes the look direction of the SAR instrument. FInally also notice that this should have generated a .gflist
file in $MUD/data/station_info/insar_los.gflist
. We haven't discussed this file, it's the most important file in MudPy, it's main control file for the inversion, it's first few rows should look something like this:
#station lat lon static disp vel tsun strain Static file displacement file velocity file tsunami file strain file static sigmas(n e u) displacement sigmas(n e u) velocity sigmas(n e u) tsunami sigma strain sigmas(5 components?)
PI0000 -76.448800 -13.839300 0 0 0 0 1 /foo/bar /foo/bar /foo/bar /foo/bar /Users/dmelgar/Slip_inv/Pisco_insar/data/statics/PI0000.los 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
PI0001 -76.448800 -13.834700 0 0 0 0 1 /foo/bar /foo/bar /foo/bar /foo/bar /Users/dmelgar/Slip_inv/Pisco_insar/data/statics/PI0001.los 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
PI0002 -76.448800 -13.830200 0 0 0 0 1 /foo/bar /foo/bar /foo/bar /foo/bar /Users/dmelgar/Slip_inv/Pisco_insar/data/statics/PI0002.los 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
PI0003 -76.445200 -13.843800 0 0 0 0 1 /foo/bar /foo/bar /foo/bar /foo/bar /Users/dmelgar/Slip_inv/Pisco_insar/data/statics/PI0003.los 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
PI0004 -76.445200 -13.839300 0 0 0 0 1 /foo/bar /foo/bar /foo/bar /foo/bar /Users/dmelgar/Slip_inv/Pisco_insar/data/statics/PI0004.los 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
PI0005 -76.441500 -13.843800 0 0 0 0 1 /foo/bar /foo/bar /foo/bar /foo/bar /Users/dmelgar/Slip_inv/Pisco_insar/data/statics/PI0005.los 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
Next we'll look at computing the Green's functions and starting the inversion.
Lohman, R. B., & Simons, M. (2005). Some thoughts on the use of InSAR data to constrain models of surface deformation: Noise structure and data downsampling. Geochemistry, Geophysics, Geosystems, 6(1).