stripmapApp.py
¶Search and download ALOS data from ASF
DEM
cd ~/data/geolocation/HawaiiAlosDT598
mkdir -p DEM; cd DEM
dem.py -a stitch -b 18 21 -156 -154 -r -s 1 -c -f
stripmapApp.py
stripmapApp.py
¶Run stripmapApp.py stripmapApp.xml --end=misregistration
Coregister using geometry without refinement.
misreg/misreg_rg.xml
and misreg/misreg_az.xml
file.stripmapApp.py stripmapApp.xml --start=refined_resample
Coregister using geometry with ancillary-based refinement
misreg/misreg_rg.xml
file. # 0.8559stripmapApp.py stripmapApp.xml --start=refined_resample
%matplotlib inline
import os
import numpy as np
import datetime as dt
from scipy import interpolate
from osgeo import gdal
from matplotlib import pyplot as plt
from mintpy.objects import sensor
from mintpy.utils import readfile, utils as ut
from mintpy import tropo_pyaps3 as tropo
import pyaps3 as pa
import pysolid
from tools.simulation import iono
proj_dir = os.path.expanduser('~/data/geolocation/HawaiiAlosDT598/')
os.chdir(proj_dir)
print('Go to directory:', proj_dir)
# ALOS info
c = 299792458
f0 = 1269999750.0604727
rg_pixel_size = sensor.SENSOR_DICT['alos']['range_pixel_size']['stripmap_FBD']
Go to directory: /Users/yunjunz/data/geolocation/HawaiiAlosDT598/
# geometry: lat/lon/inc_angle
lat_file = os.path.join(proj_dir, 'geometry', 'lat.rdr')
lon_file = os.path.join(proj_dir, 'geometry', 'lon.rdr')
los_file = os.path.join(proj_dir, 'geometry', 'los.rdr')
hgt_file = os.path.join(proj_dir, 'geometry', 'z.rdr')
lat = np.nanmedian(readfile.read(lat_file)[0])
lon = np.nanmedian(readfile.read(lon_file)[0])
hgt = np.nanmedian(readfile.read(hgt_file)[0])
inc_angle = np.nanmedian(readfile.read(los_file, datasetName='incidenceAngle')[0])
az_angle = np.nanmedian(readfile.read(los_file, datasetName='azimuthAngle')[0])
## Ionosphere
# time:
tmid = dt.datetime.fromisoformat('2011-03-06T20:32:43.389354')
tmin = tmid.hour * 60 + tmid.minute + tmid.second / 60
# geometry: lat/lon/inc_angle at the ionosphere shell
iono_height = 450e3
iono_inc_angle = iono.incidence_angle_ground2iono(inc_angle, iono_height=iono_height)
iono_lat, iono_lon = iono.lalo_ground2iono(lat, lon, inc_angle=inc_angle, az_angle=az_angle)
# geometry: JHR GIM grid
lats = np.arange(-90, 90, step=1) + 0.5
lons = np.arange(-180, 180, step=1) + 0.5
mins = np.arange(0, 24*60, step=15) + 7.5
time_ind = np.argmin(np.abs(mins - tmin))
# calculation/interpolation
tec_vals = []
fnames = [os.path.join(proj_dir, 'GIM', x) for x in ['jpli0190.11i.nc', 'jpli0650.11i.nc']]
for fname in fnames:
ds = gdal.Open(fname, gdal.GA_ReadOnly)
tec_map = gdal.Open(ds.GetSubDatasets()[0][0]).ReadAsArray()
tec_val = interpolate.interp2d(lons, lats, tec_map[time_ind, :, :], kind='linear')(iono_lon, iono_lat)[0]
tec_vals.append(tec_val)
tec_diff = np.diff(tec_vals)[0]
iono_off = iono.vtec2range_delay(tec_diff, iono_inc_angle, f0)
print('ionospheric range offset: {:.3f} m / {:.3f} pixel'.format(iono_off, iono_off / rg_pixel_size))
ionospheric range offset: 5.811 m / 0.620 pixel
## Solid earth tide
dt_obj0 = dt.datetime(2011, 1, 19, 20, 32, 44)
dt_obj1 = dt.datetime(2011, 3, 6, 20, 32, 44) + dt.timedelta(seconds=15)
step_sec = 15
# calculate
(dt_out,
tide_e,
tide_n,
tide_u) = pysolid.calc_solid_earth_tides_point(lat, lon, dt_obj0, dt_obj1, step_sec=step_sec, verbose=False)
tide_e = tide_e[-1] - tide_e[0]
tide_n = tide_n[-1] - tide_n[0]
tide_u = tide_u[-1] - tide_u[0]
tide_los = ut.enu2los(tide_e, tide_n, tide_u, inc_angle=inc_angle, az_angle=az_angle)
tide_off = np.abs(tide_los)
print('SET range offset: {:.3f} m / {:.3f} pixel'.format(tide_off, tide_off / rg_pixel_size))
PYSOLID: calculate solid Earth tides in east/north/up direction PYSOLID: lot/lon: 19.50479061273407/-154.98311663649145 degree PYSOLID: start UTC: 2011-01-19T20:32:44 PYSOLID: end UTC: 2011-03-06T20:32:59 PYSOLID: time step: 15 seconds SET range offset: 0.005 m / 0.001 pixel
## Troposphere
date_list = ['20110119', '20110306']
snwe = [15, 25, -160, -150]
grib_dir = os.path.join(os.path.expandvars('$WEATHER_DIR'), 'ERA5')
grib_files = tropo.get_grib_filenames(date_list=date_list, hour=21, model='ERA5', grib_dir=grib_dir, snwe=snwe)
grib_files = tropo.dload_grib_files(grib_files, tropo_model='ERA5', snwe=snwe)
tropo_los = np.zeros(2, dtype=np.float32)
for i, grib_file in enumerate(grib_files):
tropo_los[i] = tropo.get_delay(grib_file, tropo_model='ERA5', delay_type='comb',
dem=np.array(hgt, dtype=np.float32).reshape(1,1),
inc=np.array(inc_angle, dtype=np.float32).reshape(1,1),
lat=np.array(lat, dtype=np.float32).reshape(1,1),
lon=np.array(lon, dtype=np.float32).reshape(1,1))
tropo_off = np.abs(np.diff(tropo_los))[0]
print('tropospheric range offset: {:.3f} m / {:.3f} pixel'.format(tropo_off, tropo_off / rg_pixel_size))
------------------------------------------------------------------------------ downloading weather model data using PyAPS ... common file size: 386280 bytes number of grib files existed : 2 number of grib files to download: 0 ------------------------------------------------------------------------------ tropospheric range offset: 0.052 m / 0.005 pixel
rg_off = iono_off + tide_off + tropo_off
print('total range offset: {:.4f} m / {:.4f} pixel'.format(rg_off, rg_off / rg_pixel_size))
print('temporal baseline: {} days'.format((dt_obj1 - dt_obj0).days))
print('relative accuracy in range direction: {:.2f} pixel'.format(0.9896 - rg_off / rg_pixel_size))
print('relative accuracy in azimuth direction: {:.2f} pixel'.format(0.5258))
total range offset: 5.8680 m / 0.6263 pixel temporal baseline: 46 days relative accuracy in range direction: 0.36 pixel relative accuracy in azimuth direction: 0.53 pixel