Model limits from shapefile?

[lucianomarquetto]

Hi community,

I am new to gempy and have little knowledge of python, but have been learning both in the past months. I wonder if it is possible to define my model extent from a shapefile or even just from a rotated rectangle. My area of interest is at an angle ~30° with north-south (see image below, yellow polygon).

I am using gempy 2.2.11. Any help would be welcome - sorry if I need to specify more information here, it is the first time I post a question on the discussions. Thanks!

[javoha]

Hi @lucianomarquetto,
great to see that you are looking into gempy. To my knowledge the extent will always be rectangular in the current state of gempy. Therefore I am not sure if you can use exactly your shapefile, as it looks like the north-eastern boundary is a little shorter than the south-western one. You can however use the coordinates of your corner points to create an extent that should closely resemble your shape, for example like this:

extent = np.array([729550.0, 751500.0, 1913500.0, 1923650.0, -1800.0, 800.0])

[lucianomarquetto]

Hi Jan (@javoha)! Thank you for your response. I am enjoying learning Gempy and hope to keep doing so - it seems a very powerful tool!

Based on your response I think I might have misunderstood what regular and custom grids mean. Reading through tutorial "1.3a: Grids" I thought a custom grid would allow any polygon shape (from that tutorial: "Custom grid: Completely free XYZ values. Parameters
custom_grid (np.array) – [s0] (numpy.ndarray[float, 3]): XYZ 2D array. Axis 1 is the coordinates while axis 0 is n number of input)".

I even found another question in #604 that seems to be similar to mine in a way, and it was responded suggesting the use of a custom grid.

So I followed the steps similar to the code in that thread (my code below). The problem is, even though I do not get any error and have a message of "Active grids: ['custom']" below the "set_custom_grid", I cannot plot in 2d or compute the model. When I try to plot in 2d I get "IndexError: index 0 is out of bounds for axis 0 with size 0", and when I compute the model I get an empty array: Lithology ids [ ].

---

• Custom grid extent I'm expecting to get
  AOI_limits_custom = (np.array([[598966,6774879, -200], [628406,6755377, -200], [577976,6655121, 50], [548422,6674680, 50]]))

• Defining the grid
  grid = gp.Grid()

• Setting the custom grid
  grid = geo_model.set_custom_grid(AOI_limits_custom)

Active grids: ['custom']

• Importing data
  geo_model = gp.init_data(geo_model, path_o=data_path + '/Gempy_files/Gempy_input_orientations_troian_2020.csv', path_i=data_path + '/Gempy_files/Gempy_input_surfaces_troian_2020.csv')

• Plot points in 2D - NOT WORKING
  fig_2D = gp.plot_2d(geo_model, direction=['z'])

IndexError: index 0 is out of bounds for axis 0 with size 0

• Plot points in 3D - WORKING
  gp.plot_3d(geo_model, ve=100, plotter_type='background', show_topography=False)

• Setting the interpolator so ir can compute the model
  gp.set_interpolator(geo_model, theano_optimizer='fast_compile', verbose=[])

• Compute model - NOT WORKING
  gp.compute_model(geo_model)

Lithology ids
[ ]

---

Interestingly, when I try to plot in 3D, though, my points appear, but the custom grid extent does not seem to be correct (e.g. the northern corner of my grid should be in 598966, 6774879 but the axis ticks do not get there). Anyway, any help would be welcome!

Cheers!

[javoha]

Ok I think I see the problem, I'll try to go step by step and hope it's helpful. Maybe first of all - when you plot your model without a successful computation (or before you acutally run the model computation) gempy will default to just showing the input data with an extent that fits it - that is why you get the 3D plot.
To my knowledge a custom grid is different to what you want it to be: In gempy, a regular grid, meaning a grid where all points are ordered in a rectangular fashion, will automatically be created based on your extend and resolution. Many other applications use different meshes, where grid points are ordered in a different fashion. You can use the custom grid in gempy to use such meshes by basically passing the coordinates of all the grid points - instead of only the corner points. I think gempy will still display the result in a rectangular extent though, so I am pretty sure that a custom grid is not a solution to your problem.
Is there a reason why you need to fit the extent exactly to the shapefile? You could probably extract the surfaces calculated by gempy and cut the results to the exact extend, but that probably requires a couple of extra steps.

[lucianomarquetto]

Hi again Jan, and thanks for clarifying this issue for me - I think I understand now what the regular and custom grids are. My aim was to focus on the yellow rectangle because I am working with the top sequence of the coastal, sedimentary aquifer (green color in the first image I sent - around 150 m thick). To the northwest in the image there is a sequence of volcanic rocks from -400 up to +1000 m altitude, and to the southeast is the sea, with a different geological setting as well. Neither the marine environment nor the volcanic rocks are my target. I thought that reducing my area of interest to the yellow rectangle would be easier to populate with borehole data and also less demanding to compute the model (I am working on my notebook).

Anyway, I think I can mark this question as answered as you gave me the response to my original question. Still, If there are other ways of focusing on that extent, I would be glad if you could point out how to do that - or at least mention the function names I need so I can have a look at them.

Cheers!

Luciano.

[javoha]

You are absolutely right that it makes sense to limit the extent as much as possible to be able to create a higher resolution model. Maybe I was not super clear in my first answer: It's absolutely possible to limit yourself to the yellow rectangle - to me it just looks as if it's not exactly rectangular. If it is - just use the corner coordinates as an extend as in the example I gave above (actually taken from this example.

[lucianomarquetto]

Haha sorry I just got more confused again. Let me just ask again to be sure I understood it, then I will close the issue.

When you say "It's absolutely possible to limit yourself to the yellow rectangle" to me it seems I can have my model extent as four different [x,y] coordinates, as for example [598966, 6774880], [628406, 6755377], [577976, 6655115], [548426, 6674680] - yellow rectangle in image below, no [x] or [y] coordinates repeat in any of the vertices.

But when you refer to the Alesmodel example, the four coordinates will have repeating [x] or [y] values and thus I will have something as the red rectangle in the image below, with boundary lines being E-W or N-S, in no other azimuth.

By our discussion so far I assume that the best I can do to limit my extent is define it as the red rectangle, is that correct?

Again, sorry to take your time, and thank you for helping!

[javoha]

Ah yes, you are right. I feel like I might be missing an easier solution here (maybe @AlexanderJuestel knows), but you could just apply a 2D rotation (like this) to the x,y of your extent points and data points to rotate the whole thing in a way that lets you compute the model only in the yellow area. You would need to retransform the result if you need the real coordinates for something later though.

[lucianomarquetto]

Ahh, nice Jan, I think this solution might be what I am looking for! As I mentioned before I am still crawling in both python and gempy, so everything is new. Thanks for your patience and time helping me!

[lucianomarquetto]

Hi all!

I just thought about coming back here to share the code I used to rotate my model as Jan suggested in the comments, in case anyone like me (a beginner in python and gempy) needs it. As Jan mentioned, the coordinates of the rotated model need to be transformed back in case one needs the real coordinates for anything else, but other than that the model seems ok. Any suggestions on how to improve the code are also welcome!

---

'''Defining the extent of the model and creating a rotation matrix in order to model only the necessary area:'''

#Define the extent - in my case it matches with the yellow rectangle drawn in the images I posted earlier in this discussion:
extent = np.array([[598966, 6774880, -200], [628406, 6755377, -200], [577976, 6655115, 50], [548426, 6674680, 50]])

#Define the rotation angle in degrees
angle_deg = -25

#Convert the angle to radians
angle_rad = np.radians(angle_deg)

#Define the rotation matrix
cos = np.cos(angle_rad) sin = np.sin(angle_rad) rotation_matrix = np.array([[cos, -sin, 0], [sin, cos, 0], [0, 0, 1]])

#Apply the rotation matrix to the extent
rotated_extent = np.dot(extent, rotation_matrix)

#Extract the desired extent values
xmin = np.min(rotated_extent[:, 0]) xmax = np.max(rotated_extent[:, 0]) ymin = np.min(rotated_extent[:, 1]) ymax = np.max(rotated_extent[:, 1]) zmin = np.min(rotated_extent[:, 2]) zmax = np.max(rotated_extent[:, 2])

extent = [xmin, xmax, ymin, ymax, zmin, zmax]

'''Creating a rotation matrix for the dataset, and writing the new coordinates (rotated X and rotated Y) to the csv file. Note that here the original coordinates are the columns 'X_original' and 'Y_original', while the rotated coordinates are 'X' and 'Y' because Gempy looks for the latter when plotting the points and computing the model.'''

#Read in the CSV file with data points
df = pd.read_csv(data_path + '/Gempy_files/Gempy_input_orientations_troian_2020.csv')

#Extract the x and y coordinates as separate arrays
x_coords = df['X_original'].to_numpy() y_coords = df['Y_original'].to_numpy()

# Define the rotation matrix
rotation_matrix_points = np.array([[cos, -sin], [sin, cos]])

# Apply the rotation matrix to the x and y columns of the dataframe
rotated_points = np.dot(df[['X_original', 'Y_original']], rotation_matrix_points) df['X'] = rotated_points[:, 0] df['Y'] = rotated_points[:, 1]

#Apply the rotation matrix to the x and y columns of the dataframe
rotated_points = np.dot(df[['X_original', 'Y_original']], rotation_matrix) df['X'] = rotated_points[:, 0] df['Y'] = rotated_points[:, 1]

#Save the modified dataframe to the original csv file
df.to_csv(data_path + '/Gempy_files/Gempy_input_orientations_troian_2020.csv', index=False)

#After that can create and initialize the model.`

[javoha]

Thanks a lot for sharing - that is exactly what an open source community like ours needs. I am absolutely sure that people will find it useful!

[lucianomarquetto]

Hi all!

I am opening this discussion again instead of creating a new post as I think this is the continuation of my issue. This might be a tough question (at least it seems for me that do not understand a lot of programming).

I was able to rotate my model as suggested by @javoha , by applying a 2D matrix rotation in the model extent + surface&orientation points. Results are awesome, and this definitely helped me work only on my area of interest (I was able to focus and detail on the 3600 km² of my area instead of modeling ~9600 km²).

Nonetheless, now I want to rotate my model back, but I am not sure if just rotating back the coordinates in the the output files is enough. My model is already useful as it is, but would be great to have the real coordinates so I (or other people) could use it further.

I notice the 'gp.save_model' generates 11 files, from which 6 show coordinates (2 .npy files and 4 .csv files). Those are listed below (considering my model name as 'Final_model'):

Final_model_extent.npy
Final_model_resolution.npy
Final_model_orientations.csv
Final_model_rescaling_data.csv
Final_model_surface_points.csv
Final_model_surfaces.csv

Final_model_extent.npy and Final_model_resolution.npy are easy to replace with the original coordinates, as well as Final_model_orientations.csv and Final_model_surface_points.csv. (although those last two do contain more than just coordinates, i.e. 'X_c', 'Y_c' columns that I don't know if I need to adjust somehow).

Then we get to the Final_model_rescaling_data.csv and Final_model_surfaces.csv, which seem harder to work with. Final_model_rescaling_data.csv has two columns: 'rescaling factor' and 'centers', the latter is a field of coordinates.
Final_model_surfaces.csv has 'vertices', 'edges' and 'sfai' columns, from which only 'vertices' I think are related to coordinates.

So my main question here is, how easy would be to rotate my model back again? Do I need to worry only with coordinates or there are other fields I have to convert (somehow) as well?

Thanks!

Luciano.

[AlexanderJuestel]

@lucianomarquetto I was going over this discussion again and would like to chip in. What information do you need from the resulting GemPy model? I would not convert everything but just the results that you need. If you need the meshes, you could try to rotate them (mind the center of rotation though). I would like to implement this rotation feature in GemGIS as it may be useful to many other GemPy users and would like to provide an example for it. Maybe we could use your case study to illustrate the feature?

[AlexanderJuestel]

@lucianomarquetto I was reimplementing your code and also re-run your code and found out that the resulting rectangle (red) is quite "misplaced" compared to the original one (blue) (see image). I guess this is due to the orientation around the origin and the large coordinates that your box has.

Using the following shapely function where the polygon is rotated around its center will return the following which looks more applicable for not-so-advanced users of GemPy. We could then just extract the total_bounds from the Polygon and pass it as extent to GemPy :)

from shapely import affinity
extent_rotated_shapely = affinity.rotate(Polygon(extent), 25, 'center')

I will now to come up with some code to rotate potential interface points :)

I will keep you posted :)

[lucianomarquetto]

Hi @AlexanderJuestel!

First of all thanks for your interest in this issue - of course you can use it as a case study, I am delighted your team got interested in it! I am still writing the manuscript related to this work, but once it is published (hope within a year or so) I will let you know so the reference can be added to the text too!

Regarding the rotation I used - well, as you could notice, I am a beginner in gempy and programming in general, so I was following what I could find and that worked. Although the coordinates are really off, it seems the area kept its proportions, which is what I was looking for in the beginning. I used UTM coordinates for my area - maybe should have used decimal degrees? Also, I believe this rotation only affected x and y, right? I ask because I also made some depth maps with gemgis that I was considering to be accurate, but I guess it is better to ask before adding to the manuscript.

For the reverse rotation I was looking to get all interpolated surfaces in their original coordinates, so I could work with them in a GIS base. Basically I wanted the real x, y and z.

Anyway, would be great to see your code related to this feature once it is done - It will be great to use it in my model, and I am sure it will be very useful for other Gempy users too! For now I will try changing my code to the one you provided and see how it goes!

Thanks and cheers Alexander!

PS: Also congratulations to the team for the latest release of Gempy! I did not have much time to use it yet, but saw the doc with main changes. Thanks for your hard work providing everyone around the world with a robust-and-free way of modelling!

[AlexanderJuestel]

@lucianomarquetto, Just playing with some example data.

I wrote a function using Shapely Methods to perform the rotation of the extent, the interfaces, and the orientations. Works like a charm so far :)

Input for the extent is either a NumPy Array with four corner coordinates, a Shapely Polygon or a GeoDataFrame/Shapefile containing the bounding box. The angle to rotate the box will be calculated automatically to rotate the bounding box horizontally or vertically. A new bounding box will be created from the rotated box to get the total bounds of the box and will be returned as new extent for the modeling in GemPy.

[lucianomarquetto]

Haha, awesome @AlexanderJuestel!

I just responded the other message - wasn't expecting an update that soon! This makes me want to start modelling other areas as soon as I can!

Thanks so far!

Luciano

[AlexanderJuestel]

@lucianomarquetto

xref: #776 (reply in thread)

Just exported some data from another example model. With PyVista, you can easily rotate meshes using the following code (rotate_z rotates around the z-axis). If you use the rotation-around-the-center-method that I have implemented, this should work as a reversal of the original rotation. Afterwards, you should be able to export the mesh as tiff/zmap for further usage in QGIS/ArcGIS...

angle = 45 # angle the extent was rotated
origin = (500,500,0) # point around which the extent was rotated

mesh_rotated = mesh.rotate_z(angle, point=origin, inplace=False)

The greenish mesh is the original one, the mesh with different colors the rotated one.

[AlexanderJuestel]

@lucianomarquetto, I have created a new tutorial for GemGIS that outlines the functionality of how to rotate GemPy Input Data: https://gemgis.readthedocs.io/en/latest/getting_started/tutorial/67_rotating_gempy_input_data.html

The functionality (rotate_gempy_input_data(..)) will be added to the main branch of GemGIS via cgre-aachen/gemgis#286 but is already available in the current development branch.