-
Notifications
You must be signed in to change notification settings - Fork 0
/
plate_motion.py
201 lines (167 loc) · 8.51 KB
/
plate_motion.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
############################################################
# Program is part of MintPy #
# Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi #
# Author: Yuan-Kai Liu, Zhang Yunjun, May 2022 #
############################################################
# Recommend usage:
# from mintpy import plate_motion as pmm
#
# Reference:
# Stephenson, O. L., Liu, Y. K., Yunjun, Z., Simons, M., Rosen, P. and Xu, X., (2022),
# The Impact of Plate Motions on Long-Wavelength InSAR-Derived Velocity Fields,
# Geophys. Res. Lett. 49, e2022GL099835, doi:10.1029/2022GL099835.
import numpy as np
from skimage.transform import resize
from euler_pole import ITRF2014_PMM, EulerPole
# original mintpy dependency
from mintpy.diff import diff_file
from mintpy.objects.resample import resample
from mintpy.utils import readfile, utils as ut, writefile
####################################### Major Function ###########################################
def calc_plate_motion(geom_file, omega_cart=None, omega_sph=None, const_vel_enu=None,
pmm_enu_file=None, pmm_file=None, pmm_comp='enu2los', pmm_step=10.):
"""Estimate LOS motion due to the rigid plate motion (translation and/or rotation).
Parameters: geom_file - str, path to the input geometry file
omega_cart - list or 1D array, Cartesian representation of plate rotation
in [wx, wy, wz] (mas/yr)
omega_sph - list or 1D array, Spherical representation of plate rotation
in [lat, lon, w] (deg, deg, deg/Ma)
const_vel_enu - list or 1D array, a single-vector [ve, vn, vu] (meter/year)
simulating the rigid translation of the ground (e.g., from GNSS)
pmm_enu_file - str, path to the output plate motion in east, north, up direction
pmm_file - str, path to the output plate motion in LOS direction
pmm_comp - str, output PMM in the given component of interest
pmm_step - float, ground resolution for computing Plate rotation to ENU velocity (km)
Returns: ve/vn/vu/vlos - 2D np.ndarray, ridig plate motion in east / north / up / LOS direction
"""
# Get LOS geometry
atr_geo = ut.prepare_geo_los_geometry(geom_file, unit='deg')[2]
shape_geo = [int(atr_geo['LENGTH']), int(atr_geo['WIDTH'])]
## calc plate motion in the region
print('-'*50)
if omega_cart or omega_sph:
print('compute the rigid plate motion defined as an Euler Pole')
# construct Euler Pole object
if omega_cart is not None:
print(f'input omega_cartesian in [wx, wy, wz]: {omega_cart} [mas/yr]')
pole_obj = EulerPole(
wx=omega_cart[0],
wy=omega_cart[1],
wz=omega_cart[2],
unit='mas/yr',
)
else:
print(f'input omega_spherical in [lat, lon, w]: {omega_sph} [deg, deg, deg/Ma]')
pole_obj = EulerPole(
pole_lat=omega_sph[0],
pole_lon=omega_sph[1],
rot_rate=omega_sph[2],
unit='deg/Ma',
)
pole_obj.print_info()
# prepare the coarse grid (for the points of interest)
latc = float(atr_geo['Y_FIRST']) + float(atr_geo['Y_STEP']) * shape_geo[0] / 2
ystep = abs(int(pmm_step * 1000 / (float(atr_geo['Y_STEP']) * 108e3)))
xstep = abs(int(pmm_step * 1000 / (float(atr_geo['X_STEP']) * 108e3 * np.cos(np.deg2rad(latc)))))
ystep, xstep = max(ystep, 5), max(xstep, 5)
lats, lons = ut.get_lat_lon(atr_geo, dimension=2, ystep=ystep, xstep=xstep)
print(f'calculate plate motion on the coarse grid: size = ~{pmm_step} km, shape = {lats.shape}')
# calculate plate motion in ENU at the coarse grid
ve_low, vn_low, vu_low = pole_obj.get_velocity_enu(lats, lons, alt=0.0, ellps=True)
# for debugging purpose
debug_mode = False
if debug_mode:
from matplotlib import pyplot as plt
# calculate plate motion in ECEF (XYZ) coordinates
vx, vy, vz = pole_obj.get_velocity_xyz(lats, lons, alt=0.0, ellps=True)
# plot
fig, axs = plt.subplots(nrows=2, ncols=3, figsize=[12, 6])
vlist = [vx, vy, vz, ve_low, vn_low, vu_low]
titles = ['X', 'Y', 'Z', 'E', 'N', 'U']
for ax, data, title in zip(axs.flatten(), vlist, titles):
im = ax.imshow(data, interpolation='nearest')
fig.colorbar(im, ax=ax)
ax.set_title(title)
fig.tight_layout()
plt.show()
# resample coarse grid back to the initial fine grid
print(f'resample plate motion from corase back to original grid: {lats.shape} -> {shape_geo}'
' via skimage.transform.resize ...')
kwargs = dict(order=1, mode='edge', anti_aliasing=True, preserve_range=True)
ve = resize(ve_low, shape_geo, **kwargs)
vn = resize(vn_low, shape_geo, **kwargs)
vu = resize(vu_low, shape_geo, **kwargs)
elif const_vel_enu:
print(f'compute the rigid plate motion using a single vector (translation): {const_vel_enu}')
ve = const_vel_enu[0] * np.ones(shape_geo, dtype=np.float32)
vn = const_vel_enu[1] * np.ones(shape_geo, dtype=np.float32)
vu = const_vel_enu[2] * np.ones(shape_geo, dtype=np.float32)
else:
raise ValueError('No plate motion configuration (--om-cart/sph or --enu) found!')
# radar-code the plate motion if input geometry is in radar coordinates
atr = readfile.read_attribute(geom_file)
if 'Y_FIRST' not in atr.keys():
print('radar-coding the rigid plate motion in ENU ...')
res_obj = resample(lut_file=geom_file)
res_obj.open()
res_obj.src_meta = atr_geo
res_obj.prepare()
# resample data
box = res_obj.src_box_list[0]
ve = res_obj.run_resample(src_data=ve[box[1]:box[3], box[0]:box[2]])
vn = res_obj.run_resample(src_data=vn[box[1]:box[3], box[0]:box[2]])
vu = res_obj.run_resample(src_data=vu[box[1]:box[3], box[0]:box[2]])
## project Plate motion from ENU to direction of interest, e.g. LOS or az
c0, c1 = pmm_comp.split('2')
print(f'project the ridig plate motion from {c0.upper()} onto {c1.upper()} direction')
los_inc_angle = readfile.read(geom_file, datasetName='incidenceAngle')[0]
los_az_angle = readfile.read(geom_file, datasetName='azimuthAngle')[0]
unit_vec = ut.get_unit_vector4component_of_interest(los_inc_angle, los_az_angle, comp=pmm_comp)
vlos = ( ve * unit_vec[0]
+ vn * unit_vec[1]
+ vu * unit_vec[2])
# save the plate motion model velocity into HDF5 files
# metadata
atr['FILE_TYPE'] = 'velocity'
atr['DATA_TYPE'] = 'float32'
atr['UNIT'] = 'm/year'
for key in ['REF_Y', 'REF_X', 'REF_DATE']:
if key in atr.keys():
atr.pop(key)
if pmm_enu_file:
# dataset
dsDict = {'east' : ve,
'north' : vn,
'up' : vu}
# write
writefile.write(dsDict, out_file=pmm_enu_file, metadata=atr)
if pmm_file:
writefile.write(vlos, out_file=pmm_file, metadata=atr)
return ve, vn, vu, vlos
################################################################################################
def run_plate_motion(inps):
"""Calculate and/or correct for the rigid motion from tectonic plates."""
# check: --plate option (convert plate_name to omega_cart)
if inps.plate_name:
plate = ITRF2014_PMM[inps.plate_name]
inps.omega_cart = [plate.omega_x, plate.omega_y, plate.omega_z]
msg = f'get rotation parameters for {inps.plate_name} plate from Table 1 in Altamimi et al. (2017): '
msg += f'wx, wy, wz = {plate.omega_x}, {plate.omega_y}, {plate.omega_z} mas/yr'
print(msg)
# calculate plate motion
calc_plate_motion(
geom_file=inps.geom_file,
omega_cart=inps.omega_cart,
omega_sph=inps.omega_sph,
const_vel_enu=inps.const_vel_enu,
pmm_enu_file=inps.pmm_enu_file,
pmm_file=inps.pmm_file,
pmm_comp=inps.pmm_comp,
pmm_step=inps.pmm_step,
)
# correct plate motion from input velocity
if inps.vel_file and inps.pmm_file and inps.cor_vel_file:
print('-'*50)
print('Correct input velocity for the rigid plate motion')
diff_file(inps.vel_file, [inps.pmm_file], inps.cor_vel_file)
return