#!/usr/bin/env python
# coding: utf-8

# ## Geolocation Error Budget for NISAR
# 
# Incidence angle = 42 degree

# In[1]:


get_ipython().run_line_magic('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


# ### Table IV - Range Geolocation Error Budget

# In[2]:


## 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]))


# In[3]:


## 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]))


# In[4]:


## 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]))


# In[5]:


## 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]))


# In[6]:


## 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]))


# In[7]:


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]))


# ### OPERA NI-CSLC: Relative Geolocation Accuracy (using ALOS-2 as a proxy)

# In[8]:


## 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))


# ### OPERA NI-CSLC: Absolute Geolocation Accuracy in Range Direction (w/o DEM error)

# In[9]:


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]))


# ### OPERA NI-CSLC: Absolute Geolocation Accuracy in Azimuth Direction (w/o DEM error)
# 
# We do NOT know due to the lack of TEC gradient/slope knowledge!

# In[ ]: