Tim's Brief Notes on making a DEM with ROI_PAC

Disclaimer: This is not intended to be a comprehensive guide, but should enable you to make a DEM using the current ROI_PAC set up at Oxford.

An excellent web resource is David Schmidt's website on ROI_PAC. I always have that open somewhere for reference when making a DEM or deformation interferogram with ROI_PAC.

BEFORE YOU BEGIN:

set up an alias in your ~/.cshrc file to loadroi:

alias loadroi 'source /home/sar1/ROI_PAC_DIST/SAR_CONFIG'

type 'loadroi' before beginning SAR processing with ROI_pac

This should set up the correct environment variables

STAGE 1: Reading in SAR data from CD to '.raw' format

Make a working directory (in /home/satraid1/sardems/<country>/<area>)

e.g. /home/satraid1/sardems/iran/sefidabeh

Make subdirectories with dates of image acquisitions in format yymmdd

Load Data from CD with load_cd.pl or load_cd_lc.pl

Useage: load_cd.pl /cdrom [if data on CD is in /cdrom/*/SCENE1/]
load_cd_lc.pl /cdrom [if data on CD is in /cdrom/*/scene1/]

Do this for both dates... and all frames if DEM is multiple consecutive frames

convert data to ROI_PAC '.raw' format using make_raw.pl

e.g. make_raw.pl ODR SARLEADER1996050406364765T1Of1 960504

[Note: ODR = Delft Precise Orbits; PRC = DPAF Precise Orbits; HDR = Orbits from data on CD.
These should all be available and downloaded... but not necessarily available for every date. HDR should always work!]

Check Doppler

check DOPPLER_RANGE0 values in *.raw.rsc files.
If |DOP| > 0.5 then subtract or add 1
Plot dop.out for both frames (e.g. using 'xmgrace 960505/dop.out 960504/dop.out').

STAGE 2: Make Single Look Complex images (SLCs)

Copy .proc file from /home/satraid1/sardems/resources/date2-date1.proc to working directory

Rename e.g. 960505-960504.proc and edit to change dates and location of DEM and type of orbit used.
Other options should be OK - but can edit if need to. Details

Run process_topo.pl to create slc files:

e.g. process_topo.pl 960505-960504.proc raw slcs

This takes about 1.5 hours for both slcs running on kepler... you can set up the dem file in this time [see stage 5 for how to do this]

Check (and print) 16look slcs using dgx.pl

unlimit stack #only need this first time run dgx in a window
e.g dgx.pl 960505_16rlks.slc #dgx.pl works if a .rsc file exists... otherwise run dgx manually

Check date2_date1_baseline.rsc file (in int_date2_date1 directory)

Is the baseline what you expected ?

If all is OK then can remove some large files:

Remove IMAGERY* files
Remove .raw files [but replace by e.g 'touch 960505.raw' with 0 size file - somewhere later on it checks to see if this file exists!]

STAGE 3: Make interferogram & 1st pass unwrapping attempt

Run process_topo.pl again:

e.g. process_topo.pl 960505-960504.proc slcs unwrapped

This takes a long time again (sefidabeh took ~ 3hours on kepler) and can fall over at several stages. But you can monitor files it creates along the way with dgx.pl

[NB unfortunately a bug means that phase is not displayed on 24 bit screens for *.int files. If you can find an old 8 bit screen then run dgxo.pl and it should work, or convert to rmg (BIL, 4Byte Real) files using cpx2rmg and view in ERMapper] .

Some tips for troubleshooting are on David Schmidt's webpage.

You should get an unwrapped interferogram out at the end of this process, but it is unlikely to be the best possible unwrapped interferogram!

STAGE 4: Improving the Unwrapped interferogram

Two things can generally be improved:

A. the image is likely to have blank patches that did not unwrap with icu. Targetting the unwrapper on these areas can get at them.

B. unwrapping errors: the area has been unwrapped, but icu has made a mistake.

Fixing these is tedious and highly interactive, but necessary.

A. Patches that failed to unwrap

Make new working directory called 'unwrap' of the int_* directory
Restart icu (the unwrapping program) having edited the input file

icu filters the interferogram using the goldstein-werner power spectrum filter, and then unwraps in patches of ~1000 lines moving down the image. If it cannot join neighbouring patches, it just outputs nulls.

When restarting icu, you don't need to filter the image again so can use the filtered interferogram as input to icu, and switch off the filtering flag in the input file.

e.g.
cd unwrap
cp ../icu_unw.in icu_unw1.in
textedit icu_unw1.in & # edit icu input file - see below
ln -s ../filt_ODR_960505-960504.int . #the filtered wrapped interferogram, output by icu
ln -s ../960505-960504.amp . # the amplitude image
icu icu_unw1.in # run icu

#### Lines you MUST change in icu input file ####
Interferogram File (-) = flat_ODR_960505-960504.int
-> Interferogram File (-) = filt_ODR_960505-960504.int
[using pre-filtered interferogram]

Filtering or Unwrapping (-) = Both ![Filtering, Unwrapping, Both]
-> Filtering or Unwrapping (-) = Unwrapping ![Filtering, Unwrapping, Both]
[no need to filter as well]

Filtered Interferogram File          (-) = filt_ODR_960505-960504.int
-> Filtered Interferogram File          (-) = /dev/null
Unwrapped Phase File                 (-) = filt_ODR_960505-960504_c45.phase
Unwrapped Amplitude File             (-) = filt_ODR_960505-960504_c45.amp
Flag File                            (-) = filt_ODR_960505-960504_c45.flg
Phase Sigma Correlation File         (-) = filt_ODR_960505-960504_c45.phsig
-> Unwrapped Phase File                 (-) = filt_ODR_960505-960504_c45_1.phase
-> Unwrapped Amplitude File             (-) = filt_ODR_960505-960504_c45_1.amp
-> Flag File                            (-) = filt_ODR_960505-960504_c45_1.flg
-> Phase Sigma Correlation File         (-) = filt_ODR_960505-960504_c45_1.phsig
[file name changes]

#### Lines you CAN change in icu file #####
Start and End Sample                 (-) = 201 6112
Starting Line                        (-) = 1
Number of Lines                      (-) = 6143          !number of lines to process
[Edit these to specify the region you wish to try to unwrap again - work it out from dgx]

Filter Exponent (-) = 0.55 !PS Filter strength
[This is alpha from the gold-wern filter. Increasing this will increase the % of image unwrapped, but it may also increase the number of unwrapping errors... particularly if there are steep phase gradients.]
#############################

Repeat this process, making lots of files containing sucessfully unwrapped bits of the interferograms.
Connect them together in ERMapper, adding the appropriate number of 2 pi's so the patches match.

copy appropriate template .ers files over from /home/satraid1/sardems/resources/*ers and edit them so they have the right names, lines/columns and number of bands.

Set up multi-layer ERMapper algorithm and use formula window to add multiples of 2 pi to the patches.

If there are no unwrapping errors then save as single phase file.
Also save single .amp and .phsig files [see C below]

Otherwise...

B. Fixing Unwrapping Errors

Identify erroneous areas

look for discontinuities in the unwrapped phase - esp straight lines.
rewrapping interferogram at say 1 fringe = 30 radians or other interval may help

Make ERMapper 'regions' that encompass these areas - you will need to zoom right in and adjust the polygons precisely to get this right.
Use formula window to adjust phase of each of these regions

e.g. "if inregion(r1) then i1+3*2*pi
else if inregion(r2) then i1+1*2*pi
else i1"

Again save the result as a .phase file

-> e.g. filt_ODR_960505-960504_c45_all.phase

C. Get files back into right format for starting up again

icu.pl - the script that usually calls icu - does a number of other things in addition to just unwrapping the phase.

It converts the .phase and .amp r4 files into .unw files containing phase and amplitude (rmg / BIL). It also makes a phase noise map and a mask file.

So you need to do this:

1. Make a single .amp file to match your .phase file. This should have amplitude everywhere you have been able to unwrap, and zeroes everywhere else. Do this in ERMapper using the .amp files from your new patches and the amplitude band of the original .unw file.

-> e.g. filt_ODR_960505-960504_c45_all.amp

2. Make single .phsig file in same way. The original phsig data is in the second band of phase_var_ODR.msk, which you should mv to the unwrap directory. Do not worry that the phase_var_ODR.msk file is 0 where unwrapping has failed - that gets sorted for the rest of the file later.

-> e.g. filt_ODR_960505-960504_c45_all.phsig

3. Run these commands (changing dates as appropriate) from the unwrap directory to move these composite files back to the master directory:

mv filt_ODR_960505-960504_c45_all.phase ../filt_ODR_960505-960504_c45.phase
mv filt_ODR_960505-960504_c45_all.amp ../filt_ODR_960505-960504_c45.amp
mv filt_ODR_960505-960504_c45_all.phsig ../filt_ODR_960505-960504_c45.phsig

4. Run icu2.pl in int* directory to tidy up and create new .unw file. Use syntax from int*/log file.

e.g. icu2.pl flat_ODR_960505-960504 filt_ODR_960505-960504_c45 960505-960504 phase_var_ODR filt_ODR_960505-960504 psfilt 0.55 Phase_Sigma 5 0.45

STAGE 5: Sort out course resolution DEM, make simulation, attempt match

Make sure you have your course resolution DEM ready at this stage. If using GTOPO30 then I use a latlon grid, oversampled to around 100~m resolution (i.e. 10 times as many columns and rows) , and cut out for region surrounding frame. This should be a Signed16BitInteger file with associated .rsc file. The location and file name should agree with that given in the .proc file. Example rsc file: /home/satraid1/sardems/resources/dem.rsc

From the working directory (the one containing the .proc file) restart process_topo.pl

e.g process_topo.pl 960505-960504.proc unwrapped

Unless you have a preexisting DEM that is better than GTOPO30, this will fall over when trying to match the simulated phase from the topography to the real unwrapped phase. If you have a good preexisting DEM and it does the match properly then it should run through to near the end, so skip to stage 7.

If it falls over, it will be during the execution of offset.pl... the error message will be
>PAUSE: p must be at least 2 in moment
>To resume execution, type: go
>Any other input will terminate the program.
Hit return to terminate program

STAGE 6: Manually match Simulation to radar image and restart processing

This involves warping the SIM_2rlks.hgt file into a file in radar geometry... radar_2rlks.hgt

Done by picking tie points in ERMapper and using its warping facility:

1. Copy ERMapper header files for SIM_2rlks.hgt and date2-date1_2rlks.cor to int* directory, and edit as necessary to make them the right size. The .cor file is used because it has already had 2 looks taken and therefore its amplitude band can be compared directly to the simulated amplitudes in SIM_2rlks.hgt

2. Set up and save 2 algorithms in ERMapper containing amplitude images from the SIM and .cor files (sim.alg, cor.alg)

3. Use ERMapper's "Define Ground Control Points" button to pick tie points - pick around 20 from all areas of the image if possible. It's quicker after the first 4. The "from" algorithm is sim.alg and the "to" algorithm is the cor.alg

4. Rectify the dataset in ERMapper. Use linear transformation and bilinear interpolation. Sometimes this only works on fresnel - don't know why. Important: in 'setup', change 'extents' so that they the same as upper left and lower right hand corner of the .cor file. Save file as radar_2rlks.hgt.

5. Copy date2-date1_2rlks.cor.rsc to radar_2rlks.hgt.rsc

Restart topography making... this time using process_topo2.pl, and again from the unwrapped stage. Hopefully this will take us as far as a sch.dem file!

e.g. process_topo2.pl 960505-960504.proc unwrapped

STAGE 7: Post Processing of sch file to DEM

Make post_proc directory off working directory

ln -s ../int*/sch* .

Copy multimosaic.rdf file from resources directory and edit (file names, utm zone)

Change 2nd "Post Spacing" number in sch.hdr file to negative (e.g. -40)

run 'multimosaic multimosaic.rdf' to project into utm zone from sch coordinates

Look at your lovely DEM in ERMapper - copy a header file from resources directory