gempy-project
export gempy results to pygimli/unstructured mesh
Dear gempy team,
First of all thanks a lot for the open source code you develop: it is really useful. My name is Giulia de Pasquale and I’m a postdoc researcher at CEAZA (Centro de Estudio Avanzado en Zonas Aridas) in Chile. I’m currently using your code for modelling some aquifer layer info gained from resistivity soundings (TEM) and I’m interested in interpolate the results obtained with gempy to an unstructured mesh in pygimli or gempy if it allows it (I was not able to find examples that where not cubic/parallelepiped regular mesh on the website...).
Not sure if you already faced this kind of request, but in case not what would be really helpful for me are some instruction on how to extrapolate the 3D regular grid cell-centres and interpolated values associated to those into list or numpy arrays so that I could use them for an interpolation in a grid adapted to my field geometry.
Thank you in advance Giulia de Pasquale
Hi Giulia,
if I understand you correctly all what you need is to use custom_grids. gempy can interpolate in any xyz coordinate you feed it. There is this method geo_model.set_custom_grid (https://docs.gempy.org/Model/gempy.core.model.ImplicitCoKriging.html?highlight=set_custom_grid#gempy.core.model.ImplicitCoKriging.set_custom_grid) that let you simply pass a 2D numpy array. Then when you do compute model, you will populate geo_model.solutions.custom_grid (not sure if that is the exact name) with the values corresponding to the defined custom grid. Alternatively (it does the same), you can pass the numpy array with the coordinates to the argument at in gempy.compute_model (https://docs.gempy.org/GemPy%20API/gempy.compute_model.html?highlight=compute_model#gempy.compute_model).
Also we are currently working in this library subsurface to make easier passing the data from library to library (https://softwareunderground.github.io/subsurface/external/external_examples.html). Some of the functionality are not in the pip releases (you would have to use the main branches from github) yet but it is worth it to keep an eye on this development.
Let me know if this helps
Cheers
Hi Miguel, thank you for your answer: I think costume grid is what I'm looking for: is it able to return an unstructured grid?
I'm currently trying to set this costume grid and compute my model on it but the resulting model only has the points I set to define my surfaces. Also it seems is not possible to set the topography in case of custom grid(?-I was able with no problems with regular grid through gdal)
My code (without topography) look like this:
(I first create an unstructured 3D mesh with pygimli: mesh 3d, and I have data points with values of X,Y,Z which I use to create and orienting the surfaces)
xm3d=np.array(pg.x(mesh3d.cellCenters())).tolist() ym3d=np.array(pg.y(mesh3d.cellCenters())).tolist() zm3d=np.array(pg.z(mesh3d.cellCenters())).tolist() C3D=np.array([xm3d,ym3d,zm3d]).T #this is the 2D XYZ vector
geo_model = gp.create_model('Rio_Limari_II')
geo_model.set_custom_grid(C3D)
geo_model.add_surfaces(surface_list=['surface1','surface2','surface3','surface4','surface5','basement'])
for i in range(1,Nl+1):
name='surface'+str(i)
for j in range(Nd):
geo_model.add_surface_points(X=X[j,i], Y=Y[j,i], Z=z[j,i], surface=name)
geo_model.add_orientations(X=X[j,i], Y=Y[j,i], Z=z[j,i], surface=name, pole_vector=(0, 0, 1))
gp.set_interpolator(geo_model, theano_optimizer='fast_compile', verbose=[])
gp.get_data(geo_model, 'kriging')
sol = gp.compute_model(geo_model)