Revision f24076bbf4a4a4af8e90100f344e10fb237e9b3b authored by Zhang Yunjun on 08 November 2022, 07:17:03 UTC, committed by GitHub on 08 November 2022, 07:17:03 UTC
+ .circleci/config.yml update: - use $BASH_ENV to set and share environment variable (PATH) among multiple run steps - rename to workflow/job to "unit-n-workflow-tests" + cli/save_gmt.py: fix a typo in the module import + cli/load_data: check `-t smallbaselineApp.cfg` existence and print out error msg + tsview: show the reference index/date info on the slider as the title
1 parent 64ffdef
s1ab_range_bias.py
############################################################
# Program is part of MintPy #
# Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi #
# Author: Zhang Yunjun, Apr 2022 #
############################################################
import os
import numpy as np
from matplotlib import colors, pyplot as plt, ticker
from mintpy.objects import timeseries
from mintpy.utils import plot as pp, readfile, s1_utils, writefile
####################################################################################
def estimate_s1ab_range_bias(ts_file, mask_file=None, safe_list_file=None):
"""Estimate the S1A/B range bias based on the time series file.
Parameters: ts_file - str, path of the time series range offset file
mask_file - str, path of the mask file, e.g., maskResInv.h5 file.
safe_list_file - str, path of the SAFE_files.txt file
Returns: bias_list - list of float32, median bias per subswath in meters
bias_est - 2D np.ndarray in size of (length, width) in float32, bias in meters
mask_list - list of 2D np.ndarray in size of (length, width) in bool
"""
mintpy_dir = os.path.dirname(ts_file)
print(f'Estimating S1A/B range bias from file: {ts_file}')
# read time series
ts_data = readfile.read(ts_file)[0]
length, width = ts_data.shape[-2:]
# read mask
mask_file = mask_file if mask_file else os.path.join(mintpy_dir, 'maskResInv.h5')
if mask_file and os.path.isfile(mask_file):
mask = readfile.read(mask_file)[0]
else:
mask = np.ones((length, width), dtype=np.bool_)
# estimate bias - 2D map
bias_est_poi = s1_utils.estimate_s1ab_bias(
mintpy_dir,
ts_data[:, mask],
safe_list_file=safe_list_file)[0]
if bias_est_poi is None:
print('Exit without estimating S1A/B range bias from the input time series.')
return None, None, None
bias_est = np.ones((length, width), dtype=np.float32) * np.nan
bias_est[mask] = bias_est_poi
# estimate bias - median
geom_file = os.path.join(mintpy_dir, 'inputs', 'geometryRadar.h5')
flag = readfile.read(geom_file, datasetName='height')[0] != 0
mask_list = s1_utils.get_subswath_masks(flag, cut_overlap_in_half=False)[:3]
bias_list = [np.nanmedian(bias_est[x]) for x in mask_list]
print(f'IW1 : {bias_list[0]:.3f} m')
print(f'IW2 : {bias_list[1]:.3f} m')
print(f'IW3 : {bias_list[2]:.3f} m')
return bias_list, bias_est, mask_list
def plot_s1ab_range_bias_est(bias_list, bias_est, mask_list, out_dir=None,
save_fig=False, disp_fig=True):
"""Plot the S1A/B range bias estimation results.
Parameters: bias_list - list of float, mean S1A/B range bias in meters
bias_est - 2D np.ndarray in float32, pixelwised S1A/B range bias in meters
mask_list - list of 2D np.ndarray in bool, mask array for IW1/2/3
out_dir - str, output directory for the plotted figure.
"""
vmin, vmax = 7, 14
font_size = 12
out_dir = out_dir if out_dir else os.getcwd()
## figure 1 - map
fig_size = pp.auto_figure_size(ds_shape=bias_est.shape, disp_cbar=True, print_msg=True)
fig, ax = plt.subplots(figsize=fig_size)
cmap = colors.LinearSegmentedColormap.from_list('magma_t', plt.get_cmap('magma')(np.linspace(0.3, 1.0, 100)))
im = ax.imshow(bias_est*100., cmap=cmap, vmin=vmin, vmax=vmax, interpolation='nearest')
# axis format
ax.tick_params(which='both', direction='out', bottom=True, top=False, left=True, right=False)
ax.set_xlabel('Range [pixel]', fontsize=font_size)
ax.set_ylabel('Azimuth [pixel]', fontsize=font_size)
cbar = fig.colorbar(im, ax=ax)
cbar.set_label('S1A/B range bias [cm]', fontsize=font_size)
# output
if save_fig:
out_fig = os.path.join(out_dir, 's1ab_range_bias_map.pdf')
print('save figure to file', out_fig)
plt.savefig(out_fig, bbox_inches='tight', transparent=True, dpi=300)
## figure 2 - histogram
clist = [cmap((x*100. - vmin) / (vmax - vmin)) for x in bias_list]
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=[6, 2])
for bias, mask, c in zip(bias_list, mask_list, clist):
ax.hist(bias_est[mask].flatten()*100, bins=70, range=(vmin, vmax), density=False, alpha=0.7, color=c)
ax.axvline(bias*100, color='k')
# plot median value
ax.tick_params(which='both', direction='out', bottom=True, top=True, left=True, right=True)
ax.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax.set_xlim(vmin, vmax)
ax.set_xlabel('S1A/B range bias [cm]')
ax.set_ylabel('# of pixels')
fig.tight_layout()
# output
if save_fig:
out_fig = os.path.join(out_dir, 's1ab_range_bias_hist.pdf')
print('save figure to file', out_fig)
plt.savefig(out_fig, bbox_inches='tight', transparent=True, dpi=300)
if disp_fig:
print('showing ....')
plt.show()
else:
plt.close()
return
def write_s1ab_bias_file(bias_file, bias_list, geom_file, force=False):
"""Write estimated S1A/B range bias to HDF5 file.
Parameters: bias_file - str, path to the S1A/B range bias file
bias_list - list of float, constant S1A/B range bias per IW1/2/3
geom_file - str, path to the geometry file
force - bool, overwrite existing bias file.
Returns: bias_file - str, path to the S1A/B range bias file
"""
# run or skip
if os.path.isfile(bias_file) and not force:
print(f'S1Bias file exists in: {bias_file}, skip re-writing.')
return bias_file
# get the list of masks for IW1/2/3
flag = readfile.read(geom_file, datasetName='height')[0] != 0
mask_list = s1_utils.get_subswath_masks(flag, cut_overlap_in_half=False)[:3]
bias_mat = np.zeros(flag.shape, dtype=np.float32) * np.nan
for bias, mask in zip(bias_list, mask_list):
bias_mat[mask] = bias
# write file
atr = readfile.read_attribute(geom_file)
atr['FILE_TYPE'] = 'offset'
atr['UNIT'] = 'm'
print(f'writing S1A/B range bias to file: {bias_file}')
writefile.write(bias_mat, out_file=bias_file, metadata=atr)
return bias_file
def correct_s1ab_range_bias(ts_file, bias_file, ts_cor_file=None, safe_list_file=None):
"""Correct input time series for the S1A/B range bias.
Parameters: ts_file - str, path to the range offset time series file
bias_file - str, path to the S1A/B range bias file
ts_cor_file - str, path to the corrected range offset time series file
safe_list_file - str, path to the SAFE_files.txt file
Returns: ts_cor_file - str, path to the corrected range offset time series file
"""
if not os.path.isfile(bias_file):
msg = f'No bias file found in: {bias_file}!'
msg += '\nRe-run with "--action compute" to generate it.'
raise FileNotFoundError(msg)
date_list = timeseries(ts_file).get_date_list()
num_date = len(date_list)
# date info for Sentinel-1B
mintpy_dir = os.path.dirname(os.path.dirname(bias_file))
s1b_date_list_file = s1_utils.get_s1ab_date_list_file(mintpy_dir, safe_list_file)[1]
s1b_date_list = np.loadtxt(s1b_date_list_file, dtype=str).tolist()
s1b_flag = np.array([x in s1b_date_list for x in date_list], dtype=np.bool_)
# read data
ts_data = readfile.read(ts_file)[0].reshape(num_date, -1)
bias = readfile.read(bias_file)[0].flatten()
# correct bias
mask = ts_data == 0
ts_data[s1b_flag] -= np.tile(bias.reshape(1, -1), (np.sum(s1b_flag), 1))
ts_data[mask] = 0 # Do not change zero value in the input TS file
ts_data[:, np.isnan(bias)] = np.nan # set to nan for pixels with nan in bias file
# write file
if not ts_cor_file:
ts_cor_file = f'{os.path.splitext(ts_file)[0]}_S1Bias.h5'
atr = readfile.read_attribute(ts_file)
length = int(atr['LENGTH'])
width = int(atr['WIDTH'])
writefile.write(ts_data.reshape(num_date, length, width),
out_file=ts_cor_file,
metadata=atr,
ref_file=ts_file)
return ts_cor_file
####################################################################################
def run_s1ab_range_bias(inps):
# matplotlib backend setting
if not inps.disp_fig:
plt.switch_backend('Agg')
# default bias file path
inps.bias_file = os.path.join(os.path.dirname(inps.geom_file), 'S1Bias.h5')
# calculate the S1A/B range bias
if inps.action == 'compute':
if inps.bias_method == 'hardwired':
# option 1 - use the hardwired value from section VII-A in Yunjun et al. (2022)
bias_list = [0.087, 0.106, 0.123] # m
print('Used hardwired S1A/B range bias values from Yunjun et al. (2022):')
print(f'IW1 : {bias_list[0]:.3f} m')
print(f'IW2 : {bias_list[1]:.3f} m')
print(f'IW3 : {bias_list[2]:.3f} m')
else:
# option 2 - estimate from the time series of its dataset itself
# estimate optimal (median) value for each subswath from SenDT156
bias_list, bias_est, mask_list = estimate_s1ab_range_bias(
ts_file=inps.ts_file,
mask_file=inps.mask_file,
safe_list_file=inps.safe_list_file)
# plot the estimation result
if bias_list:
plot_s1ab_range_bias_est(
bias_list,
bias_est,
mask_list,
out_dir=os.path.dirname(inps.ts_file),
save_fig=inps.save_fig,
disp_fig=inps.disp_fig)
# write S1Bias.h5 file
if bias_list:
write_s1ab_bias_file(
bias_file=inps.bias_file,
bias_list=bias_list,
geom_file=inps.geom_file,
force=inps.force)
# correct time series range offset file
elif inps.action == 'correct':
correct_s1ab_range_bias(
ts_file=inps.ts_file,
bias_file=inps.bias_file,
ts_cor_file=inps.ts_cor_file,
safe_list_file=inps.safe_list_file)
return
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...