Master’s thesis
Physical Geography and Quaternary Geology, 45 Credits
Department of Physical Geography
Differential Interferometry
and Multiple-Aperture
Interferometry for Retrieving
Three-Dimensional Measurements
of Glacial Surface Velocity
Luke Webber
Preface
This Master’s thesis is Luke Webber’s degree project in Physical Geography and Quaternary
Geology at the Department of Physical Geography, Stockholm University. The Master’s
thesis comprises 45 credits (one and a half term of full-time studies).
Supervisor has been Ian Brown at the Department of Physical Geography, Stockholm
University. Examiner has been Johan Kleman at the Department of Physical Geography,
Stockholm University.
The author is responsible for the contents of this thesis.
Stockholm, 15 September 2016
𝚫𝑩
⊥𝚫𝒇
𝑫𝑪𝒙̅
𝜸𝝈
𝜸𝑪
𝒗𝚫𝑩
⊥𝒙̅
𝑫𝑰𝑺𝑷𝑳𝑨𝑪𝑬𝑴𝑬𝑵𝑻𝝈
𝑫𝑰𝑺𝑷𝑳𝑨𝑪𝑬𝑴𝑬𝑵𝑻𝒙̅
𝑽𝑬𝑳𝑶𝑪𝑰𝑻𝒀𝑥̅
𝐷𝐼𝑆𝑃𝐿𝐴𝐶𝐸𝑀𝐸𝑁𝑇𝜎
𝐷𝐼𝑆𝑃𝐿𝐴𝐶𝐸𝑀𝐸𝑁𝑇𝚫𝑩
⊥𝚫𝒇
𝑫𝑪𝒙̅
𝜸𝝈
𝜸𝑪
𝒗Δ𝐵
⊥Δ𝑓
𝐷𝐶𝐶
𝑣𝛾
𝑡ℎ𝑒𝑟𝑚𝑎𝑙𝑥̅
𝑦import argparse
import time
import sys
import numpy as np
import scipy.fftpack as fft
import spectral.io.envi as envi
from tqdm import tqdm
def compute_interferogram(mst, slv): """
Computes along-track interferogram by splitting doppler frequencies into positive and negative frequencies, corresponding to the forwards and backwards look angles, and computing the double difference interferogram.
:param mst: complex master SLC image
:param slv: complex slave SLC image co-registered to master :return: along-track interferogram
"""
# Transform data along azimuthal (columnar) direction into frequency domain
master_doppler = fft.fft(mst, axis=0) slave_doppler = fft.fft(slv, axis=0)
frequencies = fft.fftfreq(master_doppler.shape[0]) master_doppler_fpos = master_doppler.copy()
master_doppler_fneg = master_doppler.copy() slave_doppler_fpos = slave_doppler.copy() slave_doppler_fneg = slave_doppler.copy() # Taper to zero unwanted frequencies
master_doppler_fpos[np.where(frequencies < 0)] = 0.0 + 0.0j
master_doppler_fneg[np.where(frequencies > 0)] = 0.0 + 0.0j
slave_doppler_fpos[np.where(frequencies < 0)] = 0.0 + 0.0j
slave_doppler_fneg[np.where(frequencies > 0)] = 0.0 + 0.0j
# Inverse transform back into time domain
master_forward_look = fft.ifft(master_doppler_fpos, axis=0) master_backward_look = fft.ifft(master_doppler_fneg, axis=0) slave_forward_look = fft.ifft(slave_doppler_fpos, axis=0) slave_backward_look = fft.ifft(slave_doppler_fneg, axis=0) # Compute double difference interferogram
i1 = slave_backward_look * np.conj(master_backward_look) i2 = slave_forward_look * np.conj(master_forward_look) ifg = i1 * np.conj(i2)
# Normalise signal
ifg /= np.absolute(ifg) return ifg
def main(i_master, q_master, i_slave, q_slave, i_out, q_out, bursts=None): """
Reads input file into memory, computes along-track interferogram, interpolates the interferogram and writes the results to the output file.
:param q_slave: Path to slave file containing real component.
:param i_out: Path to save completed along-track imaginary component. :param q_out: Path to save completed along-track real component. :param bursts: (Optional) Number of bursts present if TOPSAR product. """
print("[{0}] Opening files for reading...".format(time.strftime("%H:%M:%S"))) i_mst = envi.open(i_master)
q_mst = envi.open(q_master) i_slv = envi.open(i_slave) q_slv = envi.open(q_slave)
print("[{0}] Files opened".format(time.strftime("%H:%M:%S")))
print("[{0}] Copying data into arrays...".format(time.strftime("%H:%M:%S"))) master = np.array(i_mst.load()) + 1j * np.array(q_mst.load())
slave = np.array(i_slv.load()) + 1j * np.array(q_slv.load()) ifg = np.empty(master.shape, dtype=complex)
print("[{0}] Data copied into arrays".format(time.strftime("%H:%M:%S"))) if bursts:
burst_height = master.shape[0]//bursts
print("[{0}] Processing data per burst...".format(time.strftime("%H:%M:%S"))) for i in tqdm(range(0, burst_height*bursts, burst_height)):
j = i + burst_height
ifg[i:j, :] = compute_interferogram(master[i:j, :], slave[i:j, :]) else:
print("[{0}] Processing data...".format(time.strftime("%H:%M:%S"))) for i in tqdm(range(0, master.shape[1])):
ifg[:, i] = compute_interferogram(master[:, i], slave[:, i]) print("[{0}] Writing data to disk...".format(time.strftime("%H:%M:%S")))
# In order to preserve backwards compatibility with the BEAM-DIMAP format used by the
# SNAP S1TBX program output files are written with the Spectral Python module
# in big-endian ENVI format.
envi.save_image(i_out, np.real(ifg), dtype='float32', interleave='bsq', byteorder='1', metadata=i_mst.metadata, force=True)
envi.save_image(q_out, np.imag(ifg), dtype='float32', interleave='bsq', byteorder='1', metadata=q_mst.metadata, force=True)
print("[{0}] Data written to disk".format(time.strftime("%H:%M:%S"))) print("[{0}] Finished!".format(time.strftime("%H:%M:%S")))
if __name__ == "__main__":
class Parser(argparse.ArgumentParser): def error(self, message):
sys.stderr.write('error: %s\n' % message) self.print_help()
sys.exit(2)
parser = Parser(description="Along-track interferogram formation")
parser.add_argument("i_master", type=open, help="Path to master real component") parser.add_argument("q_master", type=open, help="Path to master imaginary component") parser.add_argument("i_slave", type=open, help="Path to slave real component")
parser.add_argument("q_slave", type=open, help="Path to slave imaginary component") parser.add_argument("i_out", help="Path to save output real component")
parser.add_argument("q_out", help="Path to save output imaginary component") parser.add_argument("-b", "--bursts", type=int,
help="Number of bursts if input file is TOPSAR product") np.seterr(divide='ignore', invalid='ignore')
args = parser.parse_args()
import sys import gdal import math import argparse import numpy as np import multiprocessing as mp
from scipy.ndimage.filters import generic_filter
from osgeo.gdalconst import GDT_Float32
def window_stdev(arr, size): """
Calculate standard deviation over an 2-dimensional array using a moving window. To avoid edge effects in the final result, ensure displacement observations have an outside buffer of at least 'size'/2 or greater.
:param arr: target array
:param size: size of moving window
:return: array of stardard deviations calculated pixel-by-pixel """
return generic_filter(arr, np.std, size, mode='reflect')
def rot_matrix(a_inc, d_inc, a_az, d_az): """
Generate rotation matrix.
:param a_inc: ascending geometry radar incidence angle in radians :param d_inc: descending geometry radar incidence angle in radians :param a_az: ascending geometry orbital azimuth angle in radians :param d_az: descending geometry orbital azimuth angle in radians :return: rotation matrix
"""
p = 3*math.pi/2
return np.matrix([[math.cos(a_inc), -math.sin(a_inc)*math.sin(a_az - p), -math.sin(a_inc)*math.cos(a_az - p)],
[math.cos(d_inc), -math.sin(d_inc)*math.sin(d_az - p), -math.sin(d_inc)*math.cos(d_az - p)],
[0, -math.cos(a_az - p), -math.sin(a_az - p)], [0, -math.cos(d_az - p), -math.sin(d_az - p)]])
def chunk_array(arr_size, splits): """
Determine the location of array chunks in preparation for multiprocessing. :param arr_size: size of array along axis to be chunked
:param splits: number of chunks the array is to be split into :return: location of chunk splits
"""
axis_pos = np.arange(arr_size)
return [(min(axis_pos[i*arr_size//splits:(i+1)*arr_size//splits]),
max(axis_pos[i*arr_size//splits:(i+1)*arr_size//splits])) for i in range(splits)]
def worker(worker_id, queue): """
Subprocess worker for calculating weighted least squares adjustment. :param worker_id: identifying number
:param queue: queue for passing back data :return: displacement vectors
"""
# Perform weighted least squares adjustment to derive displacement vector from LOS and
# azimuthal observations, values along 3rd axis of final displacement map represent
# up, east and north respectively.
for j in range(chunk_loc[worker_id][0], chunk_loc[worker_id][1]+1): displacement_vectors = np.array([d_a_los[i, j], d_d_los[i, j], d_a_azi[i, j], d_d_azi[i, j]]) if displacement_vectors.any():
rotation_matrix = rot_matrix(asc_incidence[i, j], dsc_incidence[i, j], asc_orbit_azi[i, j], dsc_orbit_azi[i, j]) covariance_matrix = np.asmatrix(np.diag([stdev_d_a_los[i, j],
stdev_d_d_los[i, j], stdev_d_a_azi[i, j], stdev_d_d_azi[i, j]])) try:
a = rotation_matrix.T @ covariance_matrix.I @ rotation_matrix
b = a.I @ rotation_matrix.T @ covariance_matrix.I @ displacement_vectors
displacement_map[i, j] = np.asarray(b)[0] except np.linalg.linalg.LinAlgError:
# Catch singular matrix error, singular matrices are invertible, this
# typically occurs because all values of standard deviation are zero
# within the covariance matrix, the cause of this is due to areas of
# no-data having zero variance, as there is no data, displacements
# along all three dimensions are set to zero.
displacement_map[i, j] = [0.0, 0.0, 0.0] else:
# All measured values of displacement in the LOS and along-direction were zero,
# hence displacements along all three dimensions are set to zero.
displacement_map[i, j] = [0.0, 0.0, 0.0] queue.put([worker_id, displacement_map])
def write_results(displacements, src): """
Extracts displacment components and writes results into tiff files: :param displacements: Numpy array holding results
:param src: GDAL file handle, for copying projection and location data to output files. "up_down.tiff" Up-Down component
"east_west.tiff" East-West component "north_south.tiff" North-South component """
driver = src.GetDriver()
ud = driver.Create("up_down.tiff",
displacements.shape[1], displacements.shape[0], 1, GDT_Float32) ud.SetMetadata(src.GetMetadata())
ud.SetGeoTransform(src.GetGeoTransform()) ud.SetProjection(src.GetProjection())
ud.GetRasterBand(1).WriteArray(displacements[:, :, 0], 0, 0) ud.GetRasterBand(1).FlushCache()
ew = driver.Create("east_west.tiff",
displacements.shape[1], displacements.shape[0], 1, GDT_Float32) ew.SetMetadata(src.GetMetadata())
ew.SetGeoTransform(src.GetGeoTransform()) ew.SetProjection(src.GetProjection())
ew.GetRasterBand(1).WriteArray(displacements[:, :, 1], 0, 0) ew.GetRasterBand(1).FlushCache()
ns = driver.Create("north_south.tiff",
displacements.shape[1], displacements.shape[0], 1, GDT_Float32) ns.SetMetadata(src.GetMetadata())
ns.SetGeoTransform(src.GetGeoTransform()) ns.SetProjection(src.GetProjection())
ns.GetRasterBand(1).WriteArray(displacements[:, :, 2], 0, 0) ns.GetRasterBand(1).FlushCache()
if __name__ == "__main__":
class Parser(argparse.ArgumentParser): def error(self, message):
self.print_help() sys.exit(2)
parser = Parser(description="Three-dimensional displacement field generation") parser.add_argument("asc_los", type=open,
help="Path to LOS measurements from the ascending orbital pass") parser.add_argument("dsc_los", type=open,
help="Path to LOS measurements from the descending orbital pass") parser.add_argument("asc_at", type=open,
help="Path to along-track measurements from the ascending orbital pass") parser.add_argument("dsc_at", type=open,
help="Path to along-track measurements from the descending orbital pass") parser.add_argument("asc_inc", help="Incidence angle of ascending pass")
parser.add_argument("dsc_inc", help="Incidence angle of descending pass") parser.add_argument("asc_azi", help="Orbital azimuth angle of ascending pass") parser.add_argument("dsc_azi", help="Orbital azimuth angle of descending pass") args = parser.parse_args()
# Open files for reading
d_a_los_source = gdal.Open(args.asc_los.name) d_d_los_source = gdal.Open(args.dsc_los.name) d_a_azi_source = gdal.Open(args.asc_at.name) d_d_azi_source = gdal.Open(args.dsc_at.name) # Extract data from raster bands
d_a_los = d_a_los_source.GetRasterBand(1).ReadAsArray() d_d_los = d_d_los_source.GetRasterBand(1).ReadAsArray() d_a_azi = d_a_azi_source.GetRasterBand(1).ReadAsArray() d_d_azi = d_d_azi_source.GetRasterBand(1).ReadAsArray() array_size = d_a_los.shape
# Calculate standard deviation arrays over input displacement vector files.
stdev_d_a_los = window_stdev(d_a_los, 5) stdev_d_d_los = window_stdev(d_d_los, 5) stdev_d_a_azi = window_stdev(d_a_azi, 5) stdev_d_d_azi = window_stdev(d_d_azi, 5)
# Create matrix of ascending and descending incidence/orbital azimuthal angles.
# Ideally these should be calculated on a per-pixel basis and loaded from a raster,
# rather than using an average value over the entire radar scene!
asc_incidence = np.full(array_size, math.radians(float(args.asc_inc))) dsc_incidence = np.full(array_size, math.radians(float(args.dsc_inc))) asc_orbit_azi = np.full(array_size, math.radians(float(args.asc_azi))) dsc_orbit_azi = np.full(array_size, math.radians(float(args.dsc_azi)))
# The weighted least squares calculations are independent of each other, and are calculated
# on a pixel-by-pixel basis, as such the calculations can be parallelised to run on
# multiple cores, for maximum speed! Due to the GIL lock, threading in Python offers no
# speedup, instead new processes are spawned.
displacement_map = np.empty((array_size + (3,))) # Empty array to hold final results
cores = mp.cpu_count() # Number of cores available
# Determine location of chunks for equally splitting the work between cores
chunk_loc = chunk_array(array_size[1], cores) # Queue for worker processes to return data by
result_queue = mp.Queue()
jobs = [mp.Process(target=worker, args=(i, result_queue)) for i in range(cores)] for job in jobs:
job.start()
# There is a limit on how much data can be piped between processes before Process.join()
# will hang indefinitely, instead the calling processes is blocked with a loop until
# all results have been de-queued into a temporary buffer
temp_results = []
temp_results.append(result_queue.get())
# Rebuild final displacement map from temporary buffer
for result in temp_results:
displacement_map[:, chunk_loc[result[0]][0]:chunk_loc[result[0]][1]+1] = \ result[1][:, chunk_loc[result[0]][0]:chunk_loc[result[0]][1]+1]
# Reverse displacements to true polarity, otherwise results will be backwards!
displacement_map *= -1.0
# Write results to file
