Merge branch 'master' of https://github.com/keurfonluu/FTeikPy

.github/sample.py

@@ -1,57 +1,51 @@
 import numpy
 import pyvista
-from scipy.ndimage import gaussian_filter
 
-from fteikpy import Eikonal2D
+from fteikpy import Eikonal2D, grid_to_meshio
 
 pyvista.set_plot_theme("document")
 
 
 def ray2line(ray):
     """Convert a ray array to PolyData."""
-    poly = pyvista.PolyData()
-
     nr = len(ray)
-    poly.points = numpy.column_stack((ray[:, 1], numpy.zeros(nr), -ray[:, 0])) * 1.0e-3
-    poly.lines = numpy.column_stack((numpy.full(nr - 1, 2), numpy.arange(nr - 1), numpy.arange(1, nr)))
+    points = numpy.column_stack((ray[:, 1], numpy.zeros(nr), -ray[:, 0]))
+    lines = numpy.column_stack((numpy.full(nr - 1, 2), numpy.arange(nr - 1), numpy.arange(1, nr)))
 
-    return poly
+    return pyvista.PolyData(points, lines=lines)
 
 
 # Import Marmousi velocity model
 vel = numpy.load("marmousi.npy")
-vel = gaussian_filter(vel, 5)
 
 # Calculate traveltime grid for one source point
-eik = Eikonal2D(vel, gridsize=(10.0, 10.0))
+eik = Eikonal2D(vel * 1.0e-3, gridsize=(0.01, 0.01))
+eik.smooth(0.05)
 tt = eik.solve((0.0, 0.0), nsweep=3, return_gradient=True)
 ttgrid = tt.grid.ravel()
 
 # Trace rays for 100 locations
 nrays = 100
 end_points = numpy.zeros((nrays, 2))
-end_points[:, 1] = numpy.linspace(4400.0, eik.xaxis[-1], nrays)
+end_points[:, 1] = numpy.linspace(4.4, eik.xaxis[-1], nrays)
 rays = tt.raytrace(end_points)
 trays = [tt(ray) for ray in rays]
 
 # Create mesh
-x, y, z = numpy.meshgrid(eik.xaxis * 1.0e-3, [0.0], -eik.zaxis * 1.0e-3)
-mesh = pyvista.StructuredGrid(x, y, z).cast_to_unstructured_grid()
-mesh["velocity"] = eik.grid.ravel() * 1.0e-3
-mesh["traveltime"] = ttgrid.copy()
+mesh = pyvista.from_meshio(grid_to_meshio(eik, tt))
 mesh2 = mesh.copy()
 
 # Create contour and lines
-contour = mesh.contour(isosurfaces=100, scalars="traveltime")
+contour = mesh.contour(isosurfaces=100, scalars="Traveltime")
 lines = [ray2line(ray) for ray in rays]
 
 # Initialize plotter
 p = pyvista.Plotter(window_size=(1500, 600), notebook=False)
 p.add_mesh(
     mesh,
-    scalars="velocity",
-    stitle="Velocity [km/s]",
+    scalars="Velocity",
     scalar_bar_args={
+        "title": "Velocity [km/s]",
         "height": 0.7,
         "width": 0.05,
         "position_x": 0.92,
@@ -61,7 +55,6 @@ def ray2line(ray):
         "title_font_size": 20,
         "label_font_size": 20,
         "font_family": "arial",
-        "shadow": True,
     },
 )
 p.add_mesh(
@@ -91,7 +84,6 @@ def ray2line(ray):
     zlabel="Elevation [km]",
     font_size=20,
     font_family="arial",
-    shadow=True,
 )
 p.show(
     cpos=[
@@ -112,9 +104,9 @@ def ray2line(ray):
     time.SetText(3, "Time: {:.2f} seconds".format(t))
 
     # Update isochrones
-    mesh2["traveltime"] = ttgrid.copy()
-    mesh2["traveltime"][ttgrid > t] = t
-    c = mesh2.contour(isosurfaces=int(n), scalars="traveltime")
+    mesh2["Traveltime"] = ttgrid.copy()
+    mesh2["Traveltime"][ttgrid > t] = t
+    c = mesh2.contour(isosurfaces=int(n), scalars="Traveltime")
     contour.points = c.points
     contour.lines = c.lines
 

CITATION.cff

@@ -5,5 +5,6 @@ authors:
   given-names: "Keurfon"
   orcid: "https://orcid.org/0000-0001-7927-0019"
 title: "fteikpy: Accurate Eikonal solver for Python"
+doi: 10.5281/zenodo.4269352
 url: https://github.com/keurfonluu/fteikpy
 license: BSD-3-Clause

README.rst

@@ -1,9 +1,9 @@
 fteikpy
 =======
 
-|License| |Stars| |Pyversions| |Version| |Downloads| |Code style: black| |Codacy Badge| |Codecov| |Build| |Travis|
+|License| |Stars| |Pyversions| |Version| |Downloads| |Code style: black| |Codacy Badge| |Codecov| |Build| |Travis| |DOI|
 
-**fteikpy** is a Python library that computes accurate first arrival traveltimes in 2D and 3D heterogenous isotropic velocity models. The algorithm handles properly the curvature of wavefronts close to the source which can be placed without any problem between grid points.
+**fteikpy** is a Python library that computes accurate first arrival traveltimes in 2D and 3D heterogeneous isotropic velocity models. The algorithm handles properly the curvature of wavefronts close to the source which can be placed without any problem between grid points.
 
 The code is based on `FTeik <https://github.com/Mark-Noble/FTeik-Eikonal-Solver>`__ implemented in Python and compiled `just-in-time <https://en.wikipedia.org/wiki/Just-in-time_compilation>`__ with `numba <https://numba.pydata.org/>`__.
 
@@ -54,7 +54,7 @@ Documentation
 
 Refer to the online `documentation <https://keurfonluu.github.io/fteikpy/>`__ for detailed description of the API and examples.
 
-Alternatively, the documentation can be built using `Sphinx <https://www.sphinx-doc.org/en/master/>`__
+Alternatively, the documentation can be built using `Sphinx <https://www.sphinx-doc.org/en/master/>`__:
 
 .. code:: bash
 
@@ -64,7 +64,7 @@ Alternatively, the documentation can be built using `Sphinx <https://www.sphinx-
 Usage
 -----
 
-The following example computes the traveltime grid in a 3D homogenous velocity model:
+The following example computes the traveltime grid in a 3D homogeneous velocity model:
 
 .. code-block:: python
 
@@ -115,6 +115,9 @@ Guidelines <https://github.com/keurfonluu/fteikpy/blob/master/CONTRIBUTING.rst>`
 .. |Codecov| image:: https://img.shields.io/codecov/c/github/keurfonluu/fteikpy.svg?style=flat
    :target: https://codecov.io/gh/keurfonluu/fteikpy
 
+.. |DOI| image:: https://zenodo.org/badge/DOI/10.5281/zenodo.4269352.svg?style=flat
+   :target: https://doi.org/10.5281/zenodo.4269352
+
 .. |Build| image:: https://img.shields.io/github/workflow/status/keurfonluu/fteikpy/Python%20package
    :target: https://github.com/keurfonluu/fteikpy
 

doc/source/api/index.rst

@@ -8,6 +8,7 @@ API Reference
    eikonal2d
    eikonal3d
    grids
+   io
 
 
 * :ref:`genindex`

doc/source/api/io.rst

@@ -0,0 +1,6 @@
+I/O
+===
+
+.. autofunction:: fteikpy.grid_to_meshio
+
+.. autofunction:: fteikpy.ray_to_meshio

fteikpy/__init__.py

@@ -1,6 +1,7 @@
 from .__about__ import __version__
 from ._grid import Grid2D, Grid3D, TraveltimeGrid2D, TraveltimeGrid3D
 from ._helpers import get_num_threads, set_num_threads
+from ._io import grid_to_meshio, ray_to_meshio
 from ._solver import Eikonal2D, Eikonal3D
 
 __all__ = [
@@ -12,5 +13,7 @@
     "TraveltimeGrid3D",
     "get_num_threads",
     "set_num_threads",
+    "grid_to_meshio",
+    "ray_to_meshio",
     "__version__",
 ]

fteikpy/_base.py

@@ -1,6 +1,8 @@
 from abc import ABC
 
 import numpy
+from scipy.interpolate import RegularGridInterpolator
+from scipy.ndimage import gaussian_filter
 
 from ._interp import interp2d, interp3d
 
@@ -49,6 +51,8 @@ def ndim(self):
 
 
 class BaseGrid2D(BaseGrid):
+    _ndim = 2
+
     def __call__(self, points, fill_value=numpy.nan):
         """
         Bilinear interpolation.
@@ -74,6 +78,50 @@ def __call__(self, points, fill_value=numpy.nan):
             fill_value,
         )
 
+    def resample(self, new_shape, method="linear"):
+        """
+        Resample grid.
+
+        Parameters
+        ----------
+        new_shape : array_like
+            New grid shape (nz, nx).
+        method : str ('linear' or 'nearest'), optional, default 'linear'
+            Interpolation method.
+
+        """
+        zaxis = self.zaxis
+        xaxis = self.xaxis
+        Z, X = numpy.meshgrid(
+            numpy.linspace(zaxis[0], zaxis[-1], new_shape[0]),
+            numpy.linspace(xaxis[0], xaxis[-1], new_shape[1]),
+            indexing="ij",
+        )
+
+        fn = RegularGridInterpolator(
+            points=(zaxis, xaxis), values=self._grid, method=method, bounds_error=False,
+        )
+        self._grid = fn([[z, x] for z, x in zip(Z.ravel(), X.ravel())]).reshape(
+            new_shape
+        )
+
+        self._gridsize = tuple(
+            a * b / c for a, b, c in zip(self.gridsize, self.shape, new_shape)
+        )
+
+    def smooth(self, sigma):
+        """
+        Smooth grid.
+
+        Parameters
+        ----------
+        sigma : scalar or array_like
+            Standard deviation in meters for Gaussian kernel.
+
+        """
+        sigma = numpy.full(2, sigma) if numpy.ndim(sigma) == 0 else numpy.asarray(sigma)
+        self._grid = gaussian_filter(self._grid, sigma / self._gridsize)
+
     @property
     def zaxis(self):
         """Return grid Z axis."""
@@ -86,6 +134,8 @@ def xaxis(self):
 
 
 class BaseGrid3D(BaseGrid):
+    _ndim = 3
+
     def __call__(self, points, fill_value=numpy.nan):
         """
         Trilinear interpolaton.
@@ -112,6 +162,55 @@ def __call__(self, points, fill_value=numpy.nan):
             fill_value,
         )
 
+    def resample(self, new_shape, method="linear"):
+        """
+        Resample grid.
+
+        Parameters
+        ----------
+        new_shape : array_like
+            New grid shape (nz, nx, ny).
+        method : str ('linear' or 'nearest'), optional, default 'linear'
+            Interpolation method.
+
+        """
+        zaxis = self.zaxis
+        xaxis = self.xaxis
+        yaxis = self.yaxis
+        Z, X, Y = numpy.meshgrid(
+            numpy.linspace(zaxis[0], zaxis[-1], new_shape[0]),
+            numpy.linspace(xaxis[0], xaxis[-1], new_shape[1]),
+            numpy.linspace(yaxis[0], yaxis[-1], new_shape[2]),
+            indexing="ij",
+        )
+
+        fn = RegularGridInterpolator(
+            points=(zaxis, xaxis, yaxis),
+            values=self._grid,
+            method=method,
+            bounds_error=False,
+        )
+        self._grid = fn(
+            [[z, x, y] for z, x, y in zip(Z.ravel(), X.ravel(), Y.ravel())]
+        ).reshape(new_shape)
+
+        self._gridsize = tuple(
+            a * b / c for a, b, c in zip(self.gridsize, self.shape, new_shape)
+        )
+
+    def smooth(self, sigma):
+        """
+        Smooth grid.
+
+        Parameters
+        ----------
+        sigma : scalar or array_like
+            Standard deviation in meters for Gaussian kernel.
+
+        """
+        sigma = numpy.full(3, sigma) if numpy.ndim(sigma) == 0 else numpy.asarray(sigma)
+        self._grid = gaussian_filter(self._grid, sigma / self._gridsize)
+
     @property
     def zaxis(self):
         """Return grid Z axis."""

fteikpy/_fteik/_common.py

@@ -3,16 +3,6 @@
 from .._common import jitted
 
 
-@jitted
-def first_index(x, y):
-    """Find index of first occurrence of x in array y."""
-    for i, v in enumerate(y):
-        if numpy.array_equal(v, x):
-            return i
-
-    return -1
-
-
 @jitted("f8(f8[:], f8[:], f8[:], f8[:])")
 def shrink(pcur, delta, lower, upper):
     """Stepsize shrinking factor if ray crosses interface."""