Euler deconvolution with a moving window

Euler deconvolution attempts to estimate the coordinates of simple (idealized) sources from the input potential field data. There is a strong assumption that the sources have simple geometries, like spheres, vertical pipes, vertical planes, etc. So it wouldn't be much of a surprise if the solutions aren't great when sources are complex.

Let's test the Euler deconvolution using a moving window scheme, a very common approach used in all industry software. This is implemented in geoist.pfm.euler.EulerDeconvMW.

Out:

Centers of the model spheres:
[-1000 -1000  1500]
[1000 1500 1000]
Kept Euler solutions after the moving window scheme:
[[  989.68190633  1558.68892298   985.92770758]
 [ 1028.51412445  1551.02774256  1003.47529929]
 [  990.92900292  1632.75986926  1017.88535173]
 [ -846.38605388  -932.25154989  1546.39530207]
 [-1007.01529935  -958.90553833  1458.84620162]
 [ 1010.61058572  1063.44141734  1161.90617823]
 [-1019.87728776  -999.12690475  1716.28279143]
 [ -991.96984896  -936.38986688  1457.82722913]
 [ 1122.18986393  1140.91674932  1066.8795095 ]
 [ 1013.61672211  1479.48010121  1020.27208383]
 [ 1266.43364353  1268.44969085  1107.51785907]
 [ -887.80662905  -836.62388506  1434.53535244]
 [-1052.03331726  -958.69382461  1488.40740178]
 [  963.03938175  1466.00676155  1014.28009816]
 [ 1041.74392132   813.7900599   1684.26936552]
 [ 1093.48517575  1577.71083791  1080.78498541]
 [-1016.18921626  -907.28379211  1423.30724747]
 [ -482.42934736  -786.30472215  2168.43627994]
 [-1544.30386106  -828.54173105  1459.23359444]
 [ 1038.13879618  1493.20216299   993.48642178]
 [ -142.8704577   -131.93773605  2219.3049965 ]
 [ 2109.53087795  1352.66986754  2203.6361876 ]
 [ 1227.41763046  1215.16966722   722.07230212]
 [  923.48784557  1464.19888115   979.21273841]
 [ 1469.80915257  1504.10875424  1309.18844932]
 [ 1068.36334131  1475.09179417   982.04073591]
 [-1056.96669665  -994.05790313  1376.39399058]
 [-1013.60079694 -1060.82280253  1545.92013331]
 [ 1390.48535066  -111.53610308  2492.03192256]
 [-1013.40254268 -1056.49300094  1495.7105859 ]
 [-1019.76315975 -1055.24604826  1517.35776498]
 [  -68.08398933   379.48764513  2073.84015287]
 [-1074.22116872 -1097.53380481  1534.43104553]
 [-1012.3401741  -1036.37141188  1547.53572386]
 [-1165.61264459  -942.29808427  1224.31423615]
 [ -149.51962929   827.97251027  1969.8913013 ]
 [-1007.93868833 -1058.01948708  1470.22068576]
 [ 1306.04396369  1457.36741652   819.17052122]
 [ 1130.52465623  1634.56447324  1252.39418085]
 [  446.97334971   771.84872308  2968.32225288]]
D:\MyWeb\geoistdoc\examples\tutorials\euler_moving_window.py:82: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
  plt.show()

from geoist.pfm import sphere, pftrans, euler, giutils
from geoist import gridder
from geoist.inversion import geometry
from geoist.vis import giplt
import matplotlib.pyplot as plt

# Make some synthetic magnetic data to test our Euler deconvolution.
# The regional field
inc, dec = -45, 0
# Make a model of two spheres magnetized by induction only
model = [
    geometry.Sphere(x=-1000, y=-1000, z=1500, radius=1000,
                  props={'magnetization': giutils.ang2vec(2, inc, dec)}),
    geometry.Sphere(x=1000, y=1500, z=1000, radius=1000,
                  props={'magnetization': giutils.ang2vec(1, inc, dec)})]

print("Centers of the model spheres:")
print(model[0].center)
print(model[1].center)

# Generate some magnetic data from the model
shape = (100, 100)
area = [-5000, 5000, -5000, 5000]
x, y, z = gridder.regular(area, shape, z=-150)
data = sphere.tf(x, y, z, model, inc, dec)

# We also need the derivatives of our data
xderiv = pftrans.derivx(x, y, data, shape)
yderiv = pftrans.derivy(x, y, data, shape)
zderiv = pftrans.derivz(x, y, data, shape)

# Now we can run our Euler deconv solver on a moving window over the data.
# Each window will produce an estimated point for the source.
# We use a structural index of 3 to indicate that we think the sources are
# spheres.

# Run the Euler deconvolution on moving windows to produce a set of solutions
# by running the solver on 10 x 10 windows of size 1000 x 1000 m
solver = euler.EulerDeconvMW(x, y, z, data, xderiv, yderiv, zderiv,
                             structural_index=3, windows=(10, 20),
                             size=(1000, 500))
# Use the fit() method to obtain the estimates
solver.fit()

# The estimated positions are stored as a list of [x, y, z] coordinates
# (actually a 2D numpy array)
print('Kept Euler solutions after the moving window scheme:')
print(solver.estimate_)

# Plot the solutions on top of the magnetic data. Remember that the true depths
# of the center of these sources is 1500 m and 1000 m.

plt.figure(figsize=(6, 5))
plt.title('Euler deconvolution with a moving window')
plt.contourf(y.reshape(shape), x.reshape(shape), data.reshape(shape), 30,
             cmap="RdBu_r")
plt.scatter(solver.estimate_[:, 1], solver.estimate_[:, 0],
            s=50, c=solver.estimate_[:, 2], cmap='cubehelix')
plt.colorbar(pad=0).set_label('Depth (m)')
plt.xlim(area[2:])
plt.ylim(area[:2])
plt.tight_layout()
plt.show()