Merge branch 'devel' of github.com:MIGG-NTU/PyTomoATT into devel

.github/workflows/build-docs.yml

@@ -11,10 +11,6 @@ on:
       config-path:
         required: true
         type: string
-    secrets:
-      token:
-        required: true
-
 jobs:
   deploy:
     runs-on: ubuntu-latest

.github/workflows/python-publish.yml

@@ -23,21 +23,24 @@ jobs:
     steps:
     - uses: actions/checkout@v3
     - name: Set up Python
-      uses: actions/setup-python@v3
+      uses: actions/setup-python@v5
       with:
         python-version: '3.x'
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
-        sudo pip install twine
+        pip install twine setuptools
     - name: Build package
       run: python setup.py sdist build
     - name: Publish package
       uses: pypa/gh-action-pypi-publish@27b31702a0e7fc50959f5ad993c78deac1bdfc29
       with:
         user: __token__
         password: ${{ secrets.PYPI_API_TOKEN }}
-
-  Build_docs:
-    uses: ./.github/workflows/build-docs.yml
 
+  Build_docs:
+    uses: MIGG-NTU/PyTomoATT/.github/workflows/build-docs.yml@devel
+    with:
+      config-path: .github/workflows/build-docs.yml
+    secrets: inherit
+

pytomoatt/attarray.py

@@ -11,6 +11,11 @@ def __init__(self, data_vars, coords, attrs=None) -> None:
         __slots__ = ()
         super().__init__(data_vars, coords, attrs)
 
+    @classmethod
+    def from_xarray(cls, dataset):
+        ds = cls(dataset.data_vars, dataset.coords)
+        return ds
+
     def interp_dep(self, depth:float, field:str):
         """Interpolate map view with given depth
 

pytomoatt/checkerboard.py

@@ -103,8 +103,8 @@ def checkerboard(self, period_x, period_y, period_z,
         self.dlnv = self.perturbation*pert_vel
         self.epsilon = np.abs(self.perturbation)*pert_ani
         self.phi = np.zeros_like(self.vel)
-        self.phi[np.where(self.perturbation>0)] = 120.
-        self.phi[np.where(self.perturbation<0)] = 60.
+        self.phi[np.where(self.perturbation>0)] = 135.
+        self.phi[np.where(self.perturbation<0)] = 45.
         self.xi = self.epsilon*cosd(2*self.phi)
         self.eta = self.epsilon*sind(2*self.phi)
 

pytomoatt/data.py

@@ -6,6 +6,15 @@
 
 
 class ATTData():
+    """Read data from HDF5 or ASCII file
+
+    :param fname: Path to data file
+    :type fname: str
+    :param fname_params: Path to input parameter file
+    :type fname_params: str
+    :param fname_grid: Path to grid file
+    :type fname_grid: str
+    """
     def __init__(self, fname:str,
                  fname_params:str,
                  fname_grid='OUTPUT_FILES/out_data_grid.h5'):
@@ -37,6 +46,23 @@ def read(cls, fname:str, fname_params:str,
                 fname_grid='OUTPUT_FILES/out_data_grid.h5',
                 group_name='model', dataset_name=None,
                 format='hdf5'):
+        """Read data from HDF5 or ASCII file
+
+        :param fname: Path to data file
+        :type fname: str
+        :param fname_params: Path to input parameter file
+        :type fname_params: str
+        :param fname_grid: Path to grid file
+        :type fname_grid: str
+        :param group_name: Name of the group in the HDF5 file
+        :type group_name: str
+        :param dataset_name: Name of the dataset in the HDF5 file
+        :type dataset_name: str
+        :param format: Format of the data file, defaults to 'hdf5'
+        :type format: str, optional
+        :return: An instance of ATTData
+        :rtype: ATTData
+        """
         attdata = cls(fname, fname_params, fname_grid)
         attdata.format = format
         # attdata.group_name = group_name

pytomoatt/distaz.py

@@ -3,6 +3,22 @@
 
 class DistAZ:
     """
+    DistAZ class
+
+    Calculate the distance, azimuth and back-azimuth between two points on the
+    Earth's surface.
+
+    :param lat1: Latitude of point 1
+    :type lat1: float
+    :param lon1: Longitude of point 1
+    :type lon1: float
+    :param lat2: Latitude of point 2
+    :type lat2: float
+    :param lon2: Longitude of point 2
+    :type lon2: float
+    :return: An instance of DistAZ
+    :rtype: DistAZ
+
      Subroutine to calculate the Great Circle Arc distance
         between two sets of geographic coordinates
     
@@ -44,25 +60,24 @@ def __init__(self, lat1, lon1, lat2, lon2):
         '''
 
         rad = 2. * np.pi / 360.0
-        """
-        c
-        c scolat and ecolat are the geocentric colatitudes
-        c as defined by Richter (pg. 318)
-        c
-        c Earth Flattening of 1/298.257 take from Bott (pg. 3)
-        c
-        """
+        '''
+        scolat and ecolat are the geocentric colatitudes
+        as defined by Richter (pg. 318)
+        
+        Earth Flattening of 1/298.257 take from Bott (pg. 3)
+        
+        '''
         sph = 1.0 / 298.257
 
         scolat = np.pi / 2.0 - np.arctan((1. - sph) * (1. - sph) * np.tan(lat1 * rad))
         ecolat = np.pi / 2.0 - np.arctan((1. - sph) * (1. - sph) * np.tan(lat2 * rad))
         slon = lon1 * rad
         elon = lon2 * rad
         """
-	c
-	c  a - e are as defined by Bullen (pg. 154, Sec 10.2)
-	c     These are defined for the pt. 1
-	c
+    
+	    a - e are as defined by Bullen (pg. 154, Sec 10.2)
+	    These are defined for the pt. 1
+
         """
         a = np.sin(scolat) * np.cos(slon)
         b = np.sin(scolat) * np.sin(slon)

pytomoatt/model.py

@@ -5,7 +5,7 @@
 from .io.crustmodel import CrustModel
 from .io.asciimodel import ASCIIModel
 from .attarray import Dataset
-from .utils import init_axis, acosd
+from .utils import init_axis, acosd, atand
 
 
 class ATTModel():
@@ -68,8 +68,12 @@ def to_ani(self):
         """
         self.epsilon = np.sqrt(self.eta**2+self.xi**2)
         self.phi = np.zeros_like(self.epsilon)
-        idx_non_zero = np.where(self.epsilon != 0)
-        self.phi[idx_non_zero] = 0.5*acosd(self.xi[idx_non_zero]/self.epsilon[idx_non_zero])
+        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
@@ -139,6 +143,14 @@ def smooth(self, sigma=5.0):
         sigma_all[0:2] /= 111.19
         self.vel = gaussian_filter(self.vel, sigma)
 
+    def calc_dv_avg(self):
+        """calculate anomalies relative to average velocity at each depth
+        """
+        self.dlnv = np.zeros_like(self.vel)
+        for i, _ in enumerate(self.depths):
+            avg = np.mean(self.vel[i, :, :])
+            self.dlnv[i, :, :] = 100 * (self.vel[i, :, :] - avg)/avg
+
     def calc_dv(self, ref_mod_fname: str):
         """calculate anomalies relative to another model