From 6c49eaa1e68fc2ed11b321168609d3a53b64a80e Mon Sep 17 00:00:00 2001 From: JingChen-Thu Date: Wed, 22 Jul 2026 11:54:55 +0800 Subject: [PATCH 1/6] add rotation in crust1.0 model --- pytomoatt/io/crustmodel.py | 32 +++++++++++++++++++++++--------- pytomoatt/model.py | 6 ++++-- 2 files changed, 27 insertions(+), 11 deletions(-) diff --git a/pytomoatt/io/crustmodel.py b/pytomoatt/io/crustmodel.py index f7fa937..cb5aa1e 100644 --- a/pytomoatt/io/crustmodel.py +++ b/pytomoatt/io/crustmodel.py @@ -2,6 +2,7 @@ from os.path import dirname, abspath, join from ..utils.common import init_axis from ..setuplog import SetupLog +from ..utils.rotate import rtp_rotation_reverse import pickle import sys from tqdm import tqdm @@ -40,7 +41,7 @@ def __init__(self, fname=join(dirname(dirname(abspath(__file__))), 'data', 'crus self.points_dict = pickle.load(f) self.log = SetupLog() - def griddata(self, min_max_dep, min_max_lat, min_max_lon, n_rtp, type='vp'): + def griddata(self, min_max_dep, min_max_lat, min_max_lon, n_rtp, type='vp', rotate=None): """Linearly interpolate velocity into regular grids :param min_max_dep: min and max depth, ``[min_dep, max_dep]`` @@ -63,25 +64,38 @@ def griddata(self, min_max_dep, min_max_lat, min_max_lon, n_rtp, type='vp'): else: self.log.Modellog.error(f"Velocity type {type} not supported in CRUST1.0 model") sys.exit(1) + self.dd, self.tt, self.pp, _, _, _, = init_axis( min_max_dep, min_max_lat, min_max_lon, n_rtp ) + tt_2d, pp_2d = np.meshgrid(self.tt, self.pp, indexing='ij') + + # rotate reversely, from computational grid to physical grid + if rotate is not None: + central_lat = rotate[0] + central_lon = rotate[1] + rotation_angle = rotate[2] + tt_2d, pp_2d = rtp_rotation_reverse(tt_2d, pp_2d, central_lat, central_lon, rotation_angle) + # Grid data self.log.Modellog.info('Grid data, please wait for a few minutes') vel = np.zeros(n_rtp) with tqdm(total=self.n_rtp[1] * self.n_rtp[2], desc='Gridding') as pbar: for ilat in range(self.n_rtp[1]): - new_lat = self.tt[ilat] - idx_lat_left, ratio_lat = degree_to_idx_and_ratio(new_lat) - idx_lat_right = idx_lat_left + 1 - if idx_lat_left == -1: - idx_lat_left = 0 - idx_lat_right = 1 - for ilon in range(self.n_rtp[2]): pbar.update(1) - new_lon = self.pp[ilon] + + # latitude index and ratio + new_lat = tt_2d[ilat, ilon] + idx_lat_left, ratio_lat = degree_to_idx_and_ratio(new_lat) + idx_lat_right = idx_lat_left + 1 + if idx_lat_left == -1: + idx_lat_left = 0 + idx_lat_right = 1 + + # longitude index and ratio + new_lon = pp_2d[ilat, ilon] idx_lon_left, ratio_lon = degree_to_idx_and_ratio(new_lon) idx_lon_right = idx_lon_left + 1 if idx_lon_left == -1: # between -179.5 and +179.5 diff --git a/pytomoatt/model.py b/pytomoatt/model.py index 14c5102..b41d082 100644 --- a/pytomoatt/model.py +++ b/pytomoatt/model.py @@ -102,18 +102,20 @@ def to_xarray(self): ) return dataset - def grid_data_crust1(self, type='vp'): + def grid_data_crust1(self, type='vp', rotate=None): """Grid data from CRUST1.0 model :param type: Specify velocity type of ``vp`` or ``vs``, defaults to 'vp' :type type: str, optional + :rotate: Rotation parameters [theta0, phi0, psi], defaults to None """ cm = CrustModel() self.vel = cm.griddata( self.min_max_dep, self.min_max_lat, self.min_max_lon, - self.n_rtp, type=type + self.n_rtp, type=type, + rotate=rotate ) def grid_data_ascii(self, model_fname:str, **kwargs): From d6ac9f3683b711d3060976f1a23e9ace86ef87cf Mon Sep 17 00:00:00 2001 From: JingChen-Thu Date: Wed, 22 Jul 2026 12:05:04 +0800 Subject: [PATCH 2/6] exclude zeta in Checker --- pytomoatt/checkerboard.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/pytomoatt/checkerboard.py b/pytomoatt/checkerboard.py index 589c577..309a73b 100644 --- a/pytomoatt/checkerboard.py +++ b/pytomoatt/checkerboard.py @@ -22,7 +22,10 @@ def __init__(self, model_fname:str, para_fname='input_params.yml') -> None: self.vel = f['vel'][:] self.eta = f['eta'][:] self.xi = f['xi'][:] - self.zeta = f['zeta'][:] + try: # some model may not have zeta + self.zeta = f['zeta'][:] + except: + pass self._init_axis() def _init_axis(self): From 8079307153946f4ea50c80ba7d7376ac3161b146 Mon Sep 17 00:00:00 2001 From: JingChen-Thu Date: Wed, 22 Jul 2026 12:06:09 +0800 Subject: [PATCH 3/6] exclude zeta in Checker --- pytomoatt/checkerboard.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pytomoatt/checkerboard.py b/pytomoatt/checkerboard.py index 309a73b..6724d51 100644 --- a/pytomoatt/checkerboard.py +++ b/pytomoatt/checkerboard.py @@ -25,7 +25,7 @@ def __init__(self, model_fname:str, para_fname='input_params.yml') -> None: try: # some model may not have zeta self.zeta = f['zeta'][:] except: - pass + self.zeta = np.zeros_like(self.vel) self._init_axis() def _init_axis(self): From 46dadf484285737b2f6963a48531194ba73c63de Mon Sep 17 00:00:00 2001 From: Mijian Xu Date: Wed, 22 Jul 2026 14:17:03 +0800 Subject: [PATCH 4/6] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- pytomoatt/model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pytomoatt/model.py b/pytomoatt/model.py index b41d082..ef1adc8 100644 --- a/pytomoatt/model.py +++ b/pytomoatt/model.py @@ -107,7 +107,7 @@ def grid_data_crust1(self, type='vp', rotate=None): :param type: Specify velocity type of ``vp`` or ``vs``, defaults to 'vp' :type type: str, optional - :rotate: Rotation parameters [theta0, phi0, psi], defaults to None + :param rotate: Rotation parameters [central_lat, central_lon, rotation_angle] in degrees, defaults to None """ cm = CrustModel() self.vel = cm.griddata( From e4c6effa4a4fbe66195cf7b4f9fc30f5b513aabd Mon Sep 17 00:00:00 2001 From: Mijian Xu Date: Wed, 22 Jul 2026 14:17:54 +0800 Subject: [PATCH 5/6] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- pytomoatt/io/crustmodel.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/pytomoatt/io/crustmodel.py b/pytomoatt/io/crustmodel.py index cb5aa1e..15bdc66 100644 --- a/pytomoatt/io/crustmodel.py +++ b/pytomoatt/io/crustmodel.py @@ -73,9 +73,11 @@ def griddata(self, min_max_dep, min_max_lat, min_max_lon, n_rtp, type='vp', rota # rotate reversely, from computational grid to physical grid if rotate is not None: - central_lat = rotate[0] - central_lon = rotate[1] - rotation_angle = rotate[2] + try: + central_lat, central_lon, rotation_angle = rotate + except (TypeError, ValueError): + self.log.Modellog.error("rotate must be a 3-item sequence: [central_lat, central_lon, rotation_angle]") + sys.exit(1) tt_2d, pp_2d = rtp_rotation_reverse(tt_2d, pp_2d, central_lat, central_lon, rotation_angle) # Grid data From a577245eb5dd3eaf9bcbfe02935c89051bc061e6 Mon Sep 17 00:00:00 2001 From: Mijian Xu Date: Wed, 22 Jul 2026 15:18:23 +0800 Subject: [PATCH 6/6] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- pytomoatt/checkerboard.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pytomoatt/checkerboard.py b/pytomoatt/checkerboard.py index 6724d51..a65db5c 100644 --- a/pytomoatt/checkerboard.py +++ b/pytomoatt/checkerboard.py @@ -22,9 +22,9 @@ def __init__(self, model_fname:str, para_fname='input_params.yml') -> None: self.vel = f['vel'][:] self.eta = f['eta'][:] self.xi = f['xi'][:] - try: # some model may not have zeta + if 'zeta' in f: # some model may not have zeta self.zeta = f['zeta'][:] - except: + else: self.zeta = np.zeros_like(self.vel) self._init_axis()