Incidence angle = 42 degree
%matplotlib inline
import os
import numpy as np
from matplotlib import pyplot as plt
from mintpy.utils import readfile, utils as ut
from mintpy.simulation import iono
inc_angle = 42 / 180 * np.pi
## ionosphere
# residual ionosphere based on the bias/uncertainty of GIM
# reference: Hernandez et al. (2009, Table I)
# note that we are using JHR GIM rather than JLR GIM, which from our study here at least, JHR is significantly better than JLR.
tec_jlr = np.array([0.72, 4.49], dtype=np.float32) # mean and STD in TECUm
iono_delay_L = iono.vtec2range_delay(tec_jlr, inc_angle=42, freq=iono.SAR_BAND['L'])
iono_delay_S = iono.vtec2range_delay(tec_jlr, inc_angle=42, freq=iono.SAR_BAND['S'])
print('Iono: residual GIM: {:.3f} +/- {:.3f} m'.format(iono_delay_L[0], iono_delay_L[1]))
# worse case scenario spatially based on total TEC and the map of snapshot of topTECs
top_tec = np.array([6.6, 1.8], dtype=np.float32)
top_tec_delay_L = iono.vtec2range_delay(top_tec, inc_angle=42, freq=iono.SAR_BAND['L'])
top_tec_delay_S = iono.vtec2range_delay(top_tec, inc_angle=42, freq=iono.SAR_BAND['S'])
print('Iono: topside TEC: {:.3f} +/- {:.3f} m'.format(top_tec_delay_L[0], top_tec_delay_L[1]))
Iono: residual GIM: 0.223 +/- 1.206 m Iono: topside TEC: 1.739 +/- 0.517 m
## residual troposphere
# mean: +/- 3 cm in zenith direction (Yu et al., 2021, Fig. 1a1)
# STD: +/- 5 cm in LOS direction (Fattahi & Amelung, 2015, sec 4.2.2)
tropo = np.array([0.03 / np.cos(inc_angle), 0.05], dtype=np.float32)
print('Tropo: {:.3f} +/- {:.3f} m'.format(tropo[0], tropo[1]))
Tropo: 0.040 +/- 0.050 m
## uncompensated tidal effects
OTL_ENU = np.array([0.0100, 0.0100, 0.0500], dtype=np.float32) # reference: Martens et al. (2016); Gisinger (0.01, 0.01, 0.05)
PT_ENU = np.array([0.0035, 0.0035, 0.0125], dtype=np.float32) # reference: Petit & Luzum (2010)
ATM_ENU = np.array([0.0005, 0.0005, 0.0050], dtype=np.float32) # reference: Dong et al. (2002)
OTL_LOS = ut.enu2los(OTL_ENU[0], OTL_ENU[1], OTL_ENU[2], inc_angle=inc_angle*180/np.pi, az_angle=-102)
PT_LOS = ut.enu2los( PT_ENU[0], PT_ENU[1], PT_ENU[2], inc_angle=inc_angle*180/np.pi, az_angle=-102)
ATM_LOS = ut.enu2los(ATM_ENU[0], ATM_ENU[1], ATM_ENU[2], inc_angle=inc_angle*180/np.pi, az_angle=-102)
EM_LOS = np.array([0, np.sqrt(OTL_LOS**2 + PT_LOS**2 + ATM_LOS**2)], dtype=np.float32)
print('Earth motion: {:.3f} +/- {:.3f} m'.format(EM_LOS[0], EM_LOS[1]))
Earth motion: 0.000 +/- 0.044 m
## orbital error
# reference: Peter et al. (2017)
orb_s1 = np.array([0.005, 0.030], dtype=np.float32)
print('Orbital error: {:.3f} +/- {:.3f} m'.format(orb_s1[0], orb_s1[1]))
Orbital error: 0.005 +/- 0.030 m
## DEM error
# reference: 1. NISAR handbook; 2. Sansosti et al. (2006, Eq. 30); 3. Jung et al. (2019, Eq. 8); 4. Fahrland et al. (2020)
# Copernicus global DEM accuracy: https://spacedata.copernicus.eu/web/cscda/dataset-details?articleId=394198
rg = 980e3 # m, NISAR
bperp_max = 350 # m
dem_err = 4 # m, absolute vertical accuracy
dem_LOS_rel = dem_err * bperp_max / (rg * np.sin(inc_angle)) * 1
dem_LOS_abs = dem_err / np.tan(inc_angle)
dem_LOS = np.array([dem_LOS_abs, dem_LOS_rel], dtype=np.float32)
print('DEM error: {:.3f} +/- {:.3f} m'.format(dem_LOS[0], dem_LOS[1]))
DEM error: 4.442 +/- 0.002 m
geoloc_L = np.sqrt(iono_delay_L**2 + top_tec_delay_L**2 + tropo**2 + EM_LOS**2 + orb_s1**2 + dem_LOS**2)
geoloc_L[0] = iono_delay_L[0] + top_tec_delay_L[0] + tropo[0] + EM_LOS[0] + orb_s1[0] + dem_LOS[0]
print('NISAR L-band:')
print(' Iono residual GIM: {:.3f} +/- {:.3f} m'.format(iono_delay_L[0], iono_delay_L[1]))
print(' Iono topside TEC: {:.3f} +/- {:.3f} m'.format(top_tec_delay_L[0], top_tec_delay_L[1]))
print(' Tropo : {:.3f} +/- {:.3f} m'.format(tropo[0], tropo[1]))
print(' Tidal residual : {:.3f} +/- {:.3f} m'.format(EM_LOS[0], EM_LOS[1]))
print(' Orbital error : {:.3f} +/- {:.3f} m'.format(orb_s1[0], orb_s1[1]))
print(' DEM error : {:.3f} +/- {:.3f} m'.format(dem_LOS[0], dem_LOS[1]))
print(' Overall (L-band) : {:.3f} +/- {:.3f} m'.format(geoloc_L[0], geoloc_L[1]))
geoloc_S = np.sqrt(iono_delay_S**2 + top_tec_delay_S**2 + tropo**2 + EM_LOS**2 + orb_s1**2 + dem_LOS**2)
geoloc_S[0] = iono_delay_S[0] + top_tec_delay_S[0] + tropo[0] + EM_LOS[0] + orb_s1[0] + dem_LOS[0]
print('NISAR S-band:')
print(' Iono residual GIM: {:.3f} +/- {:.3f} m'.format(iono_delay_S[0], iono_delay_S[1]))
print(' Iono topside TEC: {:.3f} +/- {:.3f} m'.format(top_tec_delay_S[0], top_tec_delay_S[1]))
print(' ...')
print(' Overall (S-band) : {:.3f} +/- {:.3f} m'.format(geoloc_S[0], geoloc_S[1]))
NISAR L-band: Iono residual GIM: 0.223 +/- 1.206 m Iono topside TEC: 1.739 +/- 0.517 m Tropo : 0.040 +/- 0.050 m Tidal residual : 0.000 +/- 0.044 m Orbital error : 0.005 +/- 0.030 m DEM error : 4.442 +/- 0.002 m Overall (L-band) : 6.449 +/- 1.314 m NISAR S-band: Iono residual GIM: 0.037 +/- 0.215 m Iono topside TEC: 0.307 +/- 0.091 m ... Overall (S-band) : 4.832 +/- 0.244 m
## range
rg_reso = 6.25 # m for NISAR-L 24 MHz
rg_rmse = 0.56 # m from ALOS-2 desc track 23 over Kyushu, Japan
## azimuth
az_reso = 7 # m for NISAR-L (NISAR handbook)
mintpy_dir = os.path.expanduser('~/data/geolocation/KyushuAlos2DT23/mintpy_offset')
ts_file = os.path.join(mintpy_dir, 'timeseriesAz.h5') # timeseriesRg.h5 or timeseriesAz.h5
# calculate az geoloc rmse
mask_file = os.path.join(mintpy_dir, 'maskResInv.h5')
ts_data = readfile.read(ts_file)[0]; ts_data = ts_data.reshape(ts_data.shape[0], -1)
mask = readfile.read(mask_file)[0].flatten()
ts_data[:, mask == 0] = np.nan
ts_med = np.nanmedian(ts_data, axis=-1)
ts_med -= np.nanmedian(ts_med)
az_rmse = ut.root_mean_sq_error(ts_med)
print('Range : {:.3f} m --> {:.2f} pixel'.format(rg_rmse, rg_rmse / rg_reso))
print('Azimuth: {:.3f} m --> {:.2f} pixel'.format(az_rmse, az_rmse / az_reso))
Range : 0.560 m --> 0.09 pixel Azimuth: 0.810 m --> 0.12 pixel
geoloc = np.sqrt(iono_delay_L**2 + top_tec_delay_L**2 + tropo**2 + EM_LOS**2 + orb_s1**2)
print('Range w/o Topo: {:.1f} m'.format(geoloc[0]))
Range w/o Topo: 1.8 m
We do NOT know due to the lack of TEC gradient/slope knowledge!