Devel

.github/workflows/build-test-conda.yml

@@ -45,7 +45,7 @@ jobs:
         run: |
           conda install hatchling 
           python -m pip install --upgrade pip
-          pip install pytest pytest-cov obspy
+          pip install pytest pytest-cov obspy pyvista
           pip install .
       - name: test
         run: |

README.md

@@ -5,7 +5,7 @@
 
 [![Python Package using Conda](https://github.com/MIGG-NTU/PyTomoATT/actions/workflows/build-test-conda.yml/badge.svg?branch=devel)](https://github.com/MIGG-NTU/PyTomoATT/actions/workflows/build-test-conda.yml)
 [![Build documentations](https://github.com/MIGG-NTU/PyTomoATT/actions/workflows/build-docs.yml/badge.svg?branch=docs)](https://migg-ntu.github.io/PyTomoATT/)
-[![codecov](https://codecov.io/gh/MIGG-NTU/PyTomoATT/branch/devel/graph/badge.svg?token=EYOV0WOA2Y)](https://codecov.io/gh/MIGG-NTU/PyTomoATT)
+[![codecov](https://codecov.io/gh/TomoATT/PyTomoATT/graph/badge.svg?token=EYOV0WOA2Y)](https://codecov.io/gh/TomoATT/PyTomoATT)
 
 ![PyPI - License](https://img.shields.io/pypi/l/pytomoatt)
 ![PyPI](https://img.shields.io/pypi/v/pytomoatt)

pyproject.toml

@@ -17,9 +17,9 @@ dynamic = ["version"]
 
 classifiers = [
     "Programming Language :: Python",
-    "Programming Language :: Python :: 3.9",
     "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
+    "Programming Language :: Python :: 3.12",
 ]
 dependencies = [
     "numpy>=1.19.0",

pytomoatt/_version.py

@@ -1 +1 @@
-__version__ = '0.2.9'
+__version__ = '0.2.10'

pytomoatt/checkerboard.py

@@ -133,6 +133,20 @@ def copy(self):
         :rtype: Checker
         """
         return copy.deepcopy(self)
+
+    def to_attmodel(self):
+        """Convert to ATTModel object
+
+        :return: ATTModel object
+        :rtype: ATTModel
+        """
+        from .model import ATTModel
+        mod = ATTModel(self.para_fname)
+        mod.vel = self.vel_pert
+        mod.xi = self.xi
+        mod.eta = self.eta
+        mod.zeta = self.zeta
+        return mod
 
     def write(self, fname):
         """Write new model to h5 file

pytomoatt/model.py

@@ -5,7 +5,7 @@
 from .io.crustmodel import CrustModel
 from .io.asciimodel import ASCIIModel
 from .attarray import Dataset
-from .utils.common import init_axis
+from .utils.common import init_axis, km2deg
 import copy
 
 
@@ -73,14 +73,7 @@ def to_ani(self):
         """Convert to anisotropic strength (epsilon) and azimuth (phi)
         """
         self.epsilon = np.sqrt(self.eta**2+self.xi**2)
-        # self.phi = np.zeros_like(self.epsilon)
         self.phi = np.rad2deg(0.5*np.arctan2(self.eta, self.xi))
-        # idx = np.where(self.xi <= 0)
-        # self.phi[idx] = 90 + 0.5*atand(self.eta[idx]/self.xi[idx])
-        # idx = np.where((self.xi > 0) & (self.eta <= 0))
-        # self.phi[idx] = 180 + 0.5*atand(self.eta[idx]/self.xi[idx])
-        # idx = np.where((self.xi > 0) & (self.eta > 0))
-        # self.phi[idx] = 0.5*atand(self.eta[idx]/self.xi[idx])
 
     def to_xarray(self):
         """Convert to xarray
@@ -140,18 +133,61 @@ def grid_data_ascii(self, model_fname:str, **kwargs):
             self.n_rtp,
         )
 
-    def smooth(self, sigma=5.0):
+    def smooth(self, sigma=5.0, unit_deg=False, smooth_ani=False, **kwargs):
         """Gaussian smooth the 3D velocity model
 
-        :param sigma: Standard division of gaussian kernel in km, defaults to 10
-        :type sigma: scalar or sequence of scalars , optional
+        :param sigma: Standard deviation for Gaussian kernel.
+                      If scalar, apply to all dimensions.
+                      If sequence of 3, apply to [depth, lat, lon].
+                      Depth is always in km.
+                      Horizontal dimensions depend on unit_deg.
+        :type sigma: scalar or sequence of scalars
+        :param unit_deg: If True, horizontal sigma is in degrees. 
+                         If False, horizontal sigma is in km.
+                         Defaults to False.
+        :type unit_deg: bool
+        :param smooth_ani: If True, also smooth anisotropic parameters (xi, eta, zeta).
+                           Defaults to False.
+        :type smooth_ani: bool
+        :param kwargs: Additional arguments passed to scipy.ndimage.gaussian_filter
+
+        Example
+        -------------------
+        To smooth with 5 km in depth and 0.2 degrees in horizontal directions:
+        >>> model.smooth(sigma=[5.0, 0.2, 0.2], unit_deg=True)
+
+        To smooth with 5 km in depth and 20 km in horizontal directions:
+        >>> model.smooth(sigma=[5.0, 20.0, 20.0], unit_deg=False)
         """
-        if isinstance(sigma, (int, float)):
-            sigma_all = np.ones(3)*sigma/self.d_rtp/2/np.pi
-        elif len(sigma) == 3:
-            sigma_all = np.array(sigma)/self.d_rtp/2/np.pi
-        sigma_all[0:2] /= 111.19
-        self.vel = gaussian_filter(self.vel, sigma)
+        if np.isscalar(sigma):
+            sigma = [sigma, sigma, sigma]
+        elif len(sigma) != 3:
+            raise ValueError('sigma should be a scalar or a sequence of three scalars')
+
+        sigma = np.array(sigma, dtype=float)
+        sigma_pixel = np.zeros(3)
+
+        # Depth direction (always km)
+        sigma_pixel[0] = sigma[0] / self.d_rtp[0]
+
+        if unit_deg:
+            # Horizontal sigma is in degrees
+            sigma_pixel[1] = sigma[1] / self.d_rtp[1]
+            sigma_pixel[2] = sigma[2] / self.d_rtp[2]
+        else:
+            # Horizontal sigma is in km
+            # Latitude
+            sigma_pixel[1] = km2deg(sigma[1]) / self.d_rtp[1]            
+            # Longitude
+            mean_lat = np.mean(self.latitudes)
+            # 1 deg lon = cos(lat) * 1 deg lat
+            # so X km = km2deg(X) deg lat = km2deg(X) / cos(lat) deg lon
+            sigma_pixel[2] = km2deg(sigma[2]) / np.cos(np.deg2rad(mean_lat)) / self.d_rtp[2]
+        self.vel = gaussian_filter(self.vel, sigma_pixel, **kwargs)
+        if smooth_ani:
+            self.xi = gaussian_filter(self.xi, sigma_pixel, **kwargs)
+            self.eta = gaussian_filter(self.eta, sigma_pixel, **kwargs)
+            self.zeta = gaussian_filter(self.zeta, sigma_pixel, **kwargs)
 
     def calc_dv_avg(self):
         """calculate anomalies relative to average velocity at each depth

pytomoatt/utils/__init__.py

@@ -0,0 +1 @@
+_EARTH_RADIUS_KM = 6371.0

pytomoatt/utils/common.py

@@ -1,6 +1,6 @@
 import numpy as np
 from scipy.interpolate import griddata
-import pandas as pd
+from . import _EARTH_RADIUS_KM
 
 def sind(deg):
     rad = np.radians(deg)
@@ -37,6 +37,34 @@ def atand(x):
     return np.degrees(rad)
 
 
+def km2deg(km):
+    """ Convert km to degree
+
+    :param km: Distance in km
+    :type km: float
+    :return: Distance in degree
+    :rtype: float
+    """
+    circum = 2*np.pi*_EARTH_RADIUS_KM
+    conv = circum / 360
+    deg = km / conv
+    return deg
+
+
+def deg2km(deg):
+    """ Convert degree to km
+
+    :param deg: Distance in degree
+    :type deg: float
+    :return: Distance in km
+    :rtype: float
+    """
+    circum = 2*np.pi*_EARTH_RADIUS_KM
+    conv = circum / 360
+    km = deg * conv
+    return km
+
+
 def WGS84_to_cartesian(dep, lat, lon):
     """
     Convert WGS84 coordinates to cartesian coordinates

test/test_create_model.py

@@ -34,6 +34,8 @@ def test_checkerboard03(self):
             lim_y=[-0.5, 0.5],
             lim_z=[10, 120]
         )
+        mod = cm.to_attmodel()
+        mod.to_ani()
 
     def test_read_model(self):
         mod = ATTModel.read(self.out_fname, para_fname=self.para_fname)
@@ -46,6 +48,13 @@ def test_read_model(self):
             start_point=[mod.min_max_lon[0], mod.min_max_lat[1]],
             end_point=[mod.min_max_lon[1], mod.min_max_lat[1]], field='vel', flat_earth=True)
 
+    def test_smooth_model(self):
+        mod = ATTModel.read(self.out_fname, para_fname=self.para_fname)
+        mod_test1 = mod.copy()
+        mod_test1.smooth(sigma=[2.0, 0.1, 0.1], unit_deg=True, smooth_ani=False)
+        mod_test2 = mod.copy()
+        mod_test2.smooth(sigma=[5.0, 20, 20], smooth_ani=True)
+        mod_test1.write('smoothed_model1.h5')
 
 if __name__ == '__main__':
     test = TestATTModel()