Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

plate_motion: add --plate option to use built-in tables #906

Merged
merged 2 commits into from
Nov 17, 2022
Merged

plate_motion: add --plate option to use built-in tables #906

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
2a5bb67
move add_attribute.py to legacy dir
yunjunz Nov 13, 2022
ebf4ca9
plate_motion: add --plate option to use built-in tables
yunjunz Nov 12, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
59 changes: 51 additions & 8 deletions src/mintpy/cli/plate_motion.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,11 +6,37 @@
############################################################


import collections
import os
import sys

from mintpy.utils.arg_utils import create_argument_parser

# ITRF2014-PMM defined in Altamimi et al. (2017)
Tag = collections.namedtuple('Tag', 'name num_site omega_x omega_y omega_z omega wrms_e wrms_n')
ITRF2014_PMM = {
'ANTA' : Tag('Antartica' , 7, -0.248, -0.324, 0.675, 0.219, 0.20, 0.16),
'ARAB' : Tag('Arabia' , 5, 1.154, -0.136, 1.444, 0.515, 0.36, 0.43),
'AUST' : Tag('Australia' , 36, 1.510, 1.182, 1.215, 0.631, 0.24, 0.20),
'EURA' : Tag('Eurasia' , 97, -0.085, -0.531, 0.770, 0.261, 0.23, 0.19),
'INDI' : Tag('India' , 3, 1.154, -0.005, 1.454, 0.516, 0.21, 0.21),
'NAZC' : Tag('Nazca' , 2, -0.333, -1.544, 1.623, 0.629, 0.13, 0.19),
'NOAM' : Tag('N. America' , 72, 0.024, -0.694, -0.063, 0.194, 0.23, 0.28),
'NUBI' : Tag('Nubia' , 24, 0.099, -0.614, 0.733, 0.267, 0.28, 0.36),
'PCFC' : Tag('Pacific' , 18, -0.409, 1.047, -2.169, 0.679, 0.36, 0.31),
'SOAM' : Tag('S. America' , 30, -0.270, -0.301, -0.140, 0.119, 0.34, 0.35),
'SOMA' : Tag('Somalia' , 3, -0.121, -0.794, 0.884, 0.332, 0.32, 0.30),
}
PMM_UNIT = {
'omega' : 'deg/Ma', # degree per megayear or one-million-year
'omega_x' : 'mas/yr', # milli-arcsecond per year
'omega_y' : 'mas/yr', # milli-arcsecond per year
'omega_z' : 'mas/yr', # milli-arcsecond per year
'wrms_e' : 'mm/yr', # milli-meter per year, weighted root mean scatter
'wrms_n' : 'mm/yr', # milli-meter per year, weighted root mean scatter
}


######################################### Usage ##############################################
REFERENCE = """reference:
Stephenson, O. L., Liu, Y. K., Yunjun, Z., Simons, M., Rosen, P. and Xu, X., (2022),
Expand All @@ -31,10 +57,14 @@
"""

EXAMPLE = """example:
# Use build-in plate motion model of Table 1 from Altamimi et al. (2017)
plate_motion.py -g inputs/geometryGeo.h5 --plate ARAB
plate_motion.py -g inputs/geometryRadar.h5 --plate EURA

# Cartesian form of Euler pole rotation in [wx, wy, wz] in unit of mas/year [milli arc second per year]
# e.g., Arabia plate in ITRF14-PMM (Table 1 in Altamimi et al., 2017)
plate_motion.py -g inputs/geometryGeo.h5 --om-cart 1.154 -0.136 1.444 -v velocity.h5
plate_motion.py -g inputs/geometryRadar.h5 --om-cart 1.154 -0.136 1.444
plate_motion.py -g inputs/geometryGeo.h5 --om-cart 1.154 -0.136 1.444 -v velocity.h5
plate_motion.py -g inputs/geometryRadar.h5 --om-cart 1.154 -0.136 1.444 --comp en2az

# Simple constant local ENU translation (based on one GNSS vector) in [ve, vn, vu] in unit of m/year
Expand Down Expand Up @@ -68,21 +98,26 @@ def create_parser(subparsers=None):
help='Input velocity file to be corrected.')
parser.add_argument('-o', '--output', dest='cor_vel_file', type=str,
help='Output velocity file after the correction, default: add "_ITRF14" suffix.')

# plate motion calculation
parser.add_argument('--comp', dest='pmm_comp', choices={'enu2los', 'en2az'}, default='enu2los',
help='Convert the ENU components into the component of interest for radar (default: %(default)s).')
parser.add_argument('--step','--pmm-step', dest='pmm_step', type=float, default=10.,
help='Ground step/resolution in km for computing PMM to ENU velocity (default: %(default)s).')

# plate motion configurations
pmm = parser.add_mutually_exclusive_group(required=True)
pmm.add_argument('--om-cart', dest='omega_cart', type=float, nargs=3, metavar=('WX', 'WY', 'WZ'), default=None,
pmm = parser.add_argument_group('plate motion model')
pmg = pmm.add_mutually_exclusive_group(required=True)
pmg.add_argument('--plate', dest='plate_name', type=str, choices=ITRF2014_PMM.keys(), default=None,
help='Tectonic plate in ITRF14 (Table 1 in Altamimi et al., 2017) as below:\n' +
'\n'.join(f'{key} for {val.name}' for key, val in ITRF2014_PMM.items()))
pmg.add_argument('--om-cart', dest='omega_cart', type=float, nargs=3, metavar=('WX', 'WY', 'WZ'), default=None,
help='Cartesian form of Euler Pole rotation (unit: mas/yr) (default: %(default)s).')
pmm.add_argument('--om-sph', dest='omega_sph', type=float, nargs=3, metavar=('LAT', 'LON', 'W'), default=None,
pmg.add_argument('--om-sph', dest='omega_sph', type=float, nargs=3, metavar=('LAT', 'LON', 'W'), default=None,
help='Spherical form of Euler Pole rotation (unit: deg, deg, deg/Ma) (default: %(default)s).')
pmm.add_argument('--enu', dest='const_vel_enu', type=float, nargs=3, metavar=('VE', 'VN', 'VU'), default=None,
pmg.add_argument('--enu', dest='const_vel_enu', type=float, nargs=3, metavar=('VE', 'VN', 'VU'), default=None,
help='Constant local ground translation (unit: m/year) (default: %(default)s).')

parser.add_argument('--step','--pmm-step', dest='pmm_step', type=float, default=10.,
help='Ground step/resolution in km for computing PMM to ENU velocity (default: %(default)s).')

return parser


Expand All @@ -102,6 +137,14 @@ def cmd_line_parse(iargs=None):
vbase = os.path.splitext(inps.vel_file)[0]
inps.cor_vel_file = os.path.abspath(f'{vbase}_ITRF14.h5')

# check: --plate option (and convert to --om-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)

return inps


Expand Down
0 src/mintpy/add_attribute.py → src/mintpy/legacy/add_attribute.py
View file
Open in desktop
File renamed without changes.
5 changes: 5 additions & 0 deletions src/mintpy/objects/euler_pole.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,8 @@
# Palo Alto. DOI: 10.4236/ojapps.2015.54016. Page 145-156.
# Navipedia, Transformations between ECEF and ENU coordinates. [Online].
# https://gssc.esa.int/navipedia/index.php/Transformations_between_ECEF_and_ENU_coordinates
# Goudarzi, M. A., Cocard, M. & Santerre, R. (2014), EPC: Matlab software to estimate Euler
# pole parameters, GPS Solutions, 18, 153-162, doi: 10.1007/s10291-013-0354-4


import numpy as np
Expand Down Expand Up @@ -178,6 +180,9 @@ def get_velocity_xyz(self, lat, lon, alt=0.0, ellps=True, print_msg=True):
poi_shape = lat.shape if isinstance(lat, np.ndarray) else None

# convert lat/lon into x/y/z
# Note: the conversion assumes either a spherical or spheroidal Earth, tests show that
# using a ellipsoid as defined in WGS84 produce results closer to the UNAVCO website
# calculator, which also uses the WGS84 ellipsoid.
if ellps:
if print_msg:
print('assume a spheroidal Earth as defined in WGS84')
Expand Down
35 changes: 5 additions & 30 deletions src/mintpy/plate_motion.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,12 +10,7 @@
# 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.
#
# To-Do List (updated 2022.10.12 Yuan-Kai Liu):
# + Use built-in PMM table ITRF2014_PMM for easier/automatic user input string input


import collections

import numpy as np
from skimage.transform import resize
Expand All @@ -25,31 +20,6 @@
from mintpy.objects.resample import resample
from mintpy.utils import readfile, utils as ut, writefile

# ITRF2014-PMM defined in Altamimi et al. (2017)
Tag = collections.namedtuple('Tag', 'name num_site omega_x omega_y omega_z omega wrms_e wrms_n')
ITRF2014_PMM = {
'ANTA' : Tag('Antartica' , 7, -0.248, -0.324, 0.675, 0.219, 0.20, 0.16),
'ARAB' : Tag('Arabia' , 5, 1.154, -0.136, 1.444, 0.515, 0.36, 0.43),
'AUST' : Tag('Australia' , 36, 1.510, 1.182, 1.215, 0.631, 0.24, 0.20),
'EURA' : Tag('Eurasia' , 97, -0.085, -0.531, 0.770, 0.261, 0.23, 0.19),
'INDI' : Tag('India' , 3, 1.154, -0.005, 1.454, 0.516, 0.21, 0.21),
'NAZC' : Tag('Nazca' , 2, -0.333, -1.544, 1.623, 0.629, 0.13, 0.19),
'NOAM' : Tag('N. America' , 72, 0.024, -0.694, -0.063, 0.194, 0.23, 0.28),
'NUBI' : Tag('Nubia' , 24, 0.099, -0.614, 0.733, 0.267, 0.28, 0.36),
'PCFC' : Tag('Pacific' , 18, -0.409, 1.047, -2.169, 0.679, 0.36, 0.31),
'SOAM' : Tag('S. America' , 30, -0.270, -0.301, -0.140, 0.119, 0.34, 0.35),
'SOMA' : Tag('Somalia' , 3, -0.121, -0.794, 0.884, 0.332, 0.32, 0.30),
}
PMM_UNIT = {
'omega' : 'deg/Ma', # degree per megayear or one-million-year
'omega_x' : 'mas/yr', # milli-arcsecond per year
'omega_y' : 'mas/yr', # milli-arcsecond per year
'omega_z' : 'mas/yr', # milli-arcsecond per year
'wrms_e' : 'mm/yr', # milli-meter per year, weighted root mean scatter
'wrms_n' : 'mm/yr', # milli-meter per year, weighted root mean scatter
}


####################################### Major Function ###########################################

def calc_plate_motion(geom_file, omega_cart=None, omega_sph=None, const_vel_enu=None,
Expand Down Expand Up @@ -143,6 +113,9 @@ def calc_plate_motion(geom_file, omega_cart=None, omega_sph=None, const_vel_enu=
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)
Expand Down Expand Up @@ -197,6 +170,7 @@ def calc_plate_motion(geom_file, omega_cart=None, omega_sph=None, const_vel_enu=
def run_plate_motion(inps):
"""Calculate and/or correct for the rigid motion from tectonic plates."""

# calculate plate motion
calc_plate_motion(
geom_file=inps.geom_file,
omega_cart=inps.omega_cart,
Expand All @@ -208,6 +182,7 @@ def run_plate_motion(inps):
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')
Expand Down
2 changes: 1 addition & 1 deletion tests/objects/euler_pole.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,8 @@
import collections
import math

from mintpy.cli.plate_motion import ITRF2014_PMM
from mintpy.objects.euler_pole import MASY2DMY, EulerPole
from mintpy.plate_motion import ITRF2014_PMM

# validation against the UNAVCO Plate Motion Calculator
# https://www.unavco.org/software/geodetic-utilities/plate-motion-calculator/plate-motion-calculator.html
Expand Down