Skip to content

Commit

Permalink
mc
Browse files Browse the repository at this point in the history
  • Loading branch information
leonfoks committed Jun 13, 2024
2 parents c128682 + 7bc9cfb commit 820e8a4
Show file tree
Hide file tree
Showing 571 changed files with 13,907 additions and 8,920 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -136,7 +136,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.12"
"version": "3.12.3"
}
},
"nbformat": 4,
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -277,7 +277,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.12"
"version": "3.12.3"
}
},
"nbformat": 4,
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -39,9 +39,9 @@
# For more information about the time domain data set, see :ref:`Time domain dataset`

# The data file name
dataFile=dataFolder + 'skytem_512_saline_clay.csv'
dataFile=dataFolder + 'skytem_saline_clay.csv'
# The EM system file name
systemFile=[dataFolder + 'SkytemHM_512.stm', dataFolder + 'SkytemLM_512.stm']
systemFile=[dataFolder + 'SkytemHM.stm', dataFolder + 'SkytemLM.stm']

#%%
# Initialize and read an EM data set
Expand Down Expand Up @@ -108,7 +108,7 @@
#%%
# Plot the misfits for a range of half space conductivities
plt.figure()
_ = tdp.plotHalfSpaceResponses(-6.0, 4.0, 200)
_ = tdp.plot_halfspace_responses(-6.0, 4.0, 200)
plt.title("Halfspace responses")

#%%
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.12"
"version": "3.12.3"
}
},
"nbformat": 4,
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -450,7 +450,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.12"
"version": "3.12.3"
}
},
"nbformat": 4,
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -132,7 +132,7 @@
# Plot the misfits for a range of half space conductivities
plt.figure()
plt.subplot(1, 2, 1)
_ = tdp.plotHalfSpaceResponses(-6.0, 4.0, 200)
_ = tdp.plot_halfspace_responses(-6.0, 4.0, 200)
plt.title("Halfspace responses")

#%%
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -38,12 +38,15 @@ def plot_2d_summary(folder, data_type, model_type):
"cmap" : 'jet'
}

import warnings
warnings.filterwarnings('error')

fig = plt.figure(figsize=(16, 8))
plt.suptitle("{} {}".format(data_type, model_type))
gs0 = fig.add_gridspec(6, 2, hspace=1.0)

true_model = Model.create_synthetic_model(model_type)
true_model.mesh.y_edges = true_model.mesh.y_edges / 4.1
true_model.mesh.y_edges = true_model.mesh.y_edges / 10.0

kwargs['vmin'] = np.log10(np.min(true_model.values))
kwargs['vmax'] = np.log10(np.max(true_model.values))
Expand Down Expand Up @@ -134,17 +137,16 @@ def plot_2d_summary(folder, data_type, model_type):
results_2d.plot_data_elevation(linewidth=0.3, ax=ax1);
results_2d.plot_elevation(linewidth=0.3, ax=ax1);


plt.show()
# plt.savefig('{}_{}.png'.format(data_type, model_type), dpi=300)
# plt.show()
plt.savefig('{}_{}.png'.format(data_type, model_type), dpi=300)

if __name__ == '__main__':
models = ['glacial', 'saline_clay', 'resistive_dolomites', 'resistive_basement', 'coastal_salt_water', 'ice_over_salt_water']

for model in models:
try:
# try:
plot_2d_summary("../../../Parallel_Inference/", "resolve", model)
except Exception as e:
print(model)
print(e)
pass
# except Exception as e:
# print(model)
# print(e)
# pass
Original file line number Diff line number Diff line change
Expand Up @@ -205,7 +205,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.12"
"version": "3.12.3"
}
},
"nbformat": 4,
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,6 @@
# The StatArray class has been built so that we may easily
# attach not only names and units, but statistical distributions too.
# We won't go into too much detail about the different distribution
# classes here so check out the :ref:`Distribution Class` for a better description.
#
# Two types of distributions can be attached to the StatArray.
#
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -134,8 +134,8 @@ def plot_2d_summary(folder, data_type, model_type):
results_2d.plot_data_elevation(linewidth=0.3, ax=ax1);
results_2d.plot_elevation(linewidth=0.3, ax=ax1);

plt.show()
# plt.savefig('{}_{}.png'.format(data_type, model_type), dpi=300)
# plt.show()
plt.savefig('{}_{}.png'.format(data_type, model_type), dpi=300)


if __name__ == '__main__':
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -320,7 +320,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.12"
"version": "3.12.3"
}
},
"nbformat": 4,
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -5,13 +5,13 @@

#%%

from geobipy import PointCloud3D
from geobipy import Point
from os.path import join
import numpy as np
import matplotlib.pyplot as plt
import h5py

nPoints = 10000
nPoints = 200

#%%
# Create a quick test example using random points
Expand All @@ -20,15 +20,15 @@
y = -np.abs((2.0 * np.random.rand(nPoints)) - 1.0)
z = x * (1.0 - x) * np.cos(np.pi * x) * np.sin(np.pi * y)

PC3D = PointCloud3D(x=x, y=y, z=z)
PC3D = Point(x=x, y=y, z=z)

#%%
# Append pointclouds together
x = np.abs((2.0 * np.random.rand(nPoints)) - 1.0)
y = np.abs((2.0 * np.random.rand(nPoints)) - 1.0)
z = x * (1.0 - x) * np.cos(np.pi * x) * np.sin(np.pi * y)

Other_PC = PointCloud3D(x=x, y=y, z=z)
Other_PC = Point(x=x, y=y, z=z)
PC3D.append(Other_PC)

#%%
Expand All @@ -41,7 +41,6 @@

Point = PC3D[50]
# Print the point to the screen
print(Point)

#%%
# Plot the locations with Height as colour
Expand All @@ -58,7 +57,7 @@

#%%
# Interpolate the points to a 2D rectilinear mesh
mesh, dum = PC3D.interpolate(0.01, 0.01, values=PC3D.z, method='mc', mask=0.03)
mesh, dum = PC3D.interpolate(0.01, 0.01, values=PC3D.z, method='sibson', mask=0.03)

# We can save that mesh to VTK
PC3D.to_vtk('pc3d.vtk')
Expand All @@ -68,44 +67,50 @@
# Grid the points using a triangulated CloughTocher, or minimum curvature interpolation

plt.figure()
plt.subplot(321)
plt.subplot(331)
PC3D.map(dx=0.01, dy=0.01, method='ct')
plt.subplot(322)
plt.subplot(332)
PC3D.map(dx=0.01, dy=0.01, method='mc')
plt.subplot(333)
PC3D.map(dx=0.01, dy=0.01, method='sibson')

plt.subplot(323)
plt.subplot(334)
PC3D.map(dx=0.01, dy=0.01, method='ct', mask=0.03)
plt.subplot(324)
plt.subplot(335)
PC3D.map(dx=0.01, dy=0.01, method='mc', mask=0.03)
plt.subplot(336)
PC3D.map(dx=0.01, dy=0.01, method='sibson', mask=0.03)
#%%
# For lots of points, these surfaces can look noisy. Using a block filter will help
PCsub = PC3D.block_median(0.005, 0.005)
plt.subplot(325)
plt.subplot(337)
PCsub.map(dx=0.01, dy=0.01, method='ct', mask=0.03)
plt.subplot(326)
plt.subplot(338)
PCsub.map(dx=0.01, dy=0.01, method='mc', mask=0.03)
plt.subplot(339)
PCsub.map(dx=0.01, dy=0.01, method='sibson', mask=0.03)


#%%
# We can perform spatial searches on the 3D point cloud

PC3D.setKdTree(nDims=2)
PC3D.set_kdtree(ndim=2)
p = PC3D.nearest((0.0,0.0), k=200, p=2, radius=0.3)

#%%
# .nearest returns the distances and indices into the point cloud of the nearest points.
# We can then obtain those points as another point cloud

pNear = PC3D[p[1]]
plt.figure()
ax1 = plt.subplot(1,2,1)
pNear.scatter2D()
plt.plot(0.0, 0.0, 'x')
plt.subplot(1,2,2, sharex=ax1, sharey=ax1)
ax, sc, cb = PC3D.scatter2D(edgecolor='k')
searchRadius = plt.Circle((0.0, 0.0), 0.3, color='b', fill=False)
ax.add_artist(searchRadius)
plt.plot(0.0, 0.0, 'x')
# pNear = PC3D[p[1]]
# plt.figure()
# ax1 = plt.subplot(1,2,1)
# pNear.scatter2D()
# plt.plot(0.0, 0.0, 'x')
# plt.subplot(1,2,2, sharex=ax1, sharey=ax1)
# ax, sc, cb = PC3D.scatter2D(edgecolor='k')
# searchRadius = plt.Circle((0.0, 0.0), 0.3, color='b', fill=False)
# ax.add_artist(searchRadius)
# plt.plot(0.0, 0.0, 'x')

#%%
# Read in the xyz co-ordinates in columns 2,3,4 from a file. Skip 1 header line.
Expand All @@ -126,9 +131,9 @@
PC3D.writeHdf(f, 'test')

with h5py.File('test.h5', 'r') as f:
PC3D1 = PointCloud3D.fromHdf(f['test'])
PC3D1 = Point.fromHdf(f['test'])

with h5py.File('test.h5', 'r') as f:
point = PointCloud3D.fromHdf(f['test'], index=0)
point = Point.fromHdf(f['test'], index=0)

plt.show()
Original file line number Diff line number Diff line change
Expand Up @@ -320,7 +320,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.12"
"version": "3.12.3"
}
},
"nbformat": 4,
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -18,10 +18,10 @@
#%%
dataFolder = "..//..//supplementary//data//"
# The data file name
dataFiles=dataFolder + 'skytem_304_saline_clay.csv'
dataFiles=dataFolder + 'skytem_saline_clay.csv'
# dataFiles = dataFolder + 'Skytem.csv'
# The EM system file name
systemFiles=[dataFolder + 'SkytemHM_304.stm', dataFolder + 'SkytemLM_304.stm']
systemFiles=[dataFolder + 'SkytemHM.stm', dataFolder + 'SkytemLM.stm']

from pathlib import Path
for f in systemFiles[:1]:
Expand Down Expand Up @@ -56,7 +56,7 @@

#%%
plt.figure(5)
ax = TD.scatter2D(s=1.0, c=TD.secondary_field[:, TD.channel_index(system=0, channel=6)], equalize=True, log=10)
ax = TD.scatter2D(c=TD.secondary_field[:, TD.channel_index(system=0, channel=6)], log=10)
plt.axis('equal')


Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@
},
"outputs": [],
"source": [
"# Make a test model with 10 layers, and increasing parameter values\nnLayers = 2\npar = StatArray(np.linspace(0.001, 0.02, nLayers), \"Conductivity\", \"$\\\\frac{S}{m}$\")\nthk = StatArray(np.full(nLayers, fill_value=10.0))\nthk[-1] = np.inf\nmesh = RectilinearMesh1D(widths = thk)\n\nmod = Model(mesh = mesh, values=par)\n# mod = Model1D(parameters=par, widths=thk)\n\nplt.figure()\nmod.plotGrid(transpose=True, flip=True)"
"# Make a test model with 10 layers, and increasing parameter values\nnLayers = 2\npar = StatArray(np.linspace(0.001, 0.02, nLayers), \"Conductivity\", \"$\\\\frac{S}{m}$\")\nthk = StatArray(np.full(nLayers, fill_value=10.0))\nthk[-1] = np.inf\nmesh = RectilinearMesh1D(widths = thk)\n\nmod = Model(mesh = mesh, values=par)\n\nplt.figure()\nmod.plotGrid(transpose=True, flip=True)"
]
},
{
Expand Down Expand Up @@ -219,7 +219,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.12"
"version": "3.12.3"
}
},
"nbformat": 4,
Expand Down
Loading

0 comments on commit 820e8a4

Please sign in to comment.