diff --git a/cookbook/gravmag_transform_deriv.py b/cookbook/gravmag_transform_deriv.py index 7591bd0f..1e1e9b64 100644 --- a/cookbook/gravmag_transform_deriv.py +++ b/cookbook/gravmag_transform_deriv.py @@ -33,40 +33,40 @@ mpl.subplot(2, 3, 1) mpl.title("x deriv (contour) + true (color map)") mpl.axis('scaled') -levels = mpl.contourf(xp, yp, gxz_true, shape, 12) +levels = mpl.contourf(yp, xp, gxz_true, shape, 12) mpl.colorbar(shrink=0.7) -mpl.contour(xp, yp, gxz, shape, 12, color='k') +mpl.contour(yp, xp, gxz, shape, 12, color='k') mpl.m2km() mpl.subplot(2, 3, 2) mpl.title("y deriv (contour) + true (color map)") mpl.axis('scaled') -levels = mpl.contourf(xp, yp, gyz_true, shape, 12) +levels = mpl.contourf(yp, xp, gyz_true, shape, 12) mpl.colorbar(shrink=0.7) -mpl.contour(xp, yp, gyz, shape, 12, color='k') +mpl.contour(yp, xp, gyz, shape, 12, color='k') mpl.m2km() mpl.subplot(2, 3, 3) mpl.title("z deriv (contour) + true (color map)") mpl.axis('scaled') -levels = mpl.contourf(xp, yp, gzz_true, shape, 8) +levels = mpl.contourf(yp, xp, gzz_true, shape, 8) mpl.colorbar(shrink=0.7) -mpl.contour(xp, yp, gzz, shape, levels, color='k') +mpl.contour(yp, xp, gzz, shape, levels, color='k') mpl.m2km() mpl.subplot(2, 3, 4) mpl.title("Difference x deriv") mpl.axis('scaled') -mpl.pcolor(xp, yp, (gxz_true - gxz), shape) +mpl.pcolor(yp, xp, (gxz_true - gxz), shape) mpl.colorbar(shrink=0.7) mpl.m2km() mpl.subplot(2, 3, 5) mpl.title("Difference y deriv") mpl.axis('scaled') -mpl.pcolor(xp, yp, (gyz_true - gyz), shape) +mpl.pcolor(yp, xp, (gyz_true - gyz), shape) mpl.colorbar(shrink=0.7) mpl.m2km() mpl.subplot(2, 3, 6) mpl.title("Difference z deriv") mpl.axis('scaled') -mpl.pcolor(xp, yp, (gzz_true - gzz), shape) +mpl.pcolor(yp, xp, (gzz_true - gzz), shape) mpl.colorbar(shrink=0.7) mpl.m2km() mpl.show() diff --git a/cookbook/grid_surfer.py b/cookbook/grid_surfer.py index 4a398e12..97530250 100644 --- a/cookbook/grid_surfer.py +++ b/cookbook/grid_surfer.py @@ -8,8 +8,8 @@ # Will download the archive and save it with the default name archive = datasets.fetch_bouguer_alps_egm() -# Load the GRD file and convert in three numpy-arrays (y, x, bouguer) -y, x, bouguer, shape = gridder.load_surfer(archive, fmt='ascii') +# Load the GRD file and convert in three numpy-arrays (x, y, bouguer) +x, y, bouguer, shape = gridder.load_surfer(archive, fmt='ascii') mpl.figure() mpl.axis('scaled') diff --git a/doc/api/gravmag.transform.rst b/doc/api/gravmag.transform.rst index 4a8da597..faebd78d 100644 --- a/doc/api/gravmag.transform.rst +++ b/doc/api/gravmag.transform.rst @@ -1,6 +1,6 @@ .. _fatiando_gravmag_transform: -Space domain processing (``fatiando.gravmag.transform``) +Potential field transformations (``fatiando.gravmag.transform``) ============================================================================ .. automodule:: fatiando.gravmag.transform diff --git a/doc/changelog.rst b/doc/changelog.rst index a2352000..324649ce 100644 --- a/doc/changelog.rst +++ b/doc/changelog.rst @@ -10,6 +10,20 @@ Version (development) **Changes**: +* **IMPORTANT BUG FIX**: ``fatiando,gridder.regular`` and many other places in + Fatiando where using the wrong convention for x, y dimensions. + x should point North and y East. Thus, a data matrix (regular grid) should + have x varying in the lines and y varying in the columns. This is **oposite** + what we had. This fix also changes the ``shape`` argument to be ``(nx, ny)`` + instead of ``(ny, nx)``. **Users should be aware of this and double check + their code.** + (`PR 196 `__) +* More stable derivatives in ``fatiando.gravamag.transform``. The horizontal + derivatives default to central finite-differences for greater stability. The + FFT based derivatives use a grid padding to avoid edge effects. + Thanks to `Matteo Niccoli `__ for suggesting + this fix. + (`PR 196 `__) * **Renamed** ``fatiando.gravmag.fourier.ansig`` to ``fatiando.gravmag.transform.tga`` (`PR 186 `__) diff --git a/fatiando/gravmag/__init__.py b/fatiando/gravmag/__init__.py index 181a727c..c26fdd62 100644 --- a/fatiando/gravmag/__init__.py +++ b/fatiando/gravmag/__init__.py @@ -43,8 +43,8 @@ * :mod:`~fatiando.gravmag.normal_gravity`: Compute normal gravity and reductions. * :mod:`~fatiando.gravmag.eqlayer`: Equivalent layer processing -* :mod:`~fatiando.gravmag.transform`: Space domain potential field - transformations, like upward continuation +* :mod:`~fatiando.gravmag.transform`: Potential field transformations, + like upward continuation, derivatives, etc * :mod:`~fatiando.gravmag.imaging`: Imaging methods for potential fields for estimating physical property distributions * :mod:`~fatiando.gravmag.tensor`: Utilities for operating on the gradient diff --git a/fatiando/gravmag/sphere.py b/fatiando/gravmag/sphere.py index bd54d449..8661bd81 100644 --- a/fatiando/gravmag/sphere.py +++ b/fatiando/gravmag/sphere.py @@ -136,42 +136,6 @@ def tf(xp, yp, zp, spheres, inc, dec, pmag=None): * tf : array The total-field anomaly - Example: - - >>> from fatiando import mesher, gridder, utils - >>> # Set the inclination and declination of the regional field - >>> inc, dec = -30, 45 - >>> # Create a sphere model - >>> model = [ - ... # One with induced magnetization - ... mesher.Sphere(1000, 1000, 600, 500, {'magnetization':5}), - ... # and one with remanent - ... mesher.Sphere(-1000, -1000, 600, 500, - ... {'magnetization':utils.ang2vec(10, 70, -5)})] - >>> # Create a regular grid at 100m height - >>> shape = (4, 4) - >>> area = (-3000, 3000, -3000, 3000) - >>> xp, yp, zp = gridder.regular(area, shape, z=-100) - >>> # Calculate the anomaly for a given regional field - >>> for t in tf(xp, yp, zp, model, inc, dec): - ... print '%15.8e' % t - 2.72951375e+01 - 3.63637351e+01 - 5.35842876e+00 - 6.65189557e-02 - 6.01998831e+01 - -1.71920499e+03 - 3.30025228e+00 - -5.13176612e+00 - 1.47871812e+00 - -3.96103758e+01 - -1.85654021e+02 - 2.18002960e+01 - -2.82713826e+00 - -1.10341542e+01 - 1.93353982e+01 - 2.03174254e+01 - """ if xp.shape != yp.shape != zp.shape: raise ValueError("Input arrays xp, yp, and zp must have same shape!") @@ -224,40 +188,6 @@ def bx(xp, yp, zp, spheres): * bx: array The x component of the magnetic induction - Example: - - >>> from fatiando import mesher, gridder, utils - >>> # Create a model formed by two spheres - >>> # The magnetization of each sphere is a vector - >>> model = [ - ... mesher.Sphere(1000, 1000, 600, 500, - ... {'magnetization':utils.ang2vec(13, -10, 28)}), - ... mesher.Sphere(-1000, -1000, 600, 500, - ... {'magnetization':utils.ang2vec(10, 70, -5)})] - >>> # Create a regular grid at 100m height - >>> shape = (4, 4) - >>> area = (-3000, 3000, -3000, 3000) - >>> xp, yp, zp = gridder.regular(area, shape, z=-100) - >>> # Calculate the bx component - >>> for b in bx(xp, yp, zp, model): - ... print '%15.8e' % b - 1.58002397e+01 - -1.76820799e+01 - -1.48049248e+01 - -5.75238567e+00 - 9.17572697e+01 - -4.94607307e+02 - -7.92213872e+01 - -4.37781621e+00 - 2.97032297e+01 - 7.36803996e+01 - -1.73332620e+03 - 1.15884125e+02 - 4.55847152e+00 - -1.31173236e+01 - -6.42912671e+01 - 2.98847909e+01 - """ if xp.shape != yp.shape != zp.shape: raise ValueError("Input arrays xp, yp, and zp must have same shape!") @@ -294,40 +224,6 @@ def by(xp, yp, zp, spheres): * by: array The y component of the magnetic induction - Example: - - >>> from fatiando import mesher, gridder, utils - >>> # Create a model formed by two spheres - >>> # The magnetization of each sphere is a vector - >>> model = [ - ... mesher.Sphere(1000, 1000, 600, 500, - ... {'magnetization':utils.ang2vec(13, -10, 28)}), - ... mesher.Sphere(-1000, -1000, 600, 500, - ... {'magnetization':utils.ang2vec(10, 70, -5)})] - >>> # Create a regular grid at 100m height - >>> shape = (4, 4) - >>> area = (-3000, 3000, -3000, 3000) - >>> xp, yp, zp = gridder.regular(area, shape, z=-100) - >>> # Calculate the by component - >>> for b in by(xp, yp, zp, model): - ... print '%15.8e' % b - 2.51394441e+01 - 5.71383698e+01 - 7.46666729e+00 - -4.53730551e+00 - 7.44792258e+00 - 8.22174414e+01 - 4.53451310e+01 - -3.06885735e+01 - -2.49929765e+01 - -8.41961087e+01 - -9.17412395e+02 - -3.18422413e+01 - -1.32728556e+01 - -3.03825859e+01 - 6.67990083e+01 - 4.21366247e+01 - """ if xp.shape != yp.shape != zp.shape: raise ValueError("Input arrays xp, yp, and zp must have same shape!") @@ -364,40 +260,6 @@ def bz(xp, yp, zp, spheres): * bz: array The z component of the magnetic induction - Example: - - >>> from fatiando import mesher, gridder, utils - >>> # Create a model formed by two spheres - >>> # The magnetization of each sphere is a vector - >>> model = [ - ... mesher.Sphere(1000, 1000, 600, 500, - ... {'magnetization':utils.ang2vec(13, -10, 28)}), - ... mesher.Sphere(-1000, -1000, 600, 500, - ... {'magnetization':utils.ang2vec(10, 70, -5)})] - >>> # Create a regular grid at 100m height - >>> shape = (4, 4) - >>> area = (-3000, 3000, -3000, 3000) - >>> xp, yp, zp = gridder.regular(area, shape, z=-100) - >>> # Calculate the bz component - >>> for b in bz(xp, yp, zp, model): - ... print '%15.8e' % b - -1.13152279e+01 - -3.24362266e+01 - -1.63235805e+01 - -4.48136597e+00 - -1.27492012e+01 - 2.89101261e+03 - -1.30263918e+01 - -9.64182996e+00 - -6.45566985e+00 - 3.32987598e+01 - -7.08905624e+02 - -5.55139945e+01 - -1.35745203e+00 - 2.91949888e+00 - -2.78345635e+01 - -1.69425703e+01 - """ if xp.shape != yp.shape != zp.shape: raise ValueError("Input arrays xp, yp, and zp must have same shape!") @@ -481,34 +343,6 @@ def gxx(xp, yp, zp, spheres, dens=None): * res : array The field calculated on xp, yp, zp - Example: - - >>> from fatiando import mesher, gridder - >>> # Create a sphere model - >>> sphere = [mesher.Sphere(0, 0, 5, 1, {'density':1.})] - >>> # Create a regular grid at 0m height - >>> shape = (4, 4) - >>> area = (-30, 30, -30, 30) - >>> xp, yp, zp = gridder.regular(area, shape, z=0) - >>> for g in gxx(xp, yp, zp, sphere): - ... print '%15.8e' % g - 1.71893959e-06 - -6.02473678e-06 - -6.02473678e-06 - 1.71893959e-06 - 1.39192195e-05 - 2.76067131e-05 - 2.76067131e-05 - 1.39192195e-05 - 1.39192195e-05 - 2.76067131e-05 - 2.76067131e-05 - 1.39192195e-05 - 1.71893959e-06 - -6.02473678e-06 - -6.02473678e-06 - 1.71893959e-06 - """ if xp.shape != yp.shape != zp.shape: raise ValueError("Input arrays xp, yp, and zp must have same shape!") @@ -551,35 +385,6 @@ def gxy(xp, yp, zp, spheres, dens=None): * res : array The field calculated on xp, yp, zp - Example: - - >>> from fatiando import mesher, gridder - >>> # Create a sphere model - >>> sphere = [mesher.Sphere(0, 0, 5, 1, {'density':1.})] - >>> # Create a regular grid at 0m height - >>> shape = (4, 4) - >>> area = (-30, 30, -30, 30) - >>> xp, yp, zp = gridder.regular(area, shape, z=0) - >>> # Calculate the gxy component - >>> for g in gxy(xp, yp, zp, sphere): - ... print '%15.8e' % g - 5.30415646e-06 - 7.47898359e-06 - -7.47898359e-06 - -5.30415646e-06 - 7.47898359e-06 - 1.10426852e-04 - -1.10426852e-04 - -7.47898359e-06 - -7.47898359e-06 - -1.10426852e-04 - 1.10426852e-04 - 7.47898359e-06 - -5.30415646e-06 - -7.47898359e-06 - 7.47898359e-06 - 5.30415646e-06 - """ if xp.shape != yp.shape != zp.shape: raise ValueError("Input arrays xp, yp, and zp must have same shape!") @@ -622,35 +427,6 @@ def gxz(xp, yp, zp, spheres, dens=None): * res : array The field calculated on xp, yp, zp - Example: - - >>> from fatiando import mesher, gridder - >>> # Create a sphere model - >>> sphere = [mesher.Sphere(0, 0, 5, 1, {'density':1.})] - >>> # Create a regular grid at 0m height - >>> shape = (4, 4) - >>> area = (-30, 30, -30, 30) - >>> xp, yp, zp = gridder.regular(area, shape, z=0) - >>> # Calculate the gxz component - >>> for g in gxz(xp, yp, zp, sphere): - ... print '%15.8e' % g - 8.84026077e-07 - 1.24649726e-06 - -1.24649726e-06 - -8.84026077e-07 - 3.73949179e-06 - 5.52134262e-05 - -5.52134262e-05 - -3.73949179e-06 - 3.73949179e-06 - 5.52134262e-05 - -5.52134262e-05 - -3.73949179e-06 - 8.84026077e-07 - 1.24649726e-06 - -1.24649726e-06 - -8.84026077e-07 - """ if xp.shape != yp.shape != zp.shape: raise ValueError("Input arrays xp, yp, and zp must have same shape!") @@ -693,35 +469,6 @@ def gyy(xp, yp, zp, spheres, dens=None): * res : array The field calculated on xp, yp, zp - Example: - - >>> from fatiando import mesher, gridder - >>> # Create a sphere model - >>> sphere = [mesher.Sphere(0, 0, 5, 1, {'density':1.})] - >>> # Create a regular grid at 0m height - >>> shape = (4, 4) - >>> area = (-30, 30, -30, 30) - >>> xp, yp, zp = gridder.regular(area, shape, z=0) - >>> # Calculate the gyy component - >>> for g in gyy(xp, yp, zp, sphere): - ... print '%15.8e' % g - 1.71893959e-06 - 1.39192195e-05 - 1.39192195e-05 - 1.71893959e-06 - -6.02473678e-06 - 2.76067131e-05 - 2.76067131e-05 - -6.02473678e-06 - -6.02473678e-06 - 2.76067131e-05 - 2.76067131e-05 - -6.02473678e-06 - 1.71893959e-06 - 1.39192195e-05 - 1.39192195e-05 - 1.71893959e-06 - """ if xp.shape != yp.shape != zp.shape: raise ValueError("Input arrays xp, yp, and zp must have same shape!") @@ -764,35 +511,6 @@ def gyz(xp, yp, zp, spheres, dens=None): * res : array The field calculated on xp, yp, zp - Example: - - >>> from fatiando import mesher, gridder - >>> # Create a sphere model - >>> sphere = [mesher.Sphere(0, 0, 5, 1, {'density':1.})] - >>> # Create a regular grid at 0m height - >>> shape = (4, 4) - >>> area = (-30, 30, -30, 30) - >>> xp, yp, zp = gridder.regular(area, shape, z=0) - >>> # Calculate the gyz component - >>> for g in gyz(xp, yp, zp, sphere): - ... print '%15.8e' % g - 8.84026077e-07 - 3.73949179e-06 - 3.73949179e-06 - 8.84026077e-07 - 1.24649726e-06 - 5.52134262e-05 - 5.52134262e-05 - 1.24649726e-06 - -1.24649726e-06 - -5.52134262e-05 - -5.52134262e-05 - -1.24649726e-06 - -8.84026077e-07 - -3.73949179e-06 - -3.73949179e-06 - -8.84026077e-07 - """ if xp.shape != yp.shape != zp.shape: raise ValueError("Input arrays xp, yp, and zp must have same shape!") @@ -835,35 +553,6 @@ def gzz(xp, yp, zp, spheres, dens=None): * res : array The field calculated on xp, yp, zp - Example: - - >>> from fatiando import mesher, gridder - >>> # Create a sphere model - >>> sphere = [mesher.Sphere(0, 0, 5, 1, {'density':1.})] - >>> # Create a regular grid at 0m height - >>> shape = (4, 4) - >>> area = (-30, 30, -30, 30) - >>> xp, yp, zp = gridder.regular(area, shape, z=0) - >>> # Calculate the gzz component - >>> for g in gzz(xp, yp, zp, sphere): - ... print '%15.8e' % g - -3.43787919e-06 - -7.89448267e-06 - -7.89448267e-06 - -3.43787919e-06 - -7.89448267e-06 - -5.52134262e-05 - -5.52134262e-05 - -7.89448267e-06 - -7.89448267e-06 - -5.52134262e-05 - -5.52134262e-05 - -7.89448267e-06 - -3.43787919e-06 - -7.89448267e-06 - -7.89448267e-06 - -3.43787919e-06 - """ if xp.shape != yp.shape != zp.shape: raise ValueError("Input arrays xp, yp, and zp must have same shape!") @@ -918,37 +607,6 @@ def kernelxx(xp, yp, zp, sphere): * res : array The function calculated on xp, yp, zp - - Example: - - >>> from fatiando import mesher, gridder - >>> # Create a sphere model - >>> sphere = mesher.Sphere(0, 0, 5, 1) - >>> # Create a regular grid at 0m height - >>> shape = (4, 4) - >>> area = (-30, 30, -30, 30) - >>> xp, yp, zp = gridder.regular(area, shape, z=0) - >>> # Calculate the function - >>> kxx = kernelxx(xp, yp, zp, sphere) - >>> for k in kxx: - ... print '%15.8e' % k - 2.57596223e-05 - -9.02852807e-05 - -9.02852807e-05 - 2.57596223e-05 - 2.08590131e-04 - 4.13707675e-04 - 4.13707675e-04 - 2.08590131e-04 - 2.08590131e-04 - 4.13707675e-04 - 4.13707675e-04 - 2.08590131e-04 - 2.57596223e-05 - -9.02852807e-05 - -9.02852807e-05 - 2.57596223e-05 - """ if xp.shape != yp.shape != zp.shape: raise ValueError("Input arrays xp, yp, and zp must have same shape!") @@ -995,37 +653,6 @@ def kernelxy(xp, yp, zp, sphere): * res : array The function calculated on xp, yp, zp - - Example: - - >>> from fatiando import mesher, gridder - >>> # Create a sphere model - >>> sphere = mesher.Sphere(0, 0, 5, 1) - >>> # Create a regular grid at 0m height - >>> shape = (4, 4) - >>> area = (-30, 30, -30, 30) - >>> xp, yp, zp = gridder.regular(area, shape, z=0) - >>> # Calculate the function - >>> kxy = kernelxy(xp, yp, zp, sphere) - >>> for k in kxy: - ... print '%15.8e' % k - 7.94868344e-05 - 1.12078279e-04 - -1.12078279e-04 - -7.94868344e-05 - 1.12078279e-04 - 1.65483070e-03 - -1.65483070e-03 - -1.12078279e-04 - -1.12078279e-04 - -1.65483070e-03 - 1.65483070e-03 - 1.12078279e-04 - -7.94868344e-05 - -1.12078279e-04 - 1.12078279e-04 - 7.94868344e-05 - """ if xp.shape != yp.shape != zp.shape: raise ValueError("Input arrays xp, yp, and zp must have same shape!") @@ -1072,37 +699,6 @@ def kernelxz(xp, yp, zp, sphere): * res : array The function calculated on xp, yp, zp - - Example: - - >>> from fatiando import mesher, gridder - >>> # Create a sphere model - >>> sphere = mesher.Sphere(0, 0, 5, 1) - >>> # Create a regular grid at 0m height - >>> shape = (4, 4) - >>> area = (-30, 30, -30, 30) - >>> xp, yp, zp = gridder.regular(area, shape, z=0) - >>> # Calculate the function - >>> kxz = kernelxz(xp, yp, zp, sphere) - >>> for k in kxz: - ... print '%15.8e' % k - 1.32478057e-05 - 1.86797132e-05 - -1.86797132e-05 - -1.32478057e-05 - 5.60391397e-05 - 8.27415349e-04 - -8.27415349e-04 - -5.60391397e-05 - 5.60391397e-05 - 8.27415349e-04 - -8.27415349e-04 - -5.60391397e-05 - 1.32478057e-05 - 1.86797132e-05 - -1.86797132e-05 - -1.32478057e-05 - """ if xp.shape != yp.shape != zp.shape: raise ValueError("Input arrays xp, yp, and zp must have same shape!") @@ -1149,37 +745,6 @@ def kernelyy(xp, yp, zp, sphere): * res : array The function calculated on xp, yp, zp - - Example: - - >>> from fatiando import mesher, gridder - >>> # Create a sphere model - >>> sphere = mesher.Sphere(0, 0, 5, 1) - >>> # Create a regular grid at 0m height - >>> shape = (4, 4) - >>> area = (-30, 30, -30, 30) - >>> xp, yp, zp = gridder.regular(area, shape, z=0) - >>> # Calculate the function - >>> kyy = kernelyy(xp, yp, zp, sphere) - >>> for k in kyy: - ... print '%15.8e' % k - 2.57596223e-05 - 2.08590131e-04 - 2.08590131e-04 - 2.57596223e-05 - -9.02852807e-05 - 4.13707675e-04 - 4.13707675e-04 - -9.02852807e-05 - -9.02852807e-05 - 4.13707675e-04 - 4.13707675e-04 - -9.02852807e-05 - 2.57596223e-05 - 2.08590131e-04 - 2.08590131e-04 - 2.57596223e-05 - """ if xp.shape != yp.shape != zp.shape: raise ValueError("Input arrays xp, yp, and zp must have same shape!") @@ -1226,37 +791,6 @@ def kernelyz(xp, yp, zp, sphere): * res : array The function calculated on xp, yp, zp - - Example: - - >>> from fatiando import mesher, gridder - >>> # Create a sphere model - >>> sphere = mesher.Sphere(0, 0, 5, 1) - >>> # Create a regular grid at 0m height - >>> shape = (4, 4) - >>> area = (-30, 30, -30, 30) - >>> xp, yp, zp = gridder.regular(area, shape, z=0) - >>> # Calculate the function - >>> kyz = kernelyz(xp, yp, zp, sphere) - >>> for k in kyz: - ... print '%15.8e' % k - 1.32478057e-05 - 5.60391397e-05 - 5.60391397e-05 - 1.32478057e-05 - 1.86797132e-05 - 8.27415349e-04 - 8.27415349e-04 - 1.86797132e-05 - -1.86797132e-05 - -8.27415349e-04 - -8.27415349e-04 - -1.86797132e-05 - -1.32478057e-05 - -5.60391397e-05 - -5.60391397e-05 - -1.32478057e-05 - """ if xp.shape != yp.shape != zp.shape: raise ValueError("Input arrays xp, yp, and zp must have same shape!") @@ -1303,37 +837,6 @@ def kernelzz(xp, yp, zp, sphere): * res : array The function calculated on xp, yp, zp - - Example: - - >>> from fatiando import mesher, gridder - >>> # Create a sphere model - >>> sphere = mesher.Sphere(0, 0, 5, 1) - >>> # Create a regular grid at 0m height - >>> shape = (4, 4) - >>> area = (-30, 30, -30, 30) - >>> xp, yp, zp = gridder.regular(area, shape, z=0) - >>> # Calculate the function - >>> kzz = kernelzz(xp, yp, zp, sphere) - >>> for k in kzz: - ... print '%15.8e' % k - -5.15192445e-05 - -1.18304851e-04 - -1.18304851e-04 - -5.15192445e-05 - -1.18304851e-04 - -8.27415349e-04 - -8.27415349e-04 - -1.18304851e-04 - -1.18304851e-04 - -8.27415349e-04 - -8.27415349e-04 - -1.18304851e-04 - -5.15192445e-05 - -1.18304851e-04 - -1.18304851e-04 - -5.15192445e-05 - """ if xp.shape != yp.shape != zp.shape: raise ValueError("Input arrays xp, yp, and zp must have same shape!") diff --git a/fatiando/gravmag/transform.py b/fatiando/gravmag/transform.py index b88a45c3..02265171 100644 --- a/fatiando/gravmag/transform.py +++ b/fatiando/gravmag/transform.py @@ -1,26 +1,28 @@ """ -Space domain potential field transformations, like upward continuation, -derivatives and total mass. +Potential field transformations, like upward continuation and derivatives. + +.. note:: Most, if not all, functions here required gridded data. **Transformations** -* :func:`~fatiando.gravmag.transform.upcontinue`: Upward continuation of the - vertical component of gravity :math:`g_z` using numerical integration +* :func:`~fatiando.gravmag.transform.upcontinue`: Upward continuation of + gridded potential field data on a level surface * :func:`~fatiando.gravmag.transform.tga`: Calculate the amplitude of the total gradient (also called the analytic signal) **Derivatives** * :func:`~fatiando.gravmag.transform.derivx`: Calculate the n-th order - derivative of a potential field in the x-direction + derivative of a potential field in the x-direction (North-South) * :func:`~fatiando.gravmag.transform.derivy`: Calculate the n-th order - derivative of a potential field in the y-direction + derivative of a potential field in the y-direction (East-West) * :func:`~fatiando.gravmag.transform.derivz`: Calculate the n-th order derivative of a potential field in the z-direction ---- """ +from __future__ import division import numpy @@ -75,7 +77,7 @@ def upcontinue(gz, height, xp, yp, dims): return gzcont -def tga(x, y, data, shape): +def tga(x, y, data, shape, method='fd'): """ Calculate the total gradient amplitude. @@ -96,8 +98,12 @@ def tga(x, y, data, shape): The x and y coordinates of the grid points * data : 1D-array The potential field at the grid points - * shape : tuple = (ny, nx) + * shape : tuple = (nx, ny) The shape of the grid + * method : string + The method used to calculate the horizontal derivatives. Options are: + ``'fd'`` for finite-difference (more stable) or ``'fft'`` for the Fast + Fourier Transform. The z derivative is always calculated by FFT. Returns: @@ -105,14 +111,14 @@ def tga(x, y, data, shape): The amplitude of the total gradient """ - dx = derivx(x, y, data, shape) - dy = derivy(x, y, data, shape) + dx = derivx(x, y, data, shape, method=method) + dy = derivy(x, y, data, shape, method=method) dz = derivz(x, y, data, shape) res = numpy.sqrt(dx ** 2 + dy ** 2 + dz ** 2) return res -def derivx(x, y, data, shape, order=1): +def derivx(x, y, data, shape, order=1, method='fd'): """ Calculate the derivative of a potential field in the x direction. @@ -130,10 +136,14 @@ def derivx(x, y, data, shape, order=1): The x and y coordinates of the grid points * data : 1D-array The potential field at the grid points - * shape : tuple = (ny, nx) + * shape : tuple = (nx, ny) The shape of the grid * order : int The order of the derivative + * method : string + The method used to calculate the derivatives. Options are: + ``'fd'`` for central finite-differences (more stable) or ``'fft'`` + for the Fast Fourier Transform. Returns: @@ -141,13 +151,30 @@ def derivx(x, y, data, shape, order=1): The derivative """ - Fx = _getfreqs(x, y, data, shape)[0].astype('complex') - # Multiply by 1j because I don't multiply it in _deriv (this way _deriv can - # be used for the z derivative as well) - return _deriv(Fx * 1j, data, shape, order) - - -def derivy(x, y, data, shape, order=1): + nx, ny = shape + assert method in ['fft', 'fd'], \ + 'Invalid method "{}".'.format(method) + if method == 'fft': + # Pad the array with the edge values to avoid instability + padded, padx, pady = _pad_data(data, shape) + kx, _ = _fftfreqs(x, y, shape, padded.shape) + deriv_ft = numpy.fft.fft2(padded)*(kx*1j)**order + deriv_pad = numpy.real(numpy.fft.ifft2(deriv_ft)) + # Remove padding from derivative + deriv = deriv_pad[padx : padx + nx, pady : pady + ny] + elif method == 'fd': + datamat = data.reshape(shape) + dx = (x.max() - x.min())/(nx - 1) + deriv = numpy.empty_like(datamat) + deriv[1:-1, :] = (datamat[2:, :] - datamat[:-2, :])/(2*dx) + deriv[0, :] = deriv[1, :] + deriv[-1, :] = deriv[-2, :] + if order > 1: + deriv = derivx(x, y, deriv, shape, order=order - 1, method='fd') + return deriv.ravel() + + +def derivy(x, y, data, shape, order=1, method='fd'): """ Calculate the derivative of a potential field in the y direction. @@ -165,10 +192,14 @@ def derivy(x, y, data, shape, order=1): The x and y coordinates of the grid points * data : 1D-array The potential field at the grid points - * shape : tuple = (ny, nx) + * shape : tuple = (nx, ny) The shape of the grid * order : int The order of the derivative + * method : string + The method used to calculate the derivatives. Options are: + ``'fd'`` for central finite-differences (more stable) or ``'fft'`` + for the Fast Fourier Transform. Returns: @@ -176,13 +207,30 @@ def derivy(x, y, data, shape, order=1): The derivative """ - Fy = _getfreqs(x, y, data, shape)[1].astype('complex') - # Multiply by 1j because I don't multiply it in _deriv (this way _deriv can - # be used for the z derivative as well) - return _deriv(Fy * 1j, data, shape, order) - - -def derivz(x, y, data, shape, order=1): + nx, ny = shape + assert method in ['fft', 'fd'], \ + 'Invalid method "{}".'.format(method) + if method == 'fft': + # Pad the array with the edge values to avoid instability + padded, padx, pady = _pad_data(data, shape) + _, ky = _fftfreqs(x, y, shape, padded.shape) + deriv_ft = numpy.fft.fft2(padded)*(ky*1j)**order + deriv_pad = numpy.real(numpy.fft.ifft2(deriv_ft)) + # Remove padding from derivative + deriv = deriv_pad[padx : padx + nx, pady : pady + ny] + elif method == 'fd': + datamat = data.reshape(shape) + dy = (y.max() - y.min())/(ny - 1) + deriv = numpy.empty_like(datamat) + deriv[:, 1:-1] = (datamat[:, 2:] - datamat[:, :-2])/(2*dy) + deriv[:, 0] = deriv[:, 1] + deriv[:, -1] = deriv[:, -2] + if order > 1: + deriv = derivy(x, y, deriv, shape, order=order - 1, method='fd') + return deriv.ravel() + + +def derivz(x, y, data, shape, order=1, method='fft'): """ Calculate the derivative of a potential field in the z direction. @@ -200,10 +248,13 @@ def derivz(x, y, data, shape, order=1): The x and y coordinates of the grid points * data : 1D-array The potential field at the grid points - * shape : tuple = (ny, nx) + * shape : tuple = (nx, ny) The shape of the grid * order : int The order of the derivative + * method : string + The method used to calculate the derivatives. Options are: + ``'fft'`` for the Fast Fourier Transform. Returns: @@ -211,27 +262,40 @@ def derivz(x, y, data, shape, order=1): The derivative """ - Fx, Fy = _getfreqs(x, y, data, shape) - freqs = numpy.sqrt(Fx ** 2 + Fy ** 2) - return _deriv(freqs, data, shape, order) - - -def _getfreqs(x, y, data, shape): + assert method == 'fft', \ + "Invalid method '{}'".format(method) + nx, ny = shape + # Pad the array with the edge values to avoid instability + padded, padx, pady = _pad_data(data, shape) + kx, ky = _fftfreqs(x, y, shape, padded.shape) + deriv_ft = numpy.fft.fft2(padded)*numpy.sqrt(kx**2 + ky**2)**order + deriv = numpy.real(numpy.fft.ifft2(deriv_ft)) + # Remove padding from derivative + return deriv[padx : padx + nx, pady : pady + ny].ravel() + + +def _pad_data(data, shape): + n = _nextpow2(numpy.max(shape)) + nx, ny = shape + padx = (n - nx)//2 + pady = (n - ny)//2 + padded = numpy.pad(data.reshape(shape), ((padx, padx), (pady, pady)), + mode='edge') + return padded, padx, pady + + +def _nextpow2(i): + buf = numpy.ceil(numpy.log(i)/numpy.log(2)) + return int(2**buf) + + +def _fftfreqs(x, y, shape, padshape): """ Get two 2D-arrays with the wave numbers in the x and y directions. """ - ny, nx = shape - dx = float(x.max() - x.min()) / float(nx - 1) - fx = numpy.fft.fftfreq(nx, dx) - dy = float(y.max() - y.min()) / float(ny - 1) - fy = numpy.fft.fftfreq(ny, dy) - return numpy.meshgrid(fx, fy) - - -def _deriv(freqs, data, shape, order): - """ - Calculate a generic derivative using the FFT. - """ - fgrid = (2. * numpy.pi) * numpy.fft.fft2(numpy.reshape(data, shape)) - deriv = numpy.real(numpy.fft.ifft2((freqs ** order) * fgrid).ravel()) - return deriv + nx, ny = shape + dx = (x.max() - x.min())/(nx - 1) + fx = 2*numpy.pi*numpy.fft.fftfreq(padshape[0], dx) + dy = (y.max() - y.min())/(ny - 1) + fy = 2*numpy.pi*numpy.fft.fftfreq(padshape[1], dy) + return numpy.meshgrid(fy, fx)[::-1] diff --git a/fatiando/gridder.py b/fatiando/gridder.py index 34945802..9d7d9678 100644 --- a/fatiando/gridder.py +++ b/fatiando/gridder.py @@ -29,7 +29,7 @@ ---- """ - +from __future__ import division import numpy import scipy.interpolate @@ -44,14 +44,6 @@ def load_surfer(fname, fmt='ascii'): http://www.goldensoftware.com/products/surfer - According to Surfer structure, x and y are horizontal and vertical - screen-based coordinates respectively. If the grid is in geographic - coordinates, x will be longitude and y latitude. If the coordinates - are cartesian, x will be the easting and y the norting coordinates. - - WARNING: This is opposite to the convention used for Fatiando. - See io_surfer.py in cookbook. - Parameters: * fname : str @@ -62,15 +54,14 @@ def load_surfer(fname, fmt='ascii'): Returns: * x : 1d-array - Value of the horizontal coordinate of each grid point. + Value of the North-South coordinate of each grid point. * y : 1d-array - Value of the vertical coordinate of each grid point. - * grd : 1d-array + Value of the East-West coordinate of each grid point. + * data : 1d-array Values of the field in each grid point. Field can be for example topography, gravity anomaly etc - * shape : tuple = (ny, nx) - The number of points in the vertical and horizontal grid dimensions, - respectively + * shape : tuple = (nx, ny) + The number of points in the x and y grid dimensions, respectively """ assert fmt in ['ascii', 'binary'], "Invalid grid format '%s'. Should be \ @@ -86,71 +77,107 @@ def load_surfer(fname, fmt='ascii'): with open(fname) as ftext: # DSAA is a Surfer ASCII GRD ID id = ftext.readline() - # Read the number of columns (nx) and rows (ny) - nx, ny = [int(s) for s in ftext.readline().split()] - # Read the min/max value of x (columns/longitue) - xmin, xmax = [float(s) for s in ftext.readline().split()] - # Read the min/max value of y(rows/latitude) + # Read the number of columns (ny) and rows (nx) + ny, nx = [int(s) for s in ftext.readline().split()] + shape = (nx, ny) + # Read the min/max value of columns/longitue (y direction) ymin, ymax = [float(s) for s in ftext.readline().split()] - # Read the min/max value of grd - zmin, zmax = [float(s) for s in ftext.readline().split()] + # Read the min/max value of rows/latitude (x direction) + xmin, xmax = [float(s) for s in ftext.readline().split()] + area = (xmin, xmax, ymin, ymax) + # Read the min/max value of grid values + datamin, datamax = [float(s) for s in ftext.readline().split()] data = numpy.fromiter((float(i) for line in ftext for i in line.split()), dtype='f') - grd = numpy.ma.masked_greater_equal(data, 1.70141e+38) - # Create x and y numpy arrays - x = numpy.linspace(xmin, xmax, nx) - y = numpy.linspace(ymin, ymax, ny) - x, y = [tmp.ravel() for tmp in numpy.meshgrid(x, y)] + data = numpy.ma.masked_greater_equal(data, 1.70141e+38) + assert numpy.allclose(datamin, data.min()) \ + and numpy.allclose(datamax, data.max()), \ + "Min and max values of grid don't match ones read from file." \ + + "Read: ({}, {}) Actual: ({}, {})".format( + datamin, datamax, data.min(), data.max()) + # Create x and y coordinate numpy arrays + x, y = regular(area, shape) if fmt == 'binary': raise NotImplementedError( "Binary file support is not implemented yet.") - return x, y, grd, (ny, nx) + return x, y, data, shape def regular(area, shape, z=None): """ - Create a regular grid. Order of the output grid is x varies first, then y. + Create a regular grid. + + The x directions is North-South and y East-West. Imagine the grid as a + matrix with x varying in the lines and y in columns. + + Returned arrays will be flattened to 1D with ``numpy.ravel``. + + .. warning:: + + As of version 0.4, the ``shape`` argument was corrected to be + ``shape = (nx, ny)`` instead of ``shape = (ny, nx)``. + Parameters: * area ``(x1, x2, y1, y2)``: Borders of the grid * shape - Shape of the regular grid, ie ``(ny, nx)``. + Shape of the regular grid, ie ``(nx, ny)``. * z Optional. z coordinate of the grid points. If given, will return an array with the value *z*. Returns: - * ``[xcoords, ycoords]`` + * ``[x, y]`` Numpy arrays with the x and y coordinates of the grid points - * ``[xcoords, ycoords, zcoords]`` + * ``[x, y, z]`` If *z* given. Numpy arrays with the x, y, and z coordinates of the grid points + Examples:: + + >>> x, y = regular((0, 10, 0, 5), (5, 3)) + >>> x + array([ 0. , 0. , 0. , 2.5, 2.5, 2.5, 5. , 5. , 5. , + 7.5, 7.5, 7.5, 10. , 10. , 10. ]) + >>> x.reshape((5, 3)) + array([[ 0. , 0. , 0. ], + [ 2.5, 2.5, 2.5], + [ 5. , 5. , 5. ], + [ 7.5, 7.5, 7.5], + [ 10. , 10. , 10. ]]) + >>> y.reshape((5, 3)) + array([[ 0. , 2.5, 5. ], + [ 0. , 2.5, 5. ], + [ 0. , 2.5, 5. ], + [ 0. , 2.5, 5. ], + [ 0. , 2.5, 5. ]]) + >>> x, y, z = regular((0, 10, 0, 5), (5, 3), z=-10) + >>> z.reshape((5, 3)) + array([[-10., -10., -10.], + [-10., -10., -10.], + [-10., -10., -10.], + [-10., -10., -10.], + [-10., -10., -10.]]) + + """ - ny, nx = shape + nx, ny = shape x1, x2, y1, y2 = area - dy, dx = spacing(area, shape) - x_range = numpy.arange(x1, x2, dx) - y_range = numpy.arange(y1, y2, dy) - # Need to make sure that the number of points in the grid is correct - # because of rounding errors in arange. Sometimes x2 and y2 are included, - # sometimes not - if len(x_range) < nx: - x_range = numpy.append(x_range, x2) - if len(y_range) < ny: - y_range = numpy.append(y_range, y2) - assert len(x_range) == nx, "Failed! x_range doesn't have nx points" - assert len(y_range) == ny, "Failed! y_range doesn't have ny points" - xcoords, ycoords = [mat.ravel() - for mat in numpy.meshgrid(x_range, y_range)] + assert x1 < x2, \ + "Invalid area dimensions {}, {}. x1 must be < x2.".format(x1, x2) + assert y1 < y2, \ + "Invalid area dimensions {}, {}. y1 must be < y2.".format(y1, y2) + xs = numpy.linspace(x1, x2, nx) + ys = numpy.linspace(y1, y2, ny) + # Must pass ys, xs in this order because meshgrid uses the first argument + # for the columns + arrays = numpy.meshgrid(ys, xs)[::-1] if z is not None: - zcoords = z * numpy.ones_like(xcoords) - return [xcoords, ycoords, zcoords] - else: - return [xcoords, ycoords] + arrays.append(z*numpy.ones(nx*ny, dtype=numpy.float)) + return [i.ravel() for i in arrays] def scatter(area, n, z=None, seed=None): @@ -173,22 +200,26 @@ def scatter(area, n, z=None, seed=None): Returns: - * ``[xcoords, ycoords]`` + * ``[x, y]`` Numpy arrays with the x and y coordinates of the points - * ``[xcoords, ycoords, zcoords]`` + * ``[x, y, z]`` If *z* given. Arrays with the x, y, and z coordinates of the points + Examples:: + + >>> x, y = scatter((0, 10, 0, 2), 4, seed=0) + >>> x + array([ 5.48813504, 7.15189366, 6.02763376, 5.44883183]) + >>> y + array([ 0.8473096 , 1.29178823, 0.87517442, 1.783546 ]) + """ x1, x2, y1, y2 = area numpy.random.seed(seed) - xcoords = numpy.random.uniform(x1, x2, n) - ycoords = numpy.random.uniform(y1, y2, n) - numpy.random.seed() + arrays = [numpy.random.uniform(x1, x2, n), numpy.random.uniform(y1, y2, n)] if z is not None: - zcoords = z * numpy.ones(n) - return [xcoords, ycoords, zcoords] - else: - return [xcoords, ycoords] + arrays.append(z*numpy.ones(n)) + return arrays def spacing(area, shape): @@ -200,19 +231,30 @@ def spacing(area, shape): * area ``(x1, x2, y1, y2)``: Borders of the grid * shape - Shape of the regular grid, ie ``(ny, nx)``. + Shape of the regular grid, ie ``(nx, ny)``. Returns: - * ``[dy, dx]`` + * ``[dx, dy]`` Spacing the y and x directions + Examples:: + + >>> spacing((0, 10, 0, 20), (11, 11)) + [1.0, 2.0] + >>> spacing((0, 10, 0, 20), (11, 21)) + [1.0, 1.0] + >>> spacing((0, 10, 0, 20), (5, 21)) + [2.5, 1.0] + >>> spacing((0, 10, 0, 20), (21, 21)) + [0.5, 1.0] + """ x1, x2, y1, y2 = area - ny, nx = shape - dx = float(x2 - x1) / float(nx - 1) - dy = float(y2 - y1) / float(ny - 1) - return [dy, dx] + nx, ny = shape + dx = (x2 - x1)/(nx - 1) + dy = (y2 - y1)/(ny - 1) + return [dx, dy] def interp(x, y, v, shape, area=None, algorithm='cubic', extrapolate=False): @@ -225,8 +267,8 @@ def interp(x, y, v, shape, area=None, algorithm='cubic', extrapolate=False): Arrays with the x and y coordinates of the data points. * v : 1D array Array with the scalar value assigned to the data points. - * shape : tuple = (ny, nx) - Shape of the interpolated regular grid, ie (ny, nx). + * shape : tuple = (nx, ny) + Shape of the interpolated regular grid, ie (nx, ny). * area : tuple = (x1, x2, y1, y2) The are where the data will be interpolated. If None, then will get the area from *x* and *y*. @@ -245,13 +287,11 @@ def interp(x, y, v, shape, area=None, algorithm='cubic', extrapolate=False): """ if algorithm not in ['cubic', 'linear', 'nearest']: raise ValueError("Invalid interpolation algorithm: " + str(algorithm)) - ny, nx = shape + nx, ny = shape if area is None: area = (x.min(), x.max(), y.min(), y.max()) x1, x2, y1, y2 = area - xs = numpy.linspace(x1, x2, nx) - ys = numpy.linspace(y1, y2, ny) - xp, yp = [i.ravel() for i in numpy.meshgrid(xs, ys)] + xp, yp = regular(area, shape) grid = interp_at(x, y, v, xp, yp, algorithm=algorithm, extrapolate=extrapolate) return [xp, yp, grid] diff --git a/test/test_gravmag_euler.py b/test/test_gravmag_euler.py index c6c11ec7..20fc1299 100644 --- a/test/test_gravmag_euler.py +++ b/test/test_gravmag_euler.py @@ -1,7 +1,7 @@ from __future__ import division import numpy as np from fatiando.gravmag.euler import Classic, ExpandingWindow, MovingWindow -from fatiando.gravmag import sphere, transform +from fatiando.gravmag import sphere from fatiando.mesher import Sphere from fatiando import utils, gridder @@ -19,17 +19,22 @@ def setup(): global model, x, y, z, inc, dec, struct_ind, field, xderiv, yderiv, \ zderiv, base, pos inc, dec = -30, 50 - pos = np.array([1000, 1000, 200]) + pos = np.array([1000, 1200, 200]) model = Sphere(pos[0], pos[1], pos[2], 1, {'magnetization': utils.ang2vec(10000, inc, dec)}) struct_ind = 3 - shape = (128, 128) - x, y, z = gridder.regular((0, 3000, 0, 3000), shape, z=-1) + shape = (200, 200) + x, y, z = gridder.regular((0, 3000, 0, 3000), shape, z=-100) base = 10 - field = utils.nt2si(sphere.tf(x, y, z, [model], inc, dec)) + base - xderiv = transform.derivx(x, y, field, shape) - yderiv = transform.derivy(x, y, field, shape) - zderiv = transform.derivz(x, y, field, shape) + field = sphere.tf(x, y, z, [model], inc, dec) + base + # Use finite difference derivatives so that these tests don't depend on the + # performance of the FFT derivatives. + xderiv = (sphere.tf(x + 1, y, z, [model], inc, dec) + - sphere.tf(x - 1, y, z, [model], inc, dec))/2 + yderiv = (sphere.tf(x, y + 1, z, [model], inc, dec) + - sphere.tf(x, y - 1, z, [model], inc, dec))/2 + zderiv = (sphere.tf(x, y, z + 1, [model], inc, dec) + - sphere.tf(x, y, z - 1, [model], inc, dec))/2 def test_euler_classic_sphere_mag(): diff --git a/test/test_gravmag_transform.py b/test/test_gravmag_transform.py new file mode 100644 index 00000000..a73a85fc --- /dev/null +++ b/test/test_gravmag_transform.py @@ -0,0 +1,150 @@ +from __future__ import division +import numpy as np +from numpy.testing import assert_allclose +from fatiando.gravmag import transform, prism +from fatiando import gridder, utils +from fatiando.mesher import Prism + + +def _trim(array, shape, d=20): + "Remove d elements from the edges of an array" + return array.reshape(shape)[d : -d, d : -d].ravel() + + +def test_secont_horizontal_derivatives_fd(): + "gravmag.transform 2nd xy derivatives by finite diff against analytical" + model = [Prism(-1000, 1000, -500, 500, 0, 2000, {'density': 100})] + shape = (300, 300) + x, y, z = gridder.regular([-10000, 10000, -10000, 10000], shape, z=-500) + derivatives = 'xx yy'.split() + grav = prism.potential(x, y, z, model) + for deriv in derivatives: + analytical = getattr(prism, 'g{}'.format(deriv))(x, y, z, model) + func = getattr(transform, 'deriv' + deriv[0]) + calculated = utils.si2eotvos(func(x, y, grav, shape, method='fd', + order=2)) + diff = np.abs(analytical - calculated) + assert np.all(diff/np.abs(analytical).max() <= 0.01), \ + "Failed for g{}. Max: {} Mean: {} STD: {}".format( + deriv, diff.max(), diff.mean(), diff.std()) + + +def test_horizontal_derivatives_fd(): + "gravmag.transform 1st xy derivatives by finite diff against analytical" + model = [Prism(-1000, 1000, -500, 500, 0, 2000, {'density': 100})] + shape = (300, 300) + x, y, z = gridder.regular([-5000, 5000, -5000, 5000], shape, z=-200) + derivatives = 'x y'.split() + grav = utils.mgal2si(prism.gz(x, y, z, model)) + for deriv in derivatives: + analytical = getattr(prism, 'g{}z'.format(deriv))(x, y, z, model) + func = getattr(transform, 'deriv' + deriv) + calculated = utils.si2eotvos(func(x, y, grav, shape, method='fd')) + diff = np.abs(analytical - calculated) + assert np.all(diff <= 0.005*np.abs(analytical).max()), \ + "Failed for g{}. Max: {} Mean: {} STD: {}".format( + deriv, diff.max(), diff.mean(), diff.std()) + + +def test_derivatives_uneven_shape(): + "gravmag.transform FFT derivatives work if grid spacing is uneven" + model = [Prism(-1000, 1000, -500, 500, 0, 2000, {'density': 100})] + shape = (150, 300) + x, y, z = gridder.regular([-10000, 10000, -10000, 10000], shape, z=-100) + grav = utils.mgal2si(prism.gz(x, y, z, model)) + analytical = prism.gzz(x, y, z, model) + calculated = utils.si2eotvos(transform.derivz(x, y, grav, shape, + method='fft')) + diff = _trim(np.abs(analytical - calculated), shape) + assert np.all(diff <= 0.005*np.abs(analytical).max()), \ + "Failed for gzz" + + +def test_gz_derivatives(): + "gravmag.transform FFT 1st derivatives of gz against analytical solutions" + model = [Prism(-1000, 1000, -500, 500, 0, 2000, {'density': 100})] + shape = (300, 300) + x, y, z = gridder.regular([-10000, 10000, -10000, 10000], shape, z=-100) + derivatives = 'x y z'.split() + grav = utils.mgal2si(prism.gz(x, y, z, model)) + for deriv in derivatives: + analytical = getattr(prism, 'g{}z'.format(deriv))(x, y, z, model) + calculated = utils.si2eotvos( + getattr(transform, 'deriv' + deriv)(x, y, grav, shape, + method='fft')) + diff = _trim(np.abs(analytical - calculated), shape) + assert np.all(diff <= 0.005*np.abs(analytical).max()), \ + "Failed for g{}z".format(deriv) + + +def test_gx_derivatives(): + "gravmag.transform FFT 1st derivatives of gx against analytical solutions" + model = [Prism(-1000, 1000, -500, 500, 0, 2000, {'density': 100})] + shape = (300, 300) + x, y, z = gridder.regular([-10000, 10000, -10000, 10000], shape, z=-100) + derivatives = 'x y z'.split() + grav = utils.mgal2si(prism.gx(x, y, z, model)) + for deriv in derivatives: + analytical = getattr(prism, 'gx{}'.format(deriv))(x, y, z, model) + calculated = utils.si2eotvos( + getattr(transform, 'deriv' + deriv)(x, y, grav, shape, + method='fft')) + diff = _trim(np.abs(analytical - calculated), shape) + assert np.all(diff <= 0.005*np.abs(analytical).max()), \ + "Failed for gx{}".format(deriv) + + +def test_gy_derivatives(): + "gravmag.transform FFT 1st derivatives of gy against analytical solutions" + model = [Prism(-1000, 1000, -500, 500, 0, 2000, {'density': 100})] + shape = (300, 300) + x, y, z = gridder.regular([-10000, 10000, -10000, 10000], shape, z=-100) + derivatives = 'x y z'.split() + grav = utils.mgal2si(prism.gy(x, y, z, model)) + for deriv in derivatives: + if deriv == 'x': + func = getattr(prism, 'g{}y'.format(deriv)) + else: + func = getattr(prism, 'gy{}'.format(deriv)) + analytical = func(x, y, z, model) + calculated = utils.si2eotvos( + getattr(transform, 'deriv' + deriv)(x, y, grav, shape, + method='fft')) + diff = _trim(np.abs(analytical - calculated), shape) + assert np.all(diff <= 0.005*np.abs(analytical).max()), \ + "Failed for gy{}".format(deriv) + + +def test_second_derivatives(): + "gravmag.transform FFT second derivatives against analytical solutions" + model = [Prism(-1000, 1000, -500, 500, 0, 2000, {'density': -200})] + shape = (300, 300) + x, y, z = gridder.regular([-10000, 10000, -10000, 10000], shape, z=-100) + derivatives = 'xx yy zz'.split() + pot = prism.potential(x, y, z, model) + for deriv in derivatives: + analytical = getattr(prism, 'g{}'.format(deriv))(x, y, z, model) + calculated = utils.si2eotvos( + getattr(transform, 'deriv' + deriv[0])(x, y, pot, shape, order=2, + method='fft')) + diff = _trim(np.abs(analytical - calculated), shape) + assert np.all(diff <= 0.005*np.abs(analytical).max()), \ + "Failed for g{}. Max: {} Mean: {} STD: {}".format( + deriv, diff.max(), diff.mean(), diff.std()) + + +def test_laplace_from_potential(): + "gravmag.transform 2nd derivatives of potential obey the Laplace equation" + model = [Prism(-1000, 1000, -500, 500, 0, 2000, {'density': 200})] + shape = (300, 300) + x, y, z = gridder.regular([-10000, 10000, -10000, 10000], shape, z=-100) + potential = prism.potential(x, y, z, model) + gxx = utils.si2eotvos(transform.derivx(x, y, potential, shape, order=2, + method='fft')) + gyy = utils.si2eotvos(transform.derivy(x, y, potential, shape, order=2, + method='fft')) + gzz = utils.si2eotvos(transform.derivz(x, y, potential, shape, order=2)) + laplace = _trim(gxx + gyy + gzz, shape) + assert np.all(np.abs(laplace) <= 1e-10), \ + "Max: {} Mean: {} STD: {}".format( + laplace.max(), laplace.mean(), laplace.std())