2d gridding

[emurphy222000-jpg]

Do gempy produce standard 2-d grids or only 3-d grids. So I simply want to grid a sparse set of depth points(Z values) vs x,y coordinates and display that. Is it possible. If so, hints?
Thanks

[emurphy222000-jpg]

Additionally, the grids displayed seem to indicate the x,y,z nodes where values will be interpolated. Where are those interpolated values for each node stored.

[javoha]

Hello @emurphy222000-jpg,

do I understand correctly that you don't want to compute a geological model? Basically you just want to interpolate your Z values on a 2D grid? Maybe you can post a short sketch of what you are exactly looking for. Sounds like classical interpolation (e.g Inverse distance, Nearest neighbor, kriging) might just do the trick, but maybe I am missing an important point here.

[emurphy222000-jpg]

Thank you so much for your reply. "Basically you just want to interpolate your Z values on a 2D grid? "yes, exactly so. It may be what I want can't be done with gempy. I would just like to know one way or another since I have spent considerable time with it. So we have a product that takes random X,Y,Z values, with Z as depth to the top of a subsurface point in a horizon. (The points are taken from tops or markers for a number of wells in a given region). What we do now is grid these points with standard gridding algorithms. A short coming of these techniques is that they do not honor faults (or fault polygons). With gempy it is possible to model faults embedded in horizons, accurately. In order to display these in our program I need to convert the output of the gempy routine to an X,Y,Z 2d grid. I had thought I might be able to use vertices and snap those nodes to an orthogonal grid. Alternatively, Is there a way to grab the X,Y,Z for the lowest Z value directly from gempy. I have not been able to find where that possibility exists as of yet.Thank you very much for your reply, it would be of great help if we could clear up some of these questions.Ed Murphy On Thursday, May 27, 2021, 03:03:26 AM CDT, Jan von Harten ***@***.***> wrote: Hello @emurphy222000-jpg, do I understand correctly that you don't want to compute a geological model? Basically you just want to interpolate your Z values on a 2D grid? Maybe you can post a short sketch of what you are exactly looking for. Sounds like classical interpolation (e.g Inverse distance, Nearest neighbor, kriging) might just do the trick, but maybe I am missing an important point here. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[javoha]

Hi again,
so I think I got your question - and I think there is no generic gempy solution. I hacked a small workaround, not nice code or anything but maybe it does what you are looking for. I used this (https://docs.gempy.org/examples/geometries/7_combination.html#sphx-glr-examples-geometries-7-combination-py) example model for the test.

def extract_domain(sol, unit):
    '''
    Extract domain coordinates from gempy model by unit name
    arguments:
        sol: Gempy solution object.
        unit: id of gempy unit to be extracted
    returns:
        dom_x, dom_y, dom_z: coordinates of domain
    '''

    # new version with rounding, definitely necessary
    rounded_lithblock = sol.lith_block.round(0)
    rounded_lithblock = rounded_lithblock.astype(int)

    # mask by array of input surfaces (by id, can be from different series)
    mask = np.isin(rounded_lithblock, unit)

    # get coordinates by mask
    #krig_lith = sol.lith_block[mask]
    dom_grid = sol.grid.values[mask]

    return dom_grid

grid = extract_domain(sol, 3) # extract single unit by id

# Following part is to find lowest z value for each unique xy coordinate

# very much hacked and not very nice, but hopefully does the job
# make zeros array as start
res = np.zeros((1,3))

# loop over unique x and y coordiantes and save minimum
for i in np.unique(grid[:,0]):
    for j in np.unique(grid[:,1]):

        mask = (grid[:,0]==i)&(grid[:,1]==j)
        
        temp1 = grid[mask]
        temp2 = np.where(temp1[:,2] == np.amin(temp1[:,2]))
        temp3 = temp1[temp2]
        # save result
        res = np.vstack((res, temp3))
# remove first zero entry
res = res[1:,:]

# Plot result as scatter to show grid
fig = plt.figure(figsize=(15,10))
ax = fig.add_subplot(111)
c = ax.scatter(res[:,0], res[:,1], c=res[:,2], s=3, cmap="viridis")
ax.set_aspect('equal')
plt.colorbar(c, shrink=0.5)
plt.show()        

This results in the following plot:

Is this what you wanted?

[emurphy222000-jpg]

Oh goodness, you are fantastic. Let me give it a run through and I will report back.Thanks so much for the great effort.Ed On Thursday, May 27, 2021, 01:47:05 PM CDT, Jan von Harten ***@***.***> wrote: Hi again, so I think I got your question - and I think there is no generig gempy solution. I hacked a small workaround, not nice code or anything but maybe it does what you are looking for. I used this (https://docs.gempy.org/examples/geometries/7_combination.html#sphx-glr-examples-geometries-7-combination-py) example model for the test. def extract_domain(sol, unit): ''' Extract domain coordinates from gempy model by unit name arguments: sol: Gempy solution object. unit: id of gempy unit to be extracted returns: dom_x, dom_y, dom_z: coordinates of domain ''' # new version with rounding, definitely necessary rounded_lithblock = sol.lith_block.round(0) rounded_lithblock = rounded_lithblock.astype(int) # mask by array of input surfaces (by id, can be from different series) mask = np.isin(rounded_lithblock, unit) # get coordinates by mask #krig_lith = sol.lith_block[mask] dom_grid = sol.grid.values[mask] return dom_grid grid = extract_domain(sol, 3) # extract single unit by id # Following part is to find lowest z value for each unique xy coordinate # very much hacked and not very nice, but hopefully does the job # make zeros array as start res = np.zeros((1,3)) # loop over unique x and y coordiantes and save minimum for i in np.unique(grid[:,0]): for j in np.unique(grid[:,1]): mask = (grid[:,0]==i)&(grid[:,1]==j) temp1 = grid[mask] temp2 = np.where(temp1[:,2] == np.amin(temp1[:,2])) temp3 = temp1[temp2] # save result res = np.vstack((res, temp3)) # remove first zero entry res = res[1:,:] # Plot result as scatter to show grid fig = plt.figure(figsize=(15,10)) ax = fig.add_subplot(111) c = ax.scatter(res[:,0], res[:,1], c=res[:,2], s=3, cmap="viridis") ax.set_aspect('equal') plt.colorbar(c, shrink=0.5) plt.show() This results in the following plot: Is this what you wanted? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[emurphy222000-jpg]

Not finished, but so far this is working out great!!!!THANK YOU. Ed On Thursday, May 27, 2021, 01:47:05 PM CDT, Jan von Harten ***@***.***> wrote: Hi again, so I think I got your question - and I think there is no generig gempy solution. I hacked a small workaround, not nice code or anything but maybe it does what you are looking for. I used this (https://docs.gempy.org/examples/geometries/7_combination.html#sphx-glr-examples-geometries-7-combination-py) example model for the test. def extract_domain(sol, unit): ''' Extract domain coordinates from gempy model by unit name arguments: sol: Gempy solution object. unit: id of gempy unit to be extracted returns: dom_x, dom_y, dom_z: coordinates of domain ''' # new version with rounding, definitely necessary rounded_lithblock = sol.lith_block.round(0) rounded_lithblock = rounded_lithblock.astype(int) # mask by array of input surfaces (by id, can be from different series) mask = np.isin(rounded_lithblock, unit) # get coordinates by mask #krig_lith = sol.lith_block[mask] dom_grid = sol.grid.values[mask] return dom_grid grid = extract_domain(sol, 3) # extract single unit by id # Following part is to find lowest z value for each unique xy coordinate # very much hacked and not very nice, but hopefully does the job # make zeros array as start res = np.zeros((1,3)) # loop over unique x and y coordiantes and save minimum for i in np.unique(grid[:,0]): for j in np.unique(grid[:,1]): mask = (grid[:,0]==i)&(grid[:,1]==j) temp1 = grid[mask] temp2 = np.where(temp1[:,2] == np.amin(temp1[:,2])) temp3 = temp1[temp2] # save result res = np.vstack((res, temp3)) # remove first zero entry res = res[1:,:] # Plot result as scatter to show grid fig = plt.figure(figsize=(15,10)) ax = fig.add_subplot(111) c = ax.scatter(res[:,0], res[:,1], c=res[:,2], s=3, cmap="viridis") ax.set_aspect('equal') plt.colorbar(c, shrink=0.5) plt.show() This results in the following plot: Is this what you wanted? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[javoha]

Great to hear and you are very welcome - if this leads to any results that you can share, feel free to post it here - Would be nice to see and maybe also helps other people. Good luck with your work!

[emurphy222000-jpg]

Will definitely do that. I am still working on moving results from your algorithm into a format that our application reads. I think I am nearly there.Ed On Wednesday, June 2, 2021, 09:19:19 AM CDT, Jan von Harten ***@***.***> wrote: Great to hear and you are very welcome - if this leads to any results that you can share, feel free to post it here - Would be nice to see and maybe also helps other people. Good luck with your work! — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[emurphy222000-jpg]

Just a follow up question regarding a piece of your code. # loop over unique x and y coordiantes and save minimum for i in np.unique(grid[:,0]): for j in np.unique(grid[:,1]): mask = (grid[:,0]==i)&(grid[:,1]==j) temp1 = grid[mask] temp2 = np.where(temp1[:,2] == np.amin(temp1[:,2])) temp3 = temp1[temp2] # save result res = np.vstack((res, temp3)) When I run my data through this the calculation for mask it sometimes return mask as all false and as a result temp1 is defined as an ndarray with either zero rows or zero columns. Is this expected behavior with real data or does it indicate a problem with my code?Again, thanks so much. This has been a great help.Ed On Wednesday, June 2, 2021, 09:19:19 AM CDT, Jan von Harten ***@***.***> wrote: Great to hear and you are very welcome - if this leads to any results that you can share, feel free to post it here - Would be nice to see and maybe also helps other people. Good luck with your work! — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[javoha]

Mmh this should not happen if the unit you extracted contains any grid points. What does the extracted grid look like in these cases?

[emurphy222000-jpg]

Hi Jan:Thanks again for your reply. I am sending my code and the data file. If you have time to run it, I would be grateful. You just need to change the path to the data file and it should run in jupyter. You will note that I have a trap for when the  no data found condition occurs. In that case I just put in a fake value that represents an average value. So obviously that affects the look of the grid.The resultant grid looks as follows: Thank you so much, Jan.Ed On Friday, June 4, 2021, 01:40:18 AM CDT, Jan von Harten ***@***.***> wrote: Mmh this should not happen if the unit you exrtacted contains any grid points. What does the extracted grid look like in these cases? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe. import gempy as gp import numpy as np import pandas as pd import matplotlib.pyplot as plt def extract_domain(sol, unit): ''' Extract domain coordinates from gempy model by unit name arguments: sol: Gempy solution object. unit: id of gempy unit to be extracted returns: dom_x, dom_y, dom_z: coordinates of domain ''' # new version with rounding, definitely necessary rounded_lithblock = sol.lith_block.round(0) rounded_lithblock = rounded_lithblock.astype(int) # mask by array of input surfaces (by id, can be from different series) mask = np.isin(rounded_lithblock, unit) # get coordinates by mask # krig_lith = sol.lith_block[mask] dom_grid = sol.grid.values[mask] # dom_grid = sol.grid.values return dom_grid dfs = [] dfs.append(pd.read_csv("C:\\Users\\17138\\geoLOGIC\\gridding\\gempy\\input_gempy1.csv", sep=',', names=['X', 'Y', 'Z'],header=1)) minx = min(dfs[0]['X']) maxx = max(dfs[0]['X']) miny = min(dfs[0]['Y']) maxy = max(dfs[0]['Y']) minz = min(dfs[0]['Z']) maxz = max(dfs[0]['Z']) dfs_ele = dfs[0] geo_model = gp.create_model('Tiberius') # Importing the data from csv files and settign extent and resolution geo_model = gp.init_data(geo_model, extent=[minx, maxx, miny, maxy, minz, maxz], resolution=[95, 76, 50], add_basement=True) gp.set_interpolator(geo_model, theano_optimizer='fast_compile', verbose=[]) geo_model.set_default_surfaces() # p2d = gp.plot_2d(geo_model) for index, row in dfs_ele.iterrows(): geo_model.add_surface_points(X=row['X'], Y=row['Y'], Z=row['Z'], surface='surface1') #geo_model.add_orientations(X=-11535952.56, Y=6501871.555, Z=-67.2, surface='surface1', pole_vector=(.22, .45, .3)) geo_model.add_orientations(X=759376.548, Y=7037711.6, Z=-1410.54, surface='surface1', pole_vector=(.22, .45, .3)) sol = gp.compute_model(geo_model) occurrences = np.count_nonzero(geo_model.solutions.lith_block == 1) uu = np.count_nonzero(geo_model.solutions.lith_block == 2) grid = extract_domain(sol, 1) # extract single unit by id res = np.zeros((1, 3)) # loop over unique x and y coordiantes and save minimum for i in np.unique(grid[:, 0]): jcntr = 0 lastx = 0 lasty = 0 lastx = 0. lasty = 0. for j in np.unique(grid[:, 1]): mask = (grid[:, 0] == i) & (grid[:, 1] == j) temp1 = grid[mask] trc, trc1 = temp1.shape if trc != 0 and trc1 != 0: temp2 = np.where(temp1[:, 2] == np.amin(temp1[:, 2])) temp3 = temp1[temp2] # print("jcntr: ", jcntr, " temp3: ", temp3) else: # temp3=[0.,0.,-9999999.0 temp3=[759376.548, 7037711.6, -1410.54] # temp3 = [0., 0., .0] # print("null, jcntr: ", jcntr, " temp3: ", temp3) # save result res = np.vstack((res, temp3)) jcntr = jcntr + 1 temp1 = grid[mask] # remove first zero entry res = res[1:, :] # Plot result as scatter to show grid fig = plt.figure(figsize=(15, 10)) ax = fig.add_subplot(111) c = ax.scatter(res[:, 0], res[:, 1], c=res[:, 2], s=3, cmap="viridis") ax.set_aspect('equal') plt.colorbar(c, shrink=0.5) plt.show()

[emurphy222000-jpg]

Hi JanI have updated my function to more closely follow the model7 gempy example code. I think this brings me more in line with your example. It still has about 1300 location where it does not find data. Thank you once again.Ed On Friday, June 4, 2021, 01:40:18 AM CDT, Jan von Harten ***@***.***> wrote: Mmh this should not happen if the unit you exrtacted contains any grid points. What does the extracted grid look like in these cases? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe. import gempy as gp import numpy as np import pandas as pd import matplotlib.pyplot as plt def extract_domain(sol, unit): ''' Extract domain coordinates from gempy model by unit name arguments: sol: Gempy solution object. unit: id of gempy unit to be extracted returns: dom_x, dom_y, dom_z: coordinates of domain ''' # new version with rounding, definitely necessary rounded_lithblock = sol.lith_block.round(0) rounded_lithblock = rounded_lithblock.astype(int) # mask by array of input surfaces (by id, can be from different series) mask = np.isin(rounded_lithblock, unit) # get coordinates by mask # krig_lith = sol.lith_block[mask] dom_grid = sol.grid.values[mask] # dom_grid = sol.grid.values return dom_grid pd.set_option('precision', 2) data_path="C:\\Users\\17138\\geoLOGIC\\gridding\\test\\" path_to_data = data_path geo_data = gp.create_data('combination', extent=[755589.109, 780838.701, 7034865.1, 7053841.76, -1735.9, -1271.1], resolution=[125, 50, 50], path_o=path_to_data + "input_gempy1_orientations.csv", path_i=path_to_data + "input_gempy1.csv") geo_data.get_data() gp.map_stack_to_surfaces(geo_data, {"Strat_Series1": ('surface1'), "Basement_Series":('basement')}) gp.plot_2d(geo_data, direction='y') interp_data = gp.set_interpolator(geo_data, theano_optimizer='fast_compile') # %% sol = gp.compute_model(geo_data) # %% # Displaying the result in x and y direction: # # %% gp.plot_2d(geo_data, cell_number=20, direction='y', show_data=False, show_boundaries=True) grid = extract_domain(sol, 1) # extract single unit by id # Following part is to find lowest z value for each unique xy coordinate # very much hacked and not very nice, but hopefully does the job # make zeros array as start res = np.zeros((1,3)) # loop over unique x and y coordiantes and save minimum jcntr=0 for i in np.unique(grid[:,0]): for j in np.unique(grid[:,1]): mask = (grid[:,0]==i)&(grid[:,1]==j) temp1 = grid[mask] trc, trc1=temp1.shape if trc!= 0 and trc1 != 0 : temp2 = np.where(temp1[:,2] == np.amin(temp1[:,2])) temp3 = temp1[temp2] else: print("no value found for number: ",jcntr) jcntr=jcntr+ 1 # save result res = np.vstack((res, temp3)) # remove first zero entry res = res[1:,:] fig = plt.figure(figsize=(15, 10)) ax = fig.add_subplot(111) c = ax.scatter(res[:, 0], res[:, 1], c=res[:, 2], s=3, cmap="viridis") ax.set_aspect('equal') plt.colorbar(c, shrink=0.5) plt.show()

[javoha]

Hi,
so what should not happen is that it does not find any minumum z locations for a layer (surface in Gempy), there might be gaps though, like in my example, where a layer does not exist for certain xy coordinates (in the example this is due to a unconformity, but might also be because of faulting, a pinch out, etc).
Where did you upload your data, then I could maybe testrun it - didnt fiind it in your post.

[emurphy222000-jpg]

Hi Jan. I have updated my function to more closely follow the model7 gempy example code. I think this brings me more in line with your example. It still has about 1300 location where it does not find data.  I am attaching three files, one is my program and two are csv files with points and orientations. Once you have adjusted the path to the csv files in the python program, running_eds_extract_domain_model7.py, you should have not problems running it in jupyter.Once again, Thanks. Ed On Monday, June 7, 2021, 02:21:47 AM CDT, Jan von Harten ***@***.***> wrote: Hi, so what should not happen is that it does not find any minumum z locations for a layer (surface in Gempy), there might be gaps though, like in my example, where a layer does not exist for certain xy coordinates (in the example this is due to a unconformity, but might also be because of faulting, a pinch out, etc). Where did you upload your data, then I could maybe testrun it - didnt fiind it in your post. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe. import gempy as gp import numpy as np import pandas as pd import matplotlib.pyplot as plt def extract_domain(sol, unit): ''' Extract domain coordinates from gempy model by unit name arguments: sol: Gempy solution object. unit: id of gempy unit to be extracted returns: dom_x, dom_y, dom_z: coordinates of domain ''' # new version with rounding, definitely necessary rounded_lithblock = sol.lith_block.round(0) rounded_lithblock = rounded_lithblock.astype(int) # mask by array of input surfaces (by id, can be from different series) mask = np.isin(rounded_lithblock, unit) # get coordinates by mask # krig_lith = sol.lith_block[mask] dom_grid = sol.grid.values[mask] # dom_grid = sol.grid.values return dom_grid pd.set_option('precision', 2) data_path="C:\\Users\\17138\\geoLOGIC\\gridding\\test\\" path_to_data = data_path geo_data = gp.create_data('combination', extent=[755589.109, 780838.701, 7034865.1, 7053841.76, -1735.9, -1271.1], resolution=[125, 50, 50], path_o=path_to_data + "input_gempy1_orientations.csv", path_i=path_to_data + "input_gempy1.csv") geo_data.get_data() gp.map_stack_to_surfaces(geo_data, {"Strat_Series1": ('surface1'), "Basement_Series":('basement')}) gp.plot_2d(geo_data, direction='y') interp_data = gp.set_interpolator(geo_data, theano_optimizer='fast_compile') # %% sol = gp.compute_model(geo_data) # %% # Displaying the result in x and y direction: # # %% gp.plot_2d(geo_data, cell_number=20, direction='y', show_data=False, show_boundaries=True) grid = extract_domain(sol, 1) # extract single unit by id # Following part is to find lowest z value for each unique xy coordinate # very much hacked and not very nice, but hopefully does the job # make zeros array as start res = np.zeros((1,3)) # loop over unique x and y coordiantes and save minimum jcntr=0 for i in np.unique(grid[:,0]): for j in np.unique(grid[:,1]): mask = (grid[:,0]==i)&(grid[:,1]==j) temp1 = grid[mask] trc, trc1=temp1.shape if trc!= 0 and trc1 != 0 : temp2 = np.where(temp1[:,2] == np.amin(temp1[:,2])) temp3 = temp1[temp2] else: print("no value found for number: ",jcntr) jcntr=jcntr+ 1 # save result res = np.vstack((res, temp3)) # remove first zero entry res = res[1:,:] fig = plt.figure(figsize=(15, 10)) ax = fig.add_subplot(111) c = ax.scatter(res[:, 0], res[:, 1], c=res[:, 2], s=3, cmap="viridis") ax.set_aspect('equal') plt.colorbar(c, shrink=0.5) plt.show()

[emurphy222000-jpg]

Hi JanAttached is my latest version of the program as well as two csv files containing points and orientations. Again thanks so much for your great support.Ed On Monday, June 7, 2021, 02:21:47 AM CDT, Jan von Harten ***@***.***> wrote: Hi, so what should not happen is that it does not find any minumum z locations for a layer (surface in Gempy), there might be gaps though, like in my example, where a layer does not exist for certain xy coordinates (in the example this is due to a unconformity, but might also be because of faulting, a pinch out, etc). Where did you upload your data, then I could maybe testrun it - didnt fiind it in your post. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe. import gempy as gp import numpy as np import pandas as pd import matplotlib.pyplot as plt def extract_domain(sol, unit): ''' Extract domain coordinates from gempy model by unit name arguments: sol: Gempy solution object. unit: id of gempy unit to be extracted returns: dom_x, dom_y, dom_z: coordinates of domain ''' # new version with rounding, definitely necessary rounded_lithblock = sol.lith_block.round(0) rounded_lithblock = rounded_lithblock.astype(int) # mask by array of input surfaces (by id, can be from different series) mask = np.isin(rounded_lithblock, unit) # get coordinates by mask # krig_lith = sol.lith_block[mask] dom_grid = sol.grid.values[mask] # dom_grid = sol.grid.values return dom_grid pd.set_option('precision', 2) data_path="C:\\Users\\17138\\geoLOGIC\\gridding\\test\\" path_to_data = data_path geo_data = gp.create_data('combination', extent=[755589.109, 780838.701, 7034865.1, 7053841.76, -1735.9, -1271.1], resolution=[125, 50, 50], path_o=path_to_data + "input_gempy1_orientations.csv", path_i=path_to_data + "input_gempy1.csv") geo_data.get_data() gp.map_stack_to_surfaces(geo_data, {"Strat_Series1": ('surface1'), "Basement_Series":('basement')}) gp.plot_2d(geo_data, direction='y') interp_data = gp.set_interpolator(geo_data, theano_optimizer='fast_compile') # %% sol = gp.compute_model(geo_data) # %% # Displaying the result in x and y direction: # # %% gp.plot_2d(geo_data, cell_number=20, direction='y', show_data=False, show_boundaries=True) grid = extract_domain(sol, 1) # extract single unit by id # Following part is to find lowest z value for each unique xy coordinate # very much hacked and not very nice, but hopefully does the job # make zeros array as start res = np.zeros((1,3)) # loop over unique x and y coordiantes and save minimum jcntr=0 for i in np.unique(grid[:,0]): for j in np.unique(grid[:,1]): mask = (grid[:,0]==i)&(grid[:,1]==j) temp1 = grid[mask] trc, trc1=temp1.shape if trc!= 0 and trc1 != 0 : temp2 = np.where(temp1[:,2] == np.amin(temp1[:,2])) temp3 = temp1[temp2] else: print("no value found for number: ",jcntr) jcntr=jcntr+ 1 # save result res = np.vstack((res, temp3)) # remove first zero entry res = res[1:,:] fig = plt.figure(figsize=(15, 10)) ax = fig.add_subplot(111) c = ax.scatter(res[:, 0], res[:, 1], c=res[:, 2], s=3, cmap="viridis") ax.set_aspect('equal') plt.colorbar(c, shrink=0.5) plt.show()

[javoha]

Hi, sorry I am not sure if I am doing something wrong - but I cannot see any attachements on github - Not sure how exactly you are sending them at the moment.

[emurphy222000-jpg]

Oh, I am just sending emails independently of github. Could you give me the link to your github and I will try to put the files there. Thanks. On Monday, June 7, 2021, 07:22:19 AM CDT, Jan von Harten ***@***.***> wrote: Hi, sorry I am not sure if I am doing something wrong - but I cannot see any attachements on github - Not sure how exactly you are sending them at the moment. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[emurphy222000-jpg]

Finally figured out correct way to send files. I had to change the names to text files as it would not accept csv or py files.
Thanks so much.

input_gempy1_orientations.txt
input_gempy1.txt
running_eds_extract_domain_model7.txt

[emurphy222000-jpg]

I went to github and drop them there just a few moments ago. Please let me know if you can see them now. On Monday, June 7, 2021, 07:22:19 AM CDT, Jan von Harten ***@***.***> wrote: Hi, sorry I am not sure if I am doing something wrong - but I cannot see any attachements on github - Not sure how exactly you are sending them at the moment. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[emurphy222000-jpg]

running_eds_extract_domain_model7.txt
input_gempy1_orientations.txt
input_gempy1.txt

[emurphy222000-jpg]

Just tried again put them onto github. I went away and them came back and they still seem to be there. Please let me know. On Monday, June 7, 2021, 07:22:19 AM CDT, Jan von Harten ***@***.***> wrote: Hi, sorry I am not sure if I am doing something wrong - but I cannot see any attachements on github - Not sure how exactly you are sending them at the moment. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[javoha]

Hi @emurphy222000-jpg ,

finally I got around to run your example - for me it works fine. Here is the 3D model I get from your data points (30 times vertical exaggeration):

And here the result of my code, Basically there are no minimum Z vaues for the areas that exceed the model extend.

If you need values for this central part, maybe just increase the model extent in z direction.

Cheers, Jan

[emurphy222000-jpg]

Thanks Jan for your excellent support. Unfortunate it runs because I trap for missing values. Specifically  if trc!= 0 and trc1 != 0 : in the following code. If I remove this, it crashes. for i in np.unique(grid[:,0]):    for j in np.unique(grid[:,1]):        mask = (grid[:,0]==i)&(grid[:,1]==j)                temp1 = grid[mask]        trc, trc1=temp1.shape        if trc!= 0 and trc1 != 0 :            temp2 = np.where(temp1[:,2] == np.amin(temp1[:,2]))            temp3 = temp1[temp2]        else:            print("no value found for number: ",jcntr)            jcntr=jcntr+ 1        # save result        res = np.vstack((res, temp3)) I should have thought to mention that up front. Best Regards,Ed On Monday, June 14, 2021, 03:05:23 AM CDT, Jan von Harten ***@***.***> wrote: Hi @emurphy222000-jpg , finally I got around to run your example - for me it works fine. Here is the 3D model I get from your data points (30 times vertical exaggeration): And here the result of my code, Basically there are no minimum Z vaues for the areas that exceed the model extend. If you need values for this central part, maybe just increase the model extent in z direction. Cheers, Jan — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[javoha]

Ah, I see - but then you solved it, right? Or is this not what you need? Alternatively you can also set the upper z limit to -1000 to get a continiuous surface and a min z value for each point.

[emurphy222000-jpg]

I think the problem is the way I define extents. If it is too narrow it fails, as in this case. If I extend theme.g. extent=[700000, 800000, 5000000, 12000000, -2000, -1000],Then it runs without a problem. Is there a rule of thumb I should follow in defining the extents.Thanks,Ed On Monday, June 14, 2021, 07:38:53 AM CDT, Jan von Harten ***@***.***> wrote: Ah, I see - but then you solved it, right? Or is this not what you need? Alternatively you can also set the upper z limit to -1000 to get a continiuous surface and a min z value for each point. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[javoha]

There is not really a rule of thumb - for a typical geologic model it depends on the area you want to model. I your case you were cutting of the top of the strucutre, so if you need to capture this you should extend upwards. Hope it is coming together for you now.
Cheers
Jan

[emurphy222000-jpg]

Yes, your help is made a big difference. Thanks so much Jan. On Monday, June 21, 2021, 03:41:47 AM CDT, Jan von Harten ***@***.***> wrote: There is not really a rule of thumb - for a typical geologic model it depends on the area you want to model. I your case you were cutting of the top of the strucutre, so if you need to capture this you should extend upwards. Hope it is coming together for you now. Cheers Jan — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.