diff --git a/.gitignore b/.gitignore index 2a0cb4779..c579dfb5c 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,3 @@ - - .vscode .idea .DS_Store @@ -172,3 +170,14 @@ examples/integrations/*.out /.obsidian/core-plugins-migration.json /.obsidian/hotkeys.json /.obsidian/workspace.json + +# Ignore generated tutorial files +examples/tutorials/z_other_tutorials/json_io/04_simple_layer_stack.py +examples/tutorials/z_other_tutorials/json_io/fault_scalar_field.png +examples/tutorials/z_other_tutorials/json_io/initial_model_y.png +examples/tutorials/z_other_tutorials/json_io/model_with_data.png +examples/tutorials/z_other_tutorials/json_io/multiple_series_faults_computed.json + +# Generated JSON files from examples +examples/tutorials/z_other_tutorials/json_io/combination_model.json +examples/tutorials/z_other_tutorials/json_io/combination_model_computed.json diff --git a/examples/tutorials/z_other_tutorials/json_io/01_surface_points_io.py b/examples/tutorials/z_other_tutorials/json_io/01_surface_points_io.py new file mode 100644 index 000000000..9bf1a7423 --- /dev/null +++ b/examples/tutorials/z_other_tutorials/json_io/01_surface_points_io.py @@ -0,0 +1,95 @@ +""" +Tutorial for JSON I/O operations in GemPy - Surface Points +======================================================= + +This tutorial demonstrates how to save and load surface points data using JSON files in GemPy. +""" + +import json +import numpy as np +import gempy as gp +from gempy.modules.json_io import JsonIO + +# Create a sample surface points dataset +# ------------------------------------ + +# Create some sample surface points +x = np.array([0, 1, 2, 3, 4]) +y = np.array([0, 1, 2, 3, 4]) +z = np.array([0, 1, 2, 3, 4]) +ids = np.array([0, 0, 1, 1, 2]) # Three different surfaces +nugget = np.array([0.00002, 0.00002, 0.00002, 0.00002, 0.00002]) + +# Create name to id mapping +name_id_map = {f"surface_{id}": id for id in np.unique(ids)} + +# Create a SurfacePointsTable +surface_points = gp.data.SurfacePointsTable.from_arrays( + x=x, + y=y, + z=z, + names=[f"surface_{id}" for id in ids], + nugget=nugget, + name_id_map=name_id_map +) + +# Create a JSON file with the surface points data +# --------------------------------------------- + +# Create the JSON structure +json_data = { + "metadata": { + "name": "sample_model", + "creation_date": "2024-03-19", + "last_modification_date": "2024-03-19", + "owner": "tutorial" + }, + "surface_points": [ + { + "x": float(x[i]), + "y": float(y[i]), + "z": float(z[i]), + "id": int(ids[i]), + "nugget": float(nugget[i]) + } + for i in range(len(x)) + ], + "orientations": [], + "faults": [], + "series": [], + "grid_settings": { + "regular_grid_resolution": [10, 10, 10], + "regular_grid_extent": [0, 4, 0, 4, 0, 4], + "octree_levels": None + }, + "interpolation_options": {} +} + +# Save the JSON file +with open("sample_surface_points.json", "w") as f: + json.dump(json_data, f, indent=4) + +# Load the surface points from JSON +# ------------------------------- + +# Load the model from JSON +loaded_surface_points = JsonIO._load_surface_points(json_data["surface_points"]) + +# Verify the loaded data +print("\nOriginal surface points:") +print(surface_points) +print("\nLoaded surface points:") +print(loaded_surface_points) + +# Verify the data matches +print("\nVerifying data matches:") +print(f"X coordinates match: {np.allclose(surface_points.xyz[:, 0], loaded_surface_points.xyz[:, 0])}") +print(f"Y coordinates match: {np.allclose(surface_points.xyz[:, 1], loaded_surface_points.xyz[:, 1])}") +print(f"Z coordinates match: {np.allclose(surface_points.xyz[:, 2], loaded_surface_points.xyz[:, 2])}") +print(f"IDs match: {np.array_equal(surface_points.ids, loaded_surface_points.ids)}") +print(f"Nugget values match: {np.allclose(surface_points.nugget, loaded_surface_points.nugget)}") + +# Print the name_id_maps to compare +print("\nName to ID mappings:") +print(f"Original: {surface_points.name_id_map}") +print(f"Loaded: {loaded_surface_points.name_id_map}") \ No newline at end of file diff --git a/examples/tutorials/z_other_tutorials/json_io/02_horizontal_stratigraphic.py b/examples/tutorials/z_other_tutorials/json_io/02_horizontal_stratigraphic.py new file mode 100644 index 000000000..f96197b5c --- /dev/null +++ b/examples/tutorials/z_other_tutorials/json_io/02_horizontal_stratigraphic.py @@ -0,0 +1,106 @@ +""" +Tutorial: Loading a Horizontal Stratigraphic Model using JSON I/O +=============================================================== + +This tutorial demonstrates how to load a horizontal stratigraphic model using GemPy's JSON I/O functionality. +The model consists of two horizontal layers (rock1 and rock2) with surface points and orientations. +""" + +# %% +# Import necessary libraries +import matplotlib +matplotlib.use('Agg') # Use non-interactive backend +import gempy as gp +import gempy_viewer as gpv +import numpy as np +import json +from pathlib import Path +import matplotlib.pyplot as plt +from datetime import datetime + +# %% +# Define the model data +model_data = { + "metadata": { + "name": "Horizontal Stratigraphic Model", + "creation_date": "2024-03-24", + "last_modification_date": datetime.now().strftime("%Y-%m-%d"), + "owner": "GemPy Team" + }, + "surface_points": [ + {"x": 0.0, "y": 0.0, "z": 0.0, "id": 1, "nugget": 0.0}, + {"x": 1.0, "y": 0.0, "z": 0.0, "id": 1, "nugget": 0.0}, + {"x": 0.0, "y": 1.0, "z": 0.0, "id": 1, "nugget": 0.0}, + {"x": 1.0, "y": 1.0, "z": 0.0, "id": 1, "nugget": 0.0}, + {"x": 0.0, "y": 0.0, "z": 1.0, "id": 2, "nugget": 0.0}, + {"x": 1.0, "y": 0.0, "z": 1.0, "id": 2, "nugget": 0.0}, + {"x": 0.0, "y": 1.0, "z": 1.0, "id": 2, "nugget": 0.0}, + {"x": 1.0, "y": 1.0, "z": 1.0, "id": 2, "nugget": 0.0} + ], + "orientations": [ + {"x": 0.5, "y": 0.5, "z": 0.0, "G_x": 0.0, "G_y": 0.0, "G_z": 1.0, "id": 1, "nugget": 0.0, "polarity": 1}, + {"x": 0.5, "y": 0.5, "z": 1.0, "G_x": 0.0, "G_y": 0.0, "G_z": 1.0, "id": 2, "nugget": 0.0, "polarity": 1} + ], + "series": [ + { + "name": "series1", + "surfaces": ["layer1", "layer2"], + "structural_relation": "ERODE", + "colors": ["#ff0000", "#00ff00"] + } + ], + "grid_settings": { + "regular_grid_resolution": [10, 10, 10], + "regular_grid_extent": [0.0, 1.0, 0.0, 1.0, 0.0, 1.0], + "octree_levels": None + }, + "interpolation_options": { + "kernel_options": { + "range": 1.7, + "c_o": 10 + }, + "mesh_extraction": True, + "number_octree_levels": 1 + }, + "fault_relations": None, + "id_name_mapping": { + "name_to_id": { + "layer1": 1, + "layer2": 2 + } + } +} + +# %% +# Save the model data to a JSON file +tutorial_dir = Path(__file__).parent +json_file = tutorial_dir / "horizontal_stratigraphic.json" +with open(json_file, "w") as f: + json.dump(model_data, f, indent=4) + +# %% +# Load the model from JSON +model = gp.modules.json_io.JsonIO.load_model_from_json(str(json_file)) + +# %% +# Compute the geological model +gp.compute_model(model) + +# %% +# Plot the model +# Plot the initial geological model in the y direction without results +fig, ax = plt.subplots(figsize=(10, 6)) +gpv.plot_2d(model, direction=['y'], show_results=False, ax=ax) +plt.title("Initial Geological Model (y direction)") +plt.savefig('initial_model_y.png') +plt.close() + +# Plot the result of the model in the x and y direction with data and without boundaries +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6)) +gpv.plot_2d(model, direction=['x'], show_data=True, show_boundaries=False, ax=ax1) +ax1.set_title("Model with Data (x direction)") +gpv.plot_2d(model, direction=['y'], show_data=True, show_boundaries=False, ax=ax2) +ax2.set_title("Model with Data (y direction)") +plt.tight_layout() +plt.savefig('model_with_data.png') +plt.close() \ No newline at end of file diff --git a/examples/tutorials/z_other_tutorials/json_io/03_multiple_series_faults.py b/examples/tutorials/z_other_tutorials/json_io/03_multiple_series_faults.py new file mode 100644 index 000000000..f08b29777 --- /dev/null +++ b/examples/tutorials/z_other_tutorials/json_io/03_multiple_series_faults.py @@ -0,0 +1,252 @@ +""" +Tutorial: Loading a Model with Multiple Series and Faults using JSON I/O +===================================================================== + +This tutorial demonstrates how to load a geological model with multiple series and faults using GemPy's JSON I/O functionality. +The model consists of two layers (rock1, rock2) and a fault that offsets them. +""" + +# %% +# Import necessary libraries +import matplotlib +# matplotlib.use('Agg') # Use non-interactive backend +import gempy as gp +import gempy_viewer as gpv +import numpy as np +import json +from pathlib import Path +import matplotlib.pyplot as plt +from gempy_engine.core.data.stack_relation_type import StackRelationType +from gempy.modules.json_io.json_operations import JsonIO # Updated import path +from datetime import datetime + +# %% +# Define the model data +model_data = { + "metadata": { + "name": "Multiple Series and Faults Model", + "creation_date": "2024-03-24", # Fixed creation date + "last_modification_date": datetime.now().strftime("%Y-%m-%d"), # Current date only + "owner": "GemPy Team" + }, + "surface_points": [ + # rock1 points + {"x": 0, "y": 200, "z": 600, "id": 0, "nugget": 0.00002}, + {"x": 0, "y": 500, "z": 600, "id": 0, "nugget": 0.00002}, + {"x": 0, "y": 800, "z": 600, "id": 0, "nugget": 0.00002}, + {"x": 200, "y": 200, "z": 600, "id": 0, "nugget": 0.00002}, + {"x": 200, "y": 500, "z": 600, "id": 0, "nugget": 0.00002}, + {"x": 200, "y": 800, "z": 600, "id": 0, "nugget": 0.00002}, + {"x": 800, "y": 200, "z": 200, "id": 0, "nugget": 0.00002}, + {"x": 800, "y": 500, "z": 200, "id": 0, "nugget": 0.00002}, + {"x": 800, "y": 800, "z": 200, "id": 0, "nugget": 0.00002}, + {"x": 1000, "y": 200, "z": 200, "id": 0, "nugget": 0.00002}, + {"x": 1000, "y": 500, "z": 200, "id": 0, "nugget": 0.00002}, + {"x": 1000, "y": 800, "z": 200, "id": 0, "nugget": 0.00002}, + # rock2 points + {"x": 0, "y": 200, "z": 800, "id": 1, "nugget": 0.00002}, + {"x": 0, "y": 800, "z": 800, "id": 1, "nugget": 0.00002}, + {"x": 200, "y": 200, "z": 800, "id": 1, "nugget": 0.00002}, + {"x": 200, "y": 800, "z": 800, "id": 1, "nugget": 0.00002}, + {"x": 800, "y": 200, "z": 400, "id": 1, "nugget": 0.00002}, + {"x": 800, "y": 800, "z": 400, "id": 1, "nugget": 0.00002}, + {"x": 1000, "y": 200, "z": 400, "id": 1, "nugget": 0.00002}, + {"x": 1000, "y": 800, "z": 400, "id": 1, "nugget": 0.00002}, + # fault points + {"x": 500, "y": 500, "z": 500, "id": 2, "nugget": 0.00002}, + {"x": 450, "y": 500, "z": 600, "id": 2, "nugget": 0.00002}, + {"x": 500, "y": 200, "z": 500, "id": 2, "nugget": 0.00002}, + {"x": 450, "y": 200, "z": 600, "id": 2, "nugget": 0.00002}, + {"x": 500, "y": 800, "z": 500, "id": 2, "nugget": 0.00002}, + {"x": 450, "y": 800, "z": 600, "id": 2, "nugget": 0.00002} + ], + "orientations": [ + # rock2 orientations + {"x": 100, "y": 500, "z": 800, "G_x": 0, "G_y": 0, "G_z": 1, "id": 1, "nugget": 0.00002, "polarity": 1}, + {"x": 900, "y": 500, "z": 400, "G_x": 0, "G_y": 0, "G_z": 1, "id": 1, "nugget": 0.00002, "polarity": 1}, + # rock1 orientations + {"x": 100, "y": 500, "z": 600, "G_x": 0, "G_y": 0, "G_z": 1, "id": 0, "nugget": 0.00002, "polarity": 1}, + {"x": 900, "y": 500, "z": 200, "G_x": 0, "G_y": 0, "G_z": 1, "id": 0, "nugget": 0.00002, "polarity": 1}, + # fault orientation (60-degree dip) + {"x": 500, "y": 500, "z": 500, "G_x": 0.866, "G_y": 0, "G_z": 0.5, "id": 2, "nugget": 0.00002, "polarity": 1} + ], + "series": [ + { + "name": "Fault_Series", + "surfaces": ["fault"], + "structural_relation": "FAULT", + "colors": ["#ffbe00"] + }, + { + "name": "Strat_Series", + "surfaces": ["rock2", "rock1"], + "structural_relation": "ERODE", + "colors": ["#015482", "#9f0052"] + } + ], + "grid_settings": { + "regular_grid_resolution": [50, 50, 50], + "regular_grid_extent": [0, 1000, 0, 1000, 0, 1000], + "octree_levels": None + }, + "interpolation_options": { + "kernel_options": { + "range": 1.7, + "c_o": 10 + }, + "mesh_extraction": True, + "number_octree_levels": 1 + }, + "fault_relations": [[0, 1], [0, 0]], # Fault series affects Strat_Series + "id_name_mapping": { + "name_to_id": { + "rock1": 0, + "rock2": 1, + "fault": 2 + } + } +} + +# %% +# Save the model data to a JSON file +tutorial_dir = Path(__file__).parent +json_file = tutorial_dir / "multiple_series_faults.json" +with open(json_file, "w") as f: + json.dump(model_data, f, indent=4) + +# %% +# Load the model from JSON +model = JsonIO.load_model_from_json(str(json_file)) + +# Print metadata to verify it's properly loaded +print("\nModel Metadata:") +print(f"Name: {model.meta.name}") +print(f"Creation Date: {model.meta.creation_date}") +print(f"Last Modification Date: {model.meta.last_modification_date}") +print(f"Owner: {model.meta.owner}") + +# Print structural groups +print("\nStructural Groups:") +print(model.structural_frame.structural_groups) + +# %% +# Compute the geological model +gp.compute_model(model) + +# Print scalar field values to verify they are calculated +print("\nScalar Field Values (Fault Series):") +print(f"Shape: {model.solutions.raw_arrays.scalar_field_matrix.shape}") +print(f"Min value: {model.solutions.raw_arrays.scalar_field_matrix[0].min()}") +print(f"Max value: {model.solutions.raw_arrays.scalar_field_matrix[0].max()}") +print(f"Mean value: {model.solutions.raw_arrays.scalar_field_matrix[0].mean()}") + +# %% +# Save the computed model to a new JSON file +computed_json_file = tutorial_dir / "multiple_series_faults_computed.json" +JsonIO.save_model_to_json(model, str(computed_json_file)) +print(f"\nSaved computed model to: {computed_json_file}") + +# %% +# Load the computed model back to verify metadata is preserved +reloaded_model = JsonIO.load_model_from_json(str(computed_json_file)) +print("\nReloaded Model Metadata:") +print(f"Name: {reloaded_model.meta.name}") +print(f"Creation Date: {reloaded_model.meta.creation_date}") +print(f"Last Modification Date: {reloaded_model.meta.last_modification_date}") +print(f"Owner: {reloaded_model.meta.owner}") + +# %% +# Create plots with proper configuration +# Plot 1: Cross-section in Y direction (XZ plane) +fig = plt.figure(figsize=(10, 8)) +ax = fig.add_subplot(111) +gpv.plot_2d( + model, + cell_number=25, # Middle of the model + direction='y', + show_data=True, + show_boundaries=True, + show_results=True, + ax=ax +) +plt.title("Geological Model - Y Direction (XZ plane)") +plt.savefig('model_y_direction.png', dpi=300, bbox_inches='tight') +plt.close() + +# Plot 2: Cross-section in X direction (YZ plane) +fig = plt.figure(figsize=(10, 8)) +ax = fig.add_subplot(111) +gpv.plot_2d( + model, + cell_number=25, # Middle of the model + direction='x', + show_data=True, + show_boundaries=True, + show_results=True, + ax=ax +) +plt.title("Geological Model - X Direction (YZ plane)") +plt.savefig('model_x_direction.png', dpi=300, bbox_inches='tight') +plt.close() + +# Plot 3: Scalar field of the fault +plt.figure(figsize=(10, 8)) +ax = plt.gca() + +# Get scalar field values and reshape to grid +scalar_field = model.solutions.raw_arrays.scalar_field_matrix[0].reshape(50, 50, 50) + +# Plot middle slice in Y direction +middle_slice = scalar_field[:, 25, :] +im = ax.imshow(middle_slice.T, + extent=[0, 1000, 0, 1000], + origin='lower', + cmap='RdBu', + aspect='equal') + +# Add colorbar +plt.colorbar(im, ax=ax, label='Scalar Field Value') + +# Plot surface points +fault_element = model.structural_frame.get_element_by_name("fault") +if fault_element and fault_element.surface_points is not None: + fault_points_coords = fault_element.surface_points.xyz + + # Filter points near the slice (Y around 500) + mask = np.abs(fault_points_coords[:, 1] - 500) < 100 + filtered_points = fault_points_coords[mask] + + if len(filtered_points) > 0: + ax.scatter(filtered_points[:, 0], filtered_points[:, 2], + c='red', s=50, label='Surface Points') + ax.legend() + +ax.set_xlabel('X') +ax.set_ylabel('Z') +ax.set_title('Fault Scalar Field - Y Direction (Middle Slice)') + +# Save plot +plt.savefig('fault_scalar_field.png', dpi=300, bbox_inches='tight') +plt.close() + +print("\nPlot saved as fault_scalar_field.png") + +# Plot 4: 3D visualization (optional) +PLOT_3D = False # Set to True to enable 3D plotting + +if PLOT_3D: + try: + import pyvista as pv + p = pv.Plotter(notebook=False, off_screen=True) + gpv.plot_3d( + model, + show_data=True, + show_surfaces=True, + show_boundaries=True, + plotter=p + ) + p.screenshot('model_3d.png', transparent_background=False) + p.close() + except Exception as e: + print(f"Could not create 3D plot: {e}") +# %% diff --git a/examples/tutorials/z_other_tutorials/json_io/04_combination_model.py b/examples/tutorials/z_other_tutorials/json_io/04_combination_model.py new file mode 100644 index 000000000..def9c3417 --- /dev/null +++ b/examples/tutorials/z_other_tutorials/json_io/04_combination_model.py @@ -0,0 +1,134 @@ +# %% +""" +Combination Model with JSON I/O +============================== + +This example demonstrates how to create a model combining faults and unconformities using GemPy's JSON I/O functionality. +The model is based on the g07_combination.py example, featuring a folded domain with an unconformity and a fault. + +Part 1: Create and save the initial model structure +Part 2: Load the model, compute it, and visualize results +""" + +import gempy as gp +import gempy_viewer as gpv +import numpy as np +import os +from gempy_engine.core.data.stack_relation_type import StackRelationType +from gempy_engine.config import AvailableBackends +from gempy.modules.json_io.json_operations import JsonIO + +# %% +# Part 1: Create and save the initial model structure +# ------------------------------------------------- + +# Define paths +data_path = 'https://raw.githubusercontent.com/cgre-aachen/gempy_data/master/' +path_to_data = data_path + "/data/input_data/jan_models/" +json_file = "combination_model.json" +computed_json_file = "combination_model_computed.json" + +# Create the model with data import +model = gp.create_geomodel( + project_name='Combination Model', + extent=[0, 2500, 0, 1000, 0, 1000], + resolution=[20, 20, 20], + refinement=6, + importer_helper=gp.data.ImporterHelper( + path_to_orientations=path_to_data + "model7_orientations.csv", + path_to_surface_points=path_to_data + "model7_surface_points.csv" + ) +) + +# Set metadata +model.meta.creation_date = "2024-03-24" +model.meta.last_modification_date = "2024-03-24" + +# Map geological series to surfaces +gp.map_stack_to_surfaces( + gempy_model=model, + mapping_object={ + "Fault_Series": ('fault',), + "Strat_Series1": ('rock3',), + "Strat_Series2": ('rock2', 'rock1'), + } +) + +# Set structural relations +model.structural_frame.structural_groups[0].structural_relation = StackRelationType.FAULT +model.structural_frame.fault_relations = np.array([ + [0, 1, 1], + [0, 0, 0], + [0, 0, 0] +]) + +# Set colors for visualization +model.structural_frame.get_element_by_name("fault").color = "#015482" +model.structural_frame.get_element_by_name("rock3").color = "#9f0052" +model.structural_frame.get_element_by_name("rock2").color = "#ffbe00" +model.structural_frame.get_element_by_name("rock1").color = "#728f02" + +# Set interpolation options +model.interpolation_options.number_octree_levels_surface = 5 + +# Save the model data to a JSON file +JsonIO.save_model_to_json(model, json_file) +print(f"\nSaved initial model to: {os.path.abspath(json_file)}") + +# %% +# Part 2: Load the model and compute +# --------------------------------- + +print("\nLoading model from JSON...") +model = JsonIO.load_model_from_json(json_file) + +print("\nModel Metadata:") +print(f"Name: {model.meta.name}") +print(f"Creation Date: {model.meta.creation_date}") +# TODO: This does not update here when running. In 03 you have a current date time setter +print(f"Last Modified: {model.meta.last_modification_date}") + +print("\nStructural Groups:") +print(model.structural_frame) + +# Compute the model +print("\nComputing the model...") +s = gp.compute_model( + gempy_model=model, + engine_config=gp.data.GemPyEngineConfig( + backend=AvailableBackends.numpy + ) +) + +# Save the computed model +# TODO: This is identical to the file saved before computation (no results stored) +JsonIO.save_model_to_json(model, computed_json_file) +print(f"\nSaved computed model to: {os.path.abspath(computed_json_file)}") + +#%% + +# Plot the results +print("\nGenerating plots...") + + +# 2D plots +gpv.plot_2d(model, direction='y', show_results=False) +gpv.plot_2d(model, direction='y', show_data=True, show_boundaries=True) +gpv.plot_2d(model, direction='x', show_data=True) + +# Plot the blocks accounting for fault blocks +gpv.plot_2d( + model=model, + override_regular_grid=model.solutions.raw_arrays.litho_faults_block, + show_data=True, kwargs_lithology={'cmap': 'Set1', 'norm': None} +) + +# 3D plot +gpv.plot_3d(model) + +print("\nDone! The model has been:") +print("1. Created and saved to:", json_file) +print("2. Loaded from JSON") +print("3. Computed") +print("4. Saved with computation results to:", computed_json_file) +print("5. Visualized with various plots") \ No newline at end of file diff --git a/examples/tutorials/z_other_tutorials/json_io/05_minimal_json.py b/examples/tutorials/z_other_tutorials/json_io/05_minimal_json.py new file mode 100644 index 000000000..507c46a3a --- /dev/null +++ b/examples/tutorials/z_other_tutorials/json_io/05_minimal_json.py @@ -0,0 +1,87 @@ +""" +Tutorial: Minimal JSON I/O with default metadata, interpolation options, and grid settings +This tutorial demonstrates how to use the JSON I/O functionality with minimal input, +relying on default metadata values, interpolation options, and grid settings. +""" + +# %% +import numpy as np +from datetime import datetime + +import gempy as gp +import gempy_viewer as gpv +import json +import pyvista as pv + +from gempy.modules.json_io.json_operations import JsonIO # Updated import path + + +# %% +# Create the corresponding minimal JSON model +model_data = { + "surface_points": [ + { + "x": 50.0, + "y": 0.0, + "z": -20.0, + "id": 0, + } + ], + "orientations": [ + { + "x": 50.0, + "y": 0.0, + "z": -20.0, + "G_x": 1.0, + "G_y": 0.0, + "G_z": 1.0, + "id": 0, + "polarity": 1 + } + ], + "grid_settings": { + "regular_grid_resolution": [100, 2, 100], + "regular_grid_extent": [0, 150, -10, 10, -100, 0] + } +} + +# %% +# Save the minimal model to JSON +with open("minimal_model.json", "w") as f: + json.dump(model_data, f, indent=4) + +# Load the model from JSON +geo_model = JsonIO.load_model_from_json("minimal_model.json") + +# Compute the geological model +gp.compute_model(geo_model) + +p2d = gpv.plot_2d(geo_model) +# %% + +# Print the model metadata (should use default values) +print("\nModel Metadata:") +print(f"Name: {geo_model.meta.name}") +print(f"Creation Date: {geo_model.meta.creation_date}") +print(f"Last Modification Date: {geo_model.meta.last_modification_date}") +print(f"Owner: {geo_model.meta.owner}") + +# Print the interpolation options (should use default values) +print("\nInterpolation Options:") +print(f"Range: {geo_model.interpolation_options.kernel_options.range}") +print(f"Mesh Extraction: {geo_model.interpolation_options.mesh_extraction}") + +# Print the grid settings (should use default values) +print("\nGrid Settings:") +print(f"Resolution: {geo_model.grid._dense_grid.resolution}") +print(f"Extent: {geo_model.grid._dense_grid.extent}") + +# Print the structural groups +print("\nStructural Groups:") +for group in geo_model.structural_frame.structural_groups: + print(group) + +# Save the loaded model to verify the metadata, interpolation options, and grid settings are preserved +JsonIO.save_model_to_json(geo_model, "minimal_model_loaded.json") + +# %% diff --git a/examples/tutorials/z_other_tutorials/json_io/05c_minimal_comparison.py b/examples/tutorials/z_other_tutorials/json_io/05c_minimal_comparison.py new file mode 100644 index 000000000..6b93b99a8 --- /dev/null +++ b/examples/tutorials/z_other_tutorials/json_io/05c_minimal_comparison.py @@ -0,0 +1,155 @@ +""" +Tutorial: Compare minimal GemPy model with JSON representation +This tutorial demonstrates how to create a minimal GemPy model and compare it with its JSON representation. +""" + +# %% +import gempy as gp +import gempy_viewer as gpv +import json +import numpy as np +from gempy.modules.json_io.json_operations import JsonIO + +# %% +# Set up the basic model +geo_model: gp.data.GeoModel = gp.create_geomodel( + project_name='Model1', + extent=[0, 150, -10, 10, -100, 0], + resolution=[100, 2, 100], + structural_frame=gp.data.StructuralFrame.initialize_default_structure() +) + +# Add a surface point +gp.add_surface_points( + geo_model=geo_model, + x=[50], + y=[0], + z=[-20], + elements_names=['surface1'] +) + +# Add an orientation +gp.add_orientations( + geo_model=geo_model, + x=[50], + y=[0], + z=[-68], + elements_names=['surface1'], + pole_vector=[[1, 0, 1]] +) + +# Set interpolation options +geo_model.interpolation_options.kernel_options.range = 10. +geo_model.interpolation_options.kernel_options.c_o = 5. +geo_model.interpolation_options.mesh_extraction = True +geo_model.interpolation_options.number_octree_levels = 2 + +geo_model.update_transform(gp.data.GlobalAnisotropy.NONE) +gp.compute_model(geo_model, engine_config=gp.data.GemPyEngineConfig()) + +# %% +# Plot the model +p2d = gpv.plot_2d(geo_model) + +# %% +# Create the corresponding minimal JSON model +json_model_data = { + "surface_points": [ + { + "x": 50.0, + "y": 0.0, + "z": -20.0, + "id": 0 + } + ], + "orientations": [ + { + "x": 50.0, + "y": 0.0, + "z": -68.0, + "G_x": 1.0, + "G_y": 0.0, + "G_z": 1.0, + "id": 0 + } + ], + "grid_settings": { + "regular_grid_resolution": [100, 2, 100], + "regular_grid_extent": [0, 150, -10, 10, -100, 0] + }, + "interpolation_options": { + "kernel_options": { + "range": 10.0, + "c_o": 5.0 + }, + "mesh_extraction": True, + "number_octree_levels": 2 + }, + "series": [ + { + "name": "default_formations", + "surfaces": ["surface1"], + "structural_relation": "ERODE" + } + ] +} + +# Save the JSON model +with open("minimal_model.json", "w") as f: + json.dump(json_model_data, f, indent=4) + +# Load the model from JSON +json_geo_model = JsonIO.load_model_from_json("minimal_model.json") + +# Compute the JSON model +gp.compute_model(json_geo_model, engine_config=gp.data.GemPyEngineConfig()) + +# %% +# Compare the models +print("\nComparing GemPy and JSON models:") +print("\n1. Surface Points:") +print("GemPy model:") +print(geo_model.surface_points_copy) +print("\nJSON model:") +print(json_geo_model.surface_points_copy) + +print("\n2. Orientations:") +print("GemPy model:") +print(geo_model.orientations_copy) +print("\nJSON model:") +print(json_geo_model.orientations_copy) + +print("\n3. Grid Settings:") +print("GemPy model:") +print(f"Resolution: {geo_model.grid._dense_grid.resolution}") +print(f"Extent: {geo_model.grid._dense_grid.extent}") +print("\nJSON model:") +print(f"Resolution: {json_geo_model.grid._dense_grid.resolution}") +print(f"Extent: {json_geo_model.grid._dense_grid.extent}") + +print("\n4. Interpolation Options:") +print("GemPy model:") +print(f"Range: {geo_model.interpolation_options.kernel_options.range}") +print(f"C_o: {geo_model.interpolation_options.kernel_options.c_o}") +print(f"Mesh Extraction: {geo_model.interpolation_options.mesh_extraction}") +print(f"Number Octree Levels: {geo_model.interpolation_options.number_octree_levels}") +print("\nJSON model:") +print(f"Range: {json_geo_model.interpolation_options.kernel_options.range}") +print(f"C_o: {json_geo_model.interpolation_options.kernel_options.c_o}") +print(f"Mesh Extraction: {json_geo_model.interpolation_options.mesh_extraction}") +print(f"Number Octree Levels: {json_geo_model.interpolation_options.number_octree_levels}") + +print("\n5. Structural Groups:") +print("GemPy model:") +for group in geo_model.structural_frame.structural_groups: + print(group) +print("\nJSON model:") +for group in json_geo_model.structural_frame.structural_groups: + print(group) + +# %% +# Plot both models side by side +p2d_gempy = gpv.plot_2d(geo_model, title="GemPy Model") +p2d_json = gpv.plot_2d(json_geo_model, title="JSON Model") + +# %% \ No newline at end of file diff --git a/examples/tutorials/z_other_tutorials/json_io/horizontal_stratigraphic.json b/examples/tutorials/z_other_tutorials/json_io/horizontal_stratigraphic.json new file mode 100644 index 000000000..8b139d49c --- /dev/null +++ b/examples/tutorials/z_other_tutorials/json_io/horizontal_stratigraphic.json @@ -0,0 +1,144 @@ +{ + "metadata": { + "name": "horizontal_stratigraphic", + "creation_date": "2024-03-19", + "last_modification_date": "2024-03-19", + "owner": "tutorial" + }, + "surface_points": [ + { + "x": 100.0, + "y": 200.0, + "z": 600.0, + "id": 1, + "nugget": 2e-05 + }, + { + "x": 500.0, + "y": 200.0, + "z": 600.0, + "id": 1, + "nugget": 2e-05 + }, + { + "x": 900.0, + "y": 200.0, + "z": 600.0, + "id": 1, + "nugget": 2e-05 + }, + { + "x": 100.0, + "y": 800.0, + "z": 600.0, + "id": 1, + "nugget": 2e-05 + }, + { + "x": 500.0, + "y": 800.0, + "z": 600.0, + "id": 1, + "nugget": 2e-05 + }, + { + "x": 900.0, + "y": 800.0, + "z": 600.0, + "id": 1, + "nugget": 2e-05 + }, + { + "x": 100.0, + "y": 200.0, + "z": 400.0, + "id": 0, + "nugget": 2e-05 + }, + { + "x": 500.0, + "y": 200.0, + "z": 400.0, + "id": 0, + "nugget": 2e-05 + }, + { + "x": 900.0, + "y": 200.0, + "z": 400.0, + "id": 0, + "nugget": 2e-05 + }, + { + "x": 100.0, + "y": 800.0, + "z": 400.0, + "id": 0, + "nugget": 2e-05 + }, + { + "x": 500.0, + "y": 800.0, + "z": 400.0, + "id": 0, + "nugget": 2e-05 + }, + { + "x": 900.0, + "y": 800.0, + "z": 400.0, + "id": 0, + "nugget": 2e-05 + } + ], + "orientations": [ + { + "x": 500.0, + "y": 500.0, + "z": 600.0, + "G_x": 0.0, + "G_y": 0.0, + "G_z": 1.0, + "id": 1, + "nugget": 0.01, + "polarity": 1 + }, + { + "x": 500.0, + "y": 500.0, + "z": 400.0, + "G_x": 0.0, + "G_y": 0.0, + "G_z": 1.0, + "id": 0, + "nugget": 0.01, + "polarity": 1 + } + ], + "series": [ + { + "name": "Strat_Series", + "surfaces": [ + "rock2", + "rock1" + ] + } + ], + "grid_settings": { + "regular_grid_resolution": [ + 10, + 10, + 10 + ], + "regular_grid_extent": [ + 0, + 1000, + 0, + 1000, + 0, + 1000 + ], + "octree_levels": null + }, + "interpolation_options": {} +} \ No newline at end of file diff --git a/examples/tutorials/z_other_tutorials/json_io/minimal_model.json b/examples/tutorials/z_other_tutorials/json_io/minimal_model.json new file mode 100644 index 000000000..8f389b28a --- /dev/null +++ b/examples/tutorials/z_other_tutorials/json_io/minimal_model.json @@ -0,0 +1,53 @@ +{ + "surface_points": [ + { + "x": 50.0, + "y": 0.0, + "z": -20.0, + "id": 0 + } + ], + "orientations": [ + { + "x": 50.0, + "y": 0.0, + "z": -68.0, + "G_x": 1.0, + "G_y": 0.0, + "G_z": 1.0, + "id": 0 + } + ], + "grid_settings": { + "regular_grid_resolution": [ + 100, + 2, + 100 + ], + "regular_grid_extent": [ + 0, + 150, + -10, + 10, + -100, + 0 + ] + }, + "interpolation_options": { + "kernel_options": { + "range": 10.0, + "c_o": 5.0 + }, + "mesh_extraction": true, + "number_octree_levels": 2 + }, + "series": [ + { + "name": "default_formations", + "surfaces": [ + "surface1" + ], + "structural_relation": "ERODE" + } + ] +} \ No newline at end of file diff --git a/examples/tutorials/z_other_tutorials/json_io/minimal_model_loaded.json b/examples/tutorials/z_other_tutorials/json_io/minimal_model_loaded.json new file mode 100644 index 000000000..24fa789c7 --- /dev/null +++ b/examples/tutorials/z_other_tutorials/json_io/minimal_model_loaded.json @@ -0,0 +1,77 @@ +{ + "metadata": { + "name": "GemPy Model", + "creation_date": "2025-03-28", + "last_modification_date": "2025-03-28", + "owner": "GemPy Modeller" + }, + "surface_points": [ + { + "x": 50.0, + "y": 0.0, + "z": -20.0, + "id": 34966260, + "nugget": 0.0 + } + ], + "orientations": [ + { + "x": 50.0, + "y": 0.0, + "z": -20.0, + "G_x": 1.0, + "G_y": 0.0, + "G_z": 1.0, + "id": 34966260, + "nugget": 0.01, + "polarity": 1 + } + ], + "series": [ + { + "name": "Strat_Series", + "surfaces": [ + "surface_0" + ], + "structural_relation": "ERODE", + "colors": [ + "#015482" + ] + } + ], + "grid_settings": { + "regular_grid_resolution": [ + 100, + 2, + 100 + ], + "regular_grid_extent": [ + 0.0, + 150.0, + -10.0, + 10.0, + -100.0, + 0.0 + ], + "octree_levels": null + }, + "interpolation_options": { + "kernel_options": { + "range": 1.7, + "c_o": 10.0 + }, + "mesh_extraction": true, + "number_octree_levels": 4 + }, + "fault_relations": [ + [ + 0 + ] + ], + "id_name_mapping": { + "name_to_id": { + "surface_0": 34966260, + "basement": 91927817 + } + } +} \ No newline at end of file diff --git a/examples/tutorials/z_other_tutorials/json_io/model_x_direction.png b/examples/tutorials/z_other_tutorials/json_io/model_x_direction.png new file mode 100644 index 000000000..d5c6228a6 Binary files /dev/null and b/examples/tutorials/z_other_tutorials/json_io/model_x_direction.png differ diff --git a/examples/tutorials/z_other_tutorials/json_io/model_y_direction.png b/examples/tutorials/z_other_tutorials/json_io/model_y_direction.png new file mode 100644 index 000000000..095b1024c Binary files /dev/null and b/examples/tutorials/z_other_tutorials/json_io/model_y_direction.png differ diff --git a/examples/tutorials/z_other_tutorials/json_io/multiple_series_faults.json b/examples/tutorials/z_other_tutorials/json_io/multiple_series_faults.json new file mode 100644 index 000000000..fd3dc9eb2 --- /dev/null +++ b/examples/tutorials/z_other_tutorials/json_io/multiple_series_faults.json @@ -0,0 +1,288 @@ +{ + "metadata": { + "name": "multiple_series_faults", + "creation_date": "2024-03-19", + "last_modification_date": "2024-03-19", + "owner": "tutorial" + }, + "surface_points": [ + { + "x": 0.0, + "y": 200.0, + "z": 600.0, + "id": 2, + "nugget": 2e-05 + }, + { + "x": 0.0, + "y": 500.0, + "z": 600.0, + "id": 2, + "nugget": 2e-05 + }, + { + "x": 0.0, + "y": 800.0, + "z": 600.0, + "id": 2, + "nugget": 2e-05 + }, + { + "x": 200.0, + "y": 200.0, + "z": 600.0, + "id": 2, + "nugget": 2e-05 + }, + { + "x": 200.0, + "y": 500.0, + "z": 600.0, + "id": 2, + "nugget": 2e-05 + }, + { + "x": 200.0, + "y": 800.0, + "z": 600.0, + "id": 2, + "nugget": 2e-05 + }, + { + "x": 800.0, + "y": 200.0, + "z": 200.0, + "id": 2, + "nugget": 2e-05 + }, + { + "x": 800.0, + "y": 500.0, + "z": 200.0, + "id": 2, + "nugget": 2e-05 + }, + { + "x": 800.0, + "y": 800.0, + "z": 200.0, + "id": 2, + "nugget": 2e-05 + }, + { + "x": 1000.0, + "y": 200.0, + "z": 200.0, + "id": 2, + "nugget": 2e-05 + }, + { + "x": 1000.0, + "y": 500.0, + "z": 200.0, + "id": 2, + "nugget": 2e-05 + }, + { + "x": 1000.0, + "y": 800.0, + "z": 200.0, + "id": 2, + "nugget": 2e-05 + }, + { + "x": 0.0, + "y": 200.0, + "z": 800.0, + "id": 1, + "nugget": 2e-05 + }, + { + "x": 0.0, + "y": 800.0, + "z": 800.0, + "id": 1, + "nugget": 2e-05 + }, + { + "x": 200.0, + "y": 200.0, + "z": 800.0, + "id": 1, + "nugget": 2e-05 + }, + { + "x": 200.0, + "y": 800.0, + "z": 800.0, + "id": 1, + "nugget": 2e-05 + }, + { + "x": 800.0, + "y": 200.0, + "z": 400.0, + "id": 1, + "nugget": 2e-05 + }, + { + "x": 800.0, + "y": 800.0, + "z": 400.0, + "id": 1, + "nugget": 2e-05 + }, + { + "x": 1000.0, + "y": 200.0, + "z": 400.0, + "id": 1, + "nugget": 2e-05 + }, + { + "x": 1000.0, + "y": 800.0, + "z": 400.0, + "id": 1, + "nugget": 2e-05 + }, + { + "x": 500.0, + "y": 500.0, + "z": 500.0, + "id": 0, + "nugget": 2e-05 + }, + { + "x": 450.0, + "y": 500.0, + "z": 600.0, + "id": 0, + "nugget": 2e-05 + }, + { + "x": 500.0, + "y": 200.0, + "z": 500.0, + "id": 0, + "nugget": 2e-05 + }, + { + "x": 450.0, + "y": 200.0, + "z": 600.0, + "id": 0, + "nugget": 2e-05 + }, + { + "x": 500.0, + "y": 800.0, + "z": 500.0, + "id": 0, + "nugget": 2e-05 + }, + { + "x": 450.0, + "y": 800.0, + "z": 600.0, + "id": 0, + "nugget": 2e-05 + } + ], + "orientations": [ + { + "x": 100.0, + "y": 500.0, + "z": 800.0, + "G_x": 0.0, + "G_y": 0.0, + "G_z": 1.0, + "id": 1, + "nugget": 0.01, + "polarity": 1 + }, + { + "x": 100.0, + "y": 500.0, + "z": 600.0, + "G_x": 0.0, + "G_y": 0.0, + "G_z": 1.0, + "id": 2, + "nugget": 0.01, + "polarity": 1 + }, + { + "x": 500.0, + "y": 500.0, + "z": 500.0, + "G_x": 0.8, + "G_y": 0.0, + "G_z": 0.6, + "id": 0, + "nugget": 0.01, + "polarity": 1 + }, + { + "x": 900.0, + "y": 500.0, + "z": 400.0, + "G_x": 0.0, + "G_y": 0.0, + "G_z": 1.0, + "id": 1, + "nugget": 0.01, + "polarity": 1 + }, + { + "x": 900.0, + "y": 500.0, + "z": 200.0, + "G_x": 0.0, + "G_y": 0.0, + "G_z": 1.0, + "id": 2, + "nugget": 0.01, + "polarity": 1 + } + ], + "series": [ + { + "name": "Fault_Series", + "surfaces": [ + "fault" + ], + "structural_relation": "FAULT", + "color": "#015482" + }, + { + "name": "Strat_Series", + "surfaces": [ + "rock2", + "rock1" + ], + "structural_relation": "ERODE", + "colors": [ + "#ffbe00", + "#9f0052" + ] + } + ], + "grid_settings": { + "regular_grid_resolution": [ + 90, + 30, + 30 + ], + "regular_grid_extent": [ + 0, + 1000, + 0, + 1000, + 0, + 1000 + ], + "octree_levels": null + }, + "interpolation_options": {} +} \ No newline at end of file diff --git a/examples/tutorials/z_other_tutorials/json_io/sample_surface_points.json b/examples/tutorials/z_other_tutorials/json_io/sample_surface_points.json new file mode 100644 index 000000000..a1f39c42d --- /dev/null +++ b/examples/tutorials/z_other_tutorials/json_io/sample_surface_points.json @@ -0,0 +1,65 @@ +{ + "metadata": { + "name": "sample_model", + "creation_date": "2024-03-19", + "last_modification_date": "2024-03-19", + "owner": "tutorial" + }, + "surface_points": [ + { + "x": 0.0, + "y": 0.0, + "z": 0.0, + "id": 0, + "nugget": 2e-05 + }, + { + "x": 1.0, + "y": 1.0, + "z": 1.0, + "id": 0, + "nugget": 2e-05 + }, + { + "x": 2.0, + "y": 2.0, + "z": 2.0, + "id": 1, + "nugget": 2e-05 + }, + { + "x": 3.0, + "y": 3.0, + "z": 3.0, + "id": 1, + "nugget": 2e-05 + }, + { + "x": 4.0, + "y": 4.0, + "z": 4.0, + "id": 2, + "nugget": 2e-05 + } + ], + "orientations": [], + "faults": [], + "series": [], + "grid_settings": { + "regular_grid_resolution": [ + 10, + 10, + 10 + ], + "regular_grid_extent": [ + 0, + 4, + 0, + 4, + 0, + 4 + ], + "octree_levels": null + }, + "interpolation_options": {} +} \ No newline at end of file diff --git a/gempy/API/map_stack_to_surfaces_API.py b/gempy/API/map_stack_to_surfaces_API.py index ba2ee8d43..c6feb31c9 100644 --- a/gempy/API/map_stack_to_surfaces_API.py +++ b/gempy/API/map_stack_to_surfaces_API.py @@ -7,7 +7,7 @@ def map_stack_to_surfaces(gempy_model: GeoModel, mapping_object: Union[dict[str, list[str]] | dict[str, tuple]], - set_series: bool = True, remove_unused_series=True) -> StructuralFrame: + set_series: bool = True, remove_unused_series=True, series_data: list = None) -> StructuralFrame: """ Map stack (series) to surfaces by reorganizing elements between groups in a GeoModel's structural frame. @@ -20,6 +20,7 @@ def map_stack_to_surfaces(gempy_model: GeoModel, mapping_object: Union[dict[str, mapping_object (Union[dict[str, list[str]] | dict[str, tuple]]): Dictionary mapping group names to element names. set_series (bool, optional): If True, creates new series for groups not present in the GeoModel. Defaults to True. remove_unused_series (bool, optional): If True, removes groups without any elements. Defaults to True. + series_data (list, optional): List of series data from JSON containing structural relations. Defaults to None. Returns: StructuralFrame: The updated StructuralFrame object. @@ -27,14 +28,21 @@ def map_stack_to_surfaces(gempy_model: GeoModel, mapping_object: Union[dict[str, structural_groups: list[StructuralGroup] = gempy_model.structural_frame.structural_groups for index, (group_name, elements) in enumerate(mapping_object.items()): - # region Create new series if needed group_already_exists = any(group.name == group_name for group in structural_groups) if set_series and not group_already_exists: + # Get structural relation from series_data if available + structural_relation = StackRelationType.ERODE # Default value + if series_data: + for series in series_data: + if series['name'] == group_name: + structural_relation = StackRelationType[series['structural_relation']] + break + new_group = StructuralGroup( name=group_name, elements=[], - structural_relation=StackRelationType.ERODE + structural_relation=structural_relation ) structural_groups.insert(index, new_group) # endregion diff --git a/gempy/modules/__init__.py b/gempy/modules/__init__.py index e69de29bb..e37cb4efb 100644 --- a/gempy/modules/__init__.py +++ b/gempy/modules/__init__.py @@ -0,0 +1,17 @@ +""" +GemPy modules package. +""" + +from . import json_io +from . import custom_implicit_functions +from . import data_manipulation +from . import grids +from . import advance_pile + +__all__ = [ + 'json_io', + 'custom_implicit_functions', + 'data_manipulation', + 'grids', + 'advance_pile' +] diff --git a/gempy/modules/json_io/__init__.py b/gempy/modules/json_io/__init__.py new file mode 100644 index 000000000..c0d3ecf87 --- /dev/null +++ b/gempy/modules/json_io/__init__.py @@ -0,0 +1,8 @@ +""" +JSON I/O module for GemPy. +This module provides functionality to load and save GemPy models to/from JSON files. +""" + +from .json_operations import JsonIO + +__all__ = ['JsonIO'] \ No newline at end of file diff --git a/gempy/modules/json_io/json_operations.py b/gempy/modules/json_io/json_operations.py new file mode 100644 index 000000000..aefc1ce61 --- /dev/null +++ b/gempy/modules/json_io/json_operations.py @@ -0,0 +1,536 @@ +""" +Module for JSON I/O operations in GemPy. +This module provides functionality to load and save GemPy models to/from JSON files. +""" + +import json +from typing import Dict, Any, Optional, List +import numpy as np +from datetime import datetime + +from .schema import SurfacePoint, Orientation, GemPyModelJson, IdNameMapping +from gempy_engine.core.data.stack_relation_type import StackRelationType + + +class JsonIO: + """Class for handling JSON I/O operations for GemPy models.""" + + @staticmethod + def _numpy_to_list(obj): + """Convert numpy arrays to lists for JSON serialization.""" + if isinstance(obj, np.ndarray): + return obj.tolist() + elif isinstance(obj, (np.integer, np.int32, np.int64)): + return int(obj) + elif isinstance(obj, (np.floating, np.float32, np.float64)): + return float(obj) + return obj + + @staticmethod + def _create_id_to_name(name_to_id: Dict[str, int]) -> Dict[int, str]: + """Create an id_to_name mapping from a name_to_id mapping.""" + return {id: name for name, id in name_to_id.items()} + + @staticmethod + def _calculate_default_range(grid_extent: List[float]) -> float: + """Calculate the default range based on the model extent (room diagonal).""" + # Extract min and max coordinates + x_min, x_max = grid_extent[0], grid_extent[1] + y_min, y_max = grid_extent[2], grid_extent[3] + z_min, z_max = grid_extent[4], grid_extent[5] + + # Calculate the room diagonal (Euclidean distance) + return np.sqrt((x_max - x_min)**2 + (y_max - y_min)**2 + (z_max - z_min)**2) + + @staticmethod + def _calculate_default_grid_settings(surface_points: List[SurfacePoint], orientations: List[Orientation]) -> Dict[str, Any]: + """Calculate default grid settings based on data points. + + Args: + surface_points: List of surface points + orientations: List of orientations + + Returns: + Dict containing grid settings with default values + """ + # Collect all x, y, z coordinates + all_x = [sp['x'] for sp in surface_points] + [ori['x'] for ori in orientations] + all_y = [sp['y'] for sp in surface_points] + [ori['y'] for ori in orientations] + all_z = [sp['z'] for sp in surface_points] + [ori['z'] for ori in orientations] + + # Calculate extents + x_min, x_max = min(all_x), max(all_x) + y_min, y_max = min(all_y), max(all_y) + z_min, z_max = min(all_z), max(all_z) + + # Calculate ranges + x_range = x_max - x_min + y_range = y_max - y_min + z_range = z_max - z_min + + # Add 10% padding to each dimension + x_padding = x_range * 0.1 + y_padding = y_range * 0.1 + z_padding = z_range * 0.1 + + return { + "regular_grid_resolution": [20, 20, 20], + "regular_grid_extent": [ + x_min - x_padding, x_max + x_padding, + y_min - y_padding, y_max + y_padding, + z_min - z_padding, z_max + z_padding + ] + } + + @staticmethod + def load_model_from_json(file_path: str): + """ + Load a GemPy model from a JSON file. + + Args: + file_path (str): Path to the JSON file + + Returns: + GeoModel: A new GemPy model instance + """ + # Import here to avoid circular imports + from gempy.core.data.geo_model import GeoModel, GeoModelMeta + from gempy.core.data.grid import Grid + from gempy.core.data.structural_frame import StructuralFrame + from gempy_engine.core.data import InterpolationOptions + from gempy.API.map_stack_to_surfaces_API import map_stack_to_surfaces + + with open(file_path, 'r') as f: + data = json.load(f) + + # Validate the JSON data against our schema + if not JsonIO._validate_json_schema(data): + raise ValueError("Invalid JSON schema") + + # Get unique surface IDs from surface points and orientations + surface_ids = set(sp['id'] for sp in data['surface_points']) + surface_ids.update(ori['id'] for ori in data['orientations']) + + # Create default series if not provided + if 'series' not in data: + data['series'] = [{ + 'name': 'Strat_Series', + 'surfaces': [f'surface_{id}' for id in sorted(surface_ids)], + 'structural_relation': 'ERODE' + }] + + # Get surface names from series data + surface_names = [] + for series in data['series']: + surface_names.extend(series['surfaces']) + + # Create ID to name mapping if not provided + if 'id_name_mapping' in data: + id_to_name = JsonIO._create_id_to_name(data['id_name_mapping']['name_to_id']) + else: + # Create mapping from series data + id_to_name = {i: name for i, name in enumerate(surface_names)} + + # Create surface points table + surface_points_table = JsonIO._load_surface_points(data['surface_points'], id_to_name) + + # Create orientations table + orientations_table = JsonIO._load_orientations(data['orientations'], id_to_name) + + # Create structural frame + structural_frame = StructuralFrame.from_data_tables(surface_points_table, orientations_table) + + # Get grid settings with defaults if not provided + grid_settings = data.get('grid_settings', JsonIO._calculate_default_grid_settings(data['surface_points'], data['orientations'])) + + # Create grid + grid = Grid( + resolution=grid_settings['regular_grid_resolution'], + extent=grid_settings['regular_grid_extent'] + ) + + # Calculate default range based on model extent + # default_range = JsonIO._calculate_default_range(grid_settings['regular_grid_extent']) + # set as fixed value + default_range = 1.7 # Match standard GemPy default + + # Create interpolation options with defaults if not provided + interpolation_options = InterpolationOptions( + range=data.get('interpolation_options', {}).get('kernel_options', {}).get('range', default_range), + c_o=data.get('interpolation_options', {}).get('kernel_options', {}).get('c_o', 10), + mesh_extraction=data.get('interpolation_options', {}).get('mesh_extraction', True), + number_octree_levels=data.get('interpolation_options', {}).get('number_octree_levels', 1) + ) + + # Create GeoModel with default metadata if not provided + current_date = datetime.now().strftime("%Y-%m-%d") + model_name = data.get('metadata', {}).get('name', "GemPy Model") + + model = GeoModel( + name=model_name, + structural_frame=structural_frame, + grid=grid, + interpolation_options=interpolation_options + ) + + # Set the metadata with proper dates and defaults + metadata = data.get('metadata', {}) + model_meta = GeoModelMeta( + name=metadata.get('name', model.meta.name), # Use model's name if available + creation_date=metadata.get('creation_date', current_date), # Set current date as default + last_modification_date=metadata.get('last_modification_date', current_date), # Set current date as default + owner=metadata.get('owner', "GemPy Modeller") # Set default owner + ) + model.meta = model_meta + + # Map series to surfaces with structural relations + mapping_object = {series['name']: series['surfaces'] for series in data['series']} + # Ensure each series has structural_relation set to ERODE by default + series_data = [] + for series in data['series']: + series_copy = series.copy() + if 'structural_relation' not in series_copy: + series_copy['structural_relation'] = 'ERODE' + series_data.append(series_copy) + map_stack_to_surfaces(model, mapping_object, series_data=series_data) + + # Set fault relations after structural groups are set up + if 'fault_relations' in data and data['fault_relations'] is not None: + fault_relations = np.array(data['fault_relations']) + if fault_relations.shape == (len(model.structural_frame.structural_groups), len(model.structural_frame.structural_groups)): + model.structural_frame.fault_relations = fault_relations + + # Set colors for each element + for series in data['series']: + if 'colors' in series: + for surface_name, color in zip(series['surfaces'], series['colors']): + element = model.structural_frame.get_element_by_name(surface_name) + if element is not None: + element.color = color + + return model + + @staticmethod + def _load_surface_points(surface_points_data: List[SurfacePoint], id_to_name: Optional[Dict[int, str]] = None): + """ + Load surface points from JSON data. + + Args: + surface_points_data (List[SurfacePoint]): List of surface point dictionaries + id_to_name (Optional[Dict[int, str]]): Optional mapping from surface IDs to names + + Returns: + SurfacePointsTable: A new SurfacePointsTable instance + + Raises: + ValueError: If the data is invalid or missing required fields + """ + # Import here to avoid circular imports + from gempy.core.data.surface_points import SurfacePointsTable + + # Validate data structure + required_fields = {'x', 'y', 'z', 'nugget', 'id'} + for i, sp in enumerate(surface_points_data): + missing_fields = required_fields - set(sp.keys()) + if missing_fields: + raise ValueError(f"Missing required fields in surface point {i}: {missing_fields}") + + # Validate data types + if not all(isinstance(sp[field], (int, float)) for field in ['x', 'y', 'z', 'nugget']): + raise ValueError(f"Invalid data type in surface point {i}. All coordinates and nugget must be numeric.") + if not isinstance(sp['id'], int): + raise ValueError(f"Invalid data type in surface point {i}. ID must be an integer.") + + # Extract coordinates and other data + x = np.array([sp['x'] for sp in surface_points_data]) + y = np.array([sp['y'] for sp in surface_points_data]) + z = np.array([sp['z'] for sp in surface_points_data]) + nugget = np.array([sp['nugget'] for sp in surface_points_data]) + + # Handle names based on whether id_to_name mapping is provided + if id_to_name is not None: + names = [id_to_name.get(sp['id'], f"surface_{sp['id']}") for sp in surface_points_data] + else: + # If no mapping provided, use surface IDs as names + names = [f"surface_{sp['id']}" for sp in surface_points_data] + + # Create SurfacePointsTable + return SurfacePointsTable.from_arrays( + x=x, + y=y, + z=z, + names=names, + nugget=nugget + ) + + @staticmethod + def _load_orientations(orientations_data: List[Orientation], id_to_name: Optional[Dict[int, str]] = None): + """ + Load orientations from JSON data. + + Args: + orientations_data (List[Orientation]): List of orientation dictionaries + id_to_name (Optional[Dict[int, str]]): Optional mapping from surface IDs to names + + Returns: + OrientationsTable: A new OrientationsTable instance + + Raises: + ValueError: If the data is invalid or missing required fields + """ + # Import here to avoid circular imports + from gempy.core.data.orientations import OrientationsTable + + # Validate data structure + required_fields = {'x', 'y', 'z', 'G_x', 'G_y', 'G_z', 'nugget', 'polarity', 'id'} + for i, ori in enumerate(orientations_data): + missing_fields = required_fields - set(ori.keys()) + if missing_fields: + raise ValueError(f"Missing required fields in orientation {i}: {missing_fields}") + + # Validate data types + if not all(isinstance(ori[field], (int, float)) for field in ['x', 'y', 'z', 'G_x', 'G_y', 'G_z', 'nugget']): + raise ValueError(f"Invalid data type in orientation {i}. All coordinates, gradients, and nugget must be numeric.") + if not isinstance(ori.get('polarity', 1), int) or ori.get('polarity', 1) not in {-1, 1}: + raise ValueError(f"Invalid polarity in orientation {i}. Must be 1 (normal) or -1 (reverse).") + if not isinstance(ori['id'], int): + raise ValueError(f"Invalid data type in orientation {i}. ID must be an integer.") + + # Extract coordinates and other data + x = np.array([ori['x'] for ori in orientations_data]) + y = np.array([ori['y'] for ori in orientations_data]) + z = np.array([ori['z'] for ori in orientations_data]) + G_x = np.array([ori['G_x'] for ori in orientations_data]) + G_y = np.array([ori['G_y'] for ori in orientations_data]) + G_z = np.array([ori['G_z'] for ori in orientations_data]) + nugget = np.array([ori['nugget'] for ori in orientations_data]) + + # Handle names based on whether id_to_name mapping is provided + if id_to_name is not None: + names = [id_to_name.get(ori['id'], f"surface_{ori['id']}") for ori in orientations_data] + else: + # If no mapping provided, use surface IDs as names + names = [f"surface_{ori['id']}" for ori in orientations_data] + + # Apply polarity to gradients + for i, ori in enumerate(orientations_data): + if ori.get('polarity', 1) == -1: + G_x[i] *= -1 + G_y[i] *= -1 + G_z[i] *= -1 + + # Create OrientationsTable + return OrientationsTable.from_arrays( + x=x, + y=y, + z=z, + G_x=G_x, + G_y=G_y, + G_z=G_z, + names=names, + nugget=nugget + ) + + @staticmethod + def save_model_to_json(model, file_path: str) -> None: + """ + Save a GemPy model to a JSON file. + + Args: + model: The GemPy model to save + file_path (str): Path where to save the JSON file + """ + # Get current date for default values + current_date = datetime.now().strftime("%Y-%m-%d") + + # Create JSON structure with metadata handling + json_data = { + "metadata": { + "name": model.meta.name if model.meta.name is not None else "GemPy Model", + "creation_date": model.meta.creation_date if model.meta.creation_date is not None else current_date, + "last_modification_date": model.meta.last_modification_date if model.meta.last_modification_date is not None else current_date, + "owner": model.meta.owner if model.meta.owner is not None else "GemPy Modeller" + }, + "surface_points": [], + "orientations": [], + "series": [], + "grid_settings": { + "regular_grid_resolution": [int(x) for x in model.grid._dense_grid.resolution], + "regular_grid_extent": [float(x) for x in model.grid._dense_grid.extent], + "octree_levels": None # TODO: Add octree levels if needed + }, + "interpolation_options": { + "kernel_options": { + "range": float(model.interpolation_options.kernel_options.range), + "c_o": float(model.interpolation_options.kernel_options.c_o) + }, + "mesh_extraction": bool(model.interpolation_options.mesh_extraction), + "number_octree_levels": int(model.interpolation_options.number_octree_levels) + }, + "fault_relations": [[int(x) for x in row] for row in model.structural_frame.fault_relations] if hasattr(model.structural_frame, 'fault_relations') else None, + "id_name_mapping": { + "name_to_id": {k: int(v) for k, v in model.structural_frame.element_name_id_map.items()} + } + } + + # Get series and surface information + for group in model.structural_frame.structural_groups: + series_entry = { + "name": str(group.name), + "surfaces": [str(element.name) for element in group.elements], + "structural_relation": str(group.structural_relation.name), + "colors": [str(element.color) for element in group.elements] + } + json_data["series"].append(series_entry) + + # Get surface points + surface_points_table = model.surface_points_copy + xyz_values = surface_points_table.xyz + ids = surface_points_table.ids + nugget_values = surface_points_table.nugget + + for i in range(len(xyz_values)): + point = { + "x": float(xyz_values[i, 0]), + "y": float(xyz_values[i, 1]), + "z": float(xyz_values[i, 2]), + "id": int(ids[i]), + "nugget": float(nugget_values[i]) + } + json_data["surface_points"].append(point) + + # Get orientations + orientations_table = model.orientations_copy + ori_xyz_values = orientations_table.xyz + ori_grads_values = orientations_table.grads + ori_ids = orientations_table.ids + ori_nugget_values = orientations_table.nugget + + for i in range(len(ori_xyz_values)): + orientation = { + "x": float(ori_xyz_values[i, 0]), + "y": float(ori_xyz_values[i, 1]), + "z": float(ori_xyz_values[i, 2]), + "G_x": float(ori_grads_values[i, 0]), + "G_y": float(ori_grads_values[i, 1]), + "G_z": float(ori_grads_values[i, 2]), + "id": int(ori_ids[i]), + "nugget": float(ori_nugget_values[i]), + "polarity": 1 # Default value, update if available + } + json_data["orientations"].append(orientation) + + # Save to file + with open(file_path, 'w') as f: + json.dump(json_data, f, indent=4) + + @staticmethod + def _validate_json_schema(data: Dict) -> None: + """Validate the JSON schema and set default values.""" + # Required top-level keys + required_keys = ['surface_points', 'orientations', 'grid_settings'] + for key in required_keys: + if key not in data: + raise ValueError(f"Missing required key: {key}") + + # Validate surface points + for point in data['surface_points']: + required_point_keys = ['x', 'y', 'z'] + for key in required_point_keys: + if key not in point: + raise ValueError(f"Missing required key in surface point: {key}") + # Set default nugget if not provided + if 'nugget' not in point: + point['nugget'] = 0.0 # Default nugget for surface points + # Set default id if not provided + if 'id' not in point: + point['id'] = 0 # Default id + + # Validate orientations + for orientation in data['orientations']: + required_orientation_keys = ['x', 'y', 'z', 'G_x', 'G_y', 'G_z'] + for key in required_orientation_keys: + if key not in orientation: + raise ValueError(f"Missing required key in orientation: {key}") + # Set default nugget if not provided + if 'nugget' not in orientation: + orientation['nugget'] = 0.01 # Default nugget for orientations + # Set default polarity if not provided + if 'polarity' not in orientation: + orientation['polarity'] = 1 + # Set default id if not provided + if 'id' not in orientation: + orientation['id'] = 0 # Default id + + # Validate grid settings + grid_settings = data['grid_settings'] + required_grid_keys = ['regular_grid_resolution', 'regular_grid_extent'] + for key in required_grid_keys: + if key not in grid_settings: + raise ValueError(f"Missing required key in grid_settings: {key}") + + # Validate series if provided + if 'series' in data: + for series in data['series']: + required_series_keys = ['name', 'surfaces'] + for key in required_series_keys: + if key not in series: + raise ValueError(f"Missing required key in series: {key}") + # Set default structural relation if not provided + if 'structural_relation' not in series: + series['structural_relation'] = "ERODE" + # Validate colors if provided + if 'colors' in series: + if not isinstance(series['colors'], list): + raise ValueError("Colors must be a list") + if not all(isinstance(color, str) and color.startswith('#') for color in series['colors']): + raise ValueError("Colors must be hex color codes starting with #") + + # Validate interpolation options if present + if 'interpolation_options' in data: + if not isinstance(data['interpolation_options'], dict): + raise ValueError("Interpolation options must be a dictionary") + if 'kernel_options' in data['interpolation_options']: + kernel_options = data['interpolation_options']['kernel_options'] + if not isinstance(kernel_options, dict): + raise ValueError("Kernel options must be a dictionary") + if 'range' in kernel_options and not isinstance(kernel_options['range'], (int, float)): + raise ValueError("Kernel range must be a number") + if 'c_o' in kernel_options and not isinstance(kernel_options['c_o'], (int, float)): + raise ValueError("Kernel c_o must be a number") + if 'mesh_extraction' in data['interpolation_options'] and not isinstance(data['interpolation_options']['mesh_extraction'], bool): + raise ValueError("Mesh extraction must be a boolean") + if 'number_octree_levels' in data['interpolation_options'] and not isinstance(data['interpolation_options']['number_octree_levels'], int): + raise ValueError("Number of octree levels must be an integer") + + # Validate and set default metadata + if 'metadata' not in data: + data['metadata'] = {} + + metadata = data['metadata'] + if not isinstance(metadata, dict): + raise ValueError("Metadata must be a dictionary") + + # Set default values for metadata + current_date = datetime.now().strftime("%Y-%m-%d") + if 'name' not in metadata or metadata['name'] is None: + metadata['name'] = "GemPy Model" + elif not isinstance(metadata['name'], str): + raise ValueError("Metadata name must be a string") + + if 'creation_date' not in metadata or metadata['creation_date'] is None: + metadata['creation_date'] = current_date + elif not isinstance(metadata['creation_date'], str): + raise ValueError("Metadata creation_date must be a string") + + if 'last_modification_date' not in metadata or metadata['last_modification_date'] is None: + metadata['last_modification_date'] = current_date + elif not isinstance(metadata['last_modification_date'], str): + raise ValueError("Metadata last_modification_date must be a string") + + if 'owner' not in metadata or metadata['owner'] is None: + metadata['owner'] = "GemPy Modeller" + elif not isinstance(metadata['owner'], str): + raise ValueError("Metadata owner must be a string") + + return True \ No newline at end of file diff --git a/gempy/modules/json_io/schema.py b/gempy/modules/json_io/schema.py new file mode 100644 index 000000000..eae4da7f7 --- /dev/null +++ b/gempy/modules/json_io/schema.py @@ -0,0 +1,70 @@ +""" +Schema definitions for JSON I/O operations in GemPy. +This module defines the expected structure of JSON files for loading and saving GemPy models. +""" + +from typing import TypedDict, List, Dict, Any, Optional, Union, Sequence, NotRequired + +class SurfacePoint(TypedDict): + x: float + y: float + z: float + id: int + nugget: float + +class Orientation(TypedDict): + x: float + y: float + z: float + G_x: float # X component of the gradient + G_y: float # Y component of the gradient + G_z: float # Z component of the gradient + id: int + nugget: float + polarity: int # 1 for normal, -1 for reverse + +class Surface(TypedDict, total=False): + name: str # Required + id: NotRequired[int] # Optional, will be auto-generated if not provided + color: NotRequired[Optional[str]] # Optional hex color code + vertices: NotRequired[Optional[List[List[float]]]] # Optional list of [x, y, z] coordinates + +class Fault(TypedDict, total=False): + name: str # Required + id: NotRequired[int] # Optional, will be auto-generated + is_active: NotRequired[bool] # Optional, defaults to True + surface: Surface + +class Series(TypedDict, total=False): + name: str # Required + surfaces: Union[List[str], List[Surface]] # Required, can be list of names or Surface objects + id: NotRequired[int] # Optional, will be auto-generated + is_active: NotRequired[bool] # Optional, defaults to True + is_fault: NotRequired[bool] # Optional, defaults to False + order_series: NotRequired[int] # Optional, will be auto-generated + faults: NotRequired[List[Fault]] # Optional, defaults to empty list + structural_relation: NotRequired[str] # Optional, defaults to "ONLAP" + colors: NotRequired[Optional[List[str]]] # Optional + +class GridSettings(TypedDict, total=False): + regular_grid_resolution: NotRequired[List[int]] # Optional, defaults to [10, 10, 10] + regular_grid_extent: NotRequired[List[float]] # Optional, auto-calculated from data + octree_levels: NotRequired[Optional[int]] # Optional + +class ModelMetadata(TypedDict, total=False): + name: NotRequired[str] # Optional, defaults to "GemPy Model" + creation_date: NotRequired[str] # Optional, defaults to current date + last_modification_date: NotRequired[str] # Optional, defaults to current date + owner: NotRequired[Optional[str]] # Optional, defaults to "GemPy Team" + +class IdNameMapping(TypedDict, total=False): + name_to_id: Dict[str, int] # Required if id_name_mapping is provided + +class GemPyModelJson(TypedDict): + surface_points: List[SurfacePoint] # Required + orientations: List[Orientation] # Required + series: List[Series] # Required but with minimal required fields + metadata: NotRequired[ModelMetadata] # Optional + grid_settings: NotRequired[Optional[GridSettings]] # Optional + interpolation_options: NotRequired[Optional[Dict[str, Any]]] # Optional + id_name_mapping: NotRequired[IdNameMapping] # Optional \ No newline at end of file diff --git a/test/test_modules/test_json_io.py b/test/test_modules/test_json_io.py new file mode 100644 index 000000000..0fc49ce7a --- /dev/null +++ b/test/test_modules/test_json_io.py @@ -0,0 +1,816 @@ +""" +Tests for JSON I/O operations in GemPy. +""" + +import json +import numpy as np +import pytest +import gempy as gp +from gempy.modules.json_io import JsonIO +from gempy_engine.core.data.stack_relation_type import StackRelationType +from gempy.core.data.importer_helper import ImporterHelper +import pandas as pd +from gempy.core.data.structural_group import FaultsRelationSpecialCase + + +@pytest.fixture +def sample_surface_points(): + """Create sample surface points data for testing.""" + x = np.array([0, 1, 2, 3, 4]) + y = np.array([0, 1, 2, 3, 4]) + z = np.array([0, 1, 2, 3, 4]) + nugget = np.array([0.00002, 0.00002, 0.00002, 0.00002, 0.00002]) + + # Create a SurfacePointsTable + surface_points = gp.data.SurfacePointsTable.from_arrays( + x=x, + y=y, + z=z, + names=["surface_0", "surface_0", "surface_1", "surface_1", "surface_2"], + nugget=nugget + ) + + return surface_points, x, y, z, nugget + + +@pytest.fixture +def sample_orientations(): + """Create sample orientation data for testing.""" + x = np.array([0.5, 1.5, 2.5, 3.5]) + y = np.array([0.5, 1.5, 2.5, 3.5]) + z = np.array([0.5, 1.5, 2.5, 3.5]) + G_x = np.array([0, 0, 0, 0]) + G_y = np.array([0, 0, 0, 0]) + G_z = np.array([1, 1, -1, 1]) # One reversed orientation + nugget = np.array([0.01, 0.01, 0.01, 0.01]) + + # Create an OrientationsTable + orientations = gp.data.OrientationsTable.from_arrays( + x=x, + y=y, + z=z, + G_x=G_x, + G_y=G_y, + G_z=G_z, + names=["surface_0", "surface_1", "surface_1", "surface_2"], + nugget=nugget + ) + + return orientations, x, y, z, G_x, G_y, G_z, nugget + + +@pytest.fixture +def sample_json_data(sample_surface_points, sample_orientations): + """Create sample JSON data for testing.""" + _, x_sp, y_sp, z_sp, nugget_sp = sample_surface_points + _, x_ori, y_ori, z_ori, G_x, G_y, G_z, nugget_ori = sample_orientations + + return { + "metadata": { + "name": "sample_model", + "creation_date": "2024-03-19", + "last_modification_date": "2024-03-19", + "owner": "tutorial" + }, + "surface_points": [ + { + "x": float(x_sp[i]), + "y": float(y_sp[i]), + "z": float(z_sp[i]), + "id": 0, # Default ID, not used in tests + "nugget": float(nugget_sp[i]) + } + for i in range(len(x_sp)) + ], + "orientations": [ + { + "x": float(x_ori[i]), + "y": float(y_ori[i]), + "z": float(z_ori[i]), + "G_x": float(G_x[i]), + "G_y": float(G_y[i]), + "G_z": float(G_z[i]), + "id": 0, # Default ID, not used in tests + "nugget": float(nugget_ori[i]), + "polarity": 1 # Always set to 1 since we're testing the raw G_z values + } + for i in range(len(x_ori)) + ], + "series": [ + { + "name": "series1", + "surfaces": ["surface_0", "surface_1"], + "structural_relation": "erode", + "colors": ["#FF0000", "#00FF00"] + }, + { + "name": "series2", + "surfaces": ["surface_2"], + "structural_relation": "erode", + "colors": ["#0000FF"] + } + ], + "grid_settings": { + "regular_grid_resolution": [10, 10, 10], + "regular_grid_extent": [0, 4, 0, 4, 0, 4], + "octree_levels": None + }, + "interpolation_options": {} + } + + +def test_surface_points_loading(): + """Test loading surface points from JSON data.""" + data = { + "surface_points": [ + { + "x": 1.0, + "y": 1.0, + "z": 1.0, + "id": 0, + "nugget": 0.0 + } + ] + } + surface_points = JsonIO._load_surface_points(data["surface_points"]) + assert len(surface_points.xyz) == 1 + assert surface_points.xyz[0, 0] == 1.0 # x coordinate + assert surface_points.xyz[0, 1] == 1.0 # y coordinate + assert surface_points.xyz[0, 2] == 1.0 # z coordinate + assert surface_points.nugget[0] == 0.0 # nugget + + +def test_orientations_loading(): + """Test loading orientations from JSON data.""" + data = { + "orientations": [ + { + "x": 1.0, + "y": 1.0, + "z": 1.0, + "G_x": 0.0, + "G_y": 0.0, + "G_z": 1.0, + "id": 0, + "nugget": 0.01, + "polarity": 1 + } + ] + } + orientations = JsonIO._load_orientations(data["orientations"]) + assert len(orientations.xyz) == 1 + assert orientations.xyz[0, 0] == 1.0 # x coordinate + assert orientations.xyz[0, 1] == 1.0 # y coordinate + assert orientations.xyz[0, 2] == 1.0 # z coordinate + assert orientations.grads[0, 0] == 0.0 # G_x + assert orientations.grads[0, 1] == 0.0 # G_y + assert orientations.grads[0, 2] == 1.0 # G_z + assert orientations.nugget[0] == 0.01 # nugget + + +def test_surface_points_saving(tmp_path): + """Test saving surface points to a JSON file.""" + # Create sample surface points + surface_points = np.array([[1.0, 1.0, 1.0, 0, 0.0]]) + + # Create a complete model data structure + data = { + "surface_points": [ + { + "x": 1.0, + "y": 1.0, + "z": 1.0, + "id": 0, + "nugget": 0.0 + } + ], + "orientations": [], + "grid_settings": { + "regular_grid_resolution": [10, 10, 10], + "regular_grid_extent": [0, 10, 0, 10, 0, 10] + } + } + + # Save to temporary file + file_path = tmp_path / "test_surface_points.json" + with open(file_path, 'w') as f: + json.dump(data, f) + + # Load and verify + loaded_model = JsonIO.load_model_from_json(str(file_path)) + assert len(loaded_model.surface_points_copy.xyz) == 1 + assert loaded_model.surface_points_copy.xyz[0, 0] == 1.0 # x coordinate + assert loaded_model.surface_points_copy.xyz[0, 1] == 1.0 # y coordinate + assert loaded_model.surface_points_copy.xyz[0, 2] == 1.0 # z coordinate + assert loaded_model.surface_points_copy.nugget[0] == 0.0 # nugget + + +def test_invalid_surface_points_data(): + """Test that invalid surface points data raises appropriate errors.""" + data = { + 'surface_points': [{'x': 0, 'y': 0}], # Missing z + 'orientations': [], + 'grid_settings': {'regular_grid_resolution': [50, 50, 50], 'regular_grid_extent': [0, 1, 0, 1, 0, 1]} + } + with pytest.raises(ValueError, match="Missing required key in surface point: z"): + JsonIO._validate_json_schema(data) + + +def test_missing_surface_points_data(): + """Test handling of missing surface points data.""" + data = { + "orientations": [] # Missing surface_points key + } + with pytest.raises(ValueError, match="Missing required key: surface_points"): + JsonIO._validate_json_schema(data) + + +def test_invalid_orientations_data(): + """Test that invalid orientations data raises appropriate errors.""" + data = { + 'surface_points': [{'x': 0, 'y': 0, 'z': 0}], + 'orientations': [{'x': 0, 'y': 0, 'z': 0}], # Missing G_x, G_y, G_z + 'grid_settings': {'regular_grid_resolution': [50, 50, 50], 'regular_grid_extent': [0, 1, 0, 1, 0, 1]} + } + with pytest.raises(ValueError, match="Missing required key in orientation: G_x"): + JsonIO._validate_json_schema(data) + + +def test_missing_orientations_data(): + """Test handling of missing orientations data.""" + data = { + "surface_points": [] # Missing orientations key + } + with pytest.raises(ValueError, match="Missing required key: orientations"): + JsonIO._validate_json_schema(data) + + +def test_invalid_orientation_polarity(): + """Test handling of invalid orientation polarity.""" + data = { + "orientations": [ + { + "x": 1.0, + "y": 1.0, + "z": 1.0, + "G_x": 0.0, + "G_y": 0.0, + "G_z": 1.0, + "id": 0, + "nugget": 0.01, + "polarity": 2 # Invalid polarity value + } + ] + } + with pytest.raises(ValueError, match="Invalid polarity in orientation"): + JsonIO._load_orientations(data["orientations"]) + + +def test_default_nugget_values(tmp_path): + """Test that default nugget values are correctly applied.""" + # Create minimal data without nugget values + data = { + "surface_points": [ + { + "x": 1.0, + "y": 1.0, + "z": 1.0, + "id": 0 + } + ], + "orientations": [ + { + "x": 1.0, + "y": 1.0, + "z": 1.0, + "G_x": 0.0, + "G_y": 0.0, + "G_z": 1.0, + "id": 0, + "polarity": 1 + } + ], + "grid_settings": { + "regular_grid_resolution": [10, 10, 10], + "regular_grid_extent": [0, 10, 0, 10, 0, 10] + } + } + + # Save data to temporary file + file_path = tmp_path / "test_default_nugget.json" + with open(file_path, 'w') as f: + json.dump(data, f) + + # Load model from JSON + model = JsonIO.load_model_from_json(str(file_path)) + + # Verify default nugget values + assert model.surface_points_copy.nugget[0] == 0.0 # Default surface point nugget + assert model.orientations_copy.nugget[0] == 0.01 # Default orientation nugget + + +def test_default_series_values(tmp_path): + """Test that default series values are correctly applied.""" + # Create minimal data without series + data = { + "surface_points": [ + { + "x": 1.0, + "y": 1.0, + "z": 1.0, + "id": 0, + "nugget": 0.0 + } + ], + "orientations": [ + { + "x": 1.0, + "y": 1.0, + "z": 1.0, + "G_x": 0.0, + "G_y": 0.0, + "G_z": 1.0, + "id": 0, + "nugget": 0.01, + "polarity": 1 + } + ], + "grid_settings": { + "regular_grid_resolution": [10, 10, 10], + "regular_grid_extent": [0, 10, 0, 10, 0, 10] + } + } + + # Save data to temporary file + file_path = tmp_path / "test_default_series.json" + with open(file_path, 'w') as f: + json.dump(data, f) + + # Load model from JSON (this will create default series) + model = JsonIO.load_model_from_json(str(file_path)) + + # Verify default series + assert len(model.structural_frame.structural_groups) == 1 + assert model.structural_frame.structural_groups[0].name == "Strat_Series" + assert model.structural_frame.structural_groups[0].structural_relation == StackRelationType.ERODE + + +def test_default_interpolation_options(): + """Test that default interpolation options are correctly applied.""" + data = { + 'surface_points': [{'x': 0, 'y': 0, 'z': 0}], + 'orientations': [{'x': 0, 'y': 0, 'z': 0, 'G_x': 0, 'G_y': 0, 'G_z': 1}], + 'grid_settings': {'regular_grid_resolution': [50, 50, 50], 'regular_grid_extent': [0, 1, 0, 1, 0, 1]}, + 'interpolation_options': {} + } + JsonIO._validate_json_schema(data) + assert data['interpolation_options'].get('number_octree_levels', 4) == 4 # Default is now 4 + assert data['interpolation_options'].get('mesh_extraction', True) is True + + +def test_default_metadata(tmp_path): + """Test that default metadata is correctly applied.""" + # Create minimal data without metadata + data = { + "surface_points": [ + { + "x": 1.0, + "y": 1.0, + "z": 1.0, + "id": 0, + "nugget": 0.0 + } + ], + "orientations": [ + { + "x": 1.0, + "y": 1.0, + "z": 1.0, + "G_x": 0.0, + "G_y": 0.0, + "G_z": 1.0, + "id": 0, + "nugget": 0.01, + "polarity": 1 + } + ], + "grid_settings": { + "regular_grid_resolution": [10, 10, 10], + "regular_grid_extent": [0, 10, 0, 10, 0, 10] + } + } + + # Save data to temporary file + file_path = tmp_path / "test_default_metadata.json" + with open(file_path, 'w') as f: + json.dump(data, f) + + # Load model from JSON + model = JsonIO.load_model_from_json(str(file_path)) + + # Verify default metadata + assert model.meta.name == "GemPy Model" # Default name + assert model.meta.owner == "GemPy Modeller" # Default owner + assert model.meta.creation_date is not None # Should be set to current date + assert model.meta.last_modification_date is not None # Should be set to current date + + +def test_optional_series_colors(): + """Test that series colors are optional and properly validated.""" + data = { + 'surface_points': [{'x': 0, 'y': 0, 'z': 0}], + 'orientations': [{'x': 0, 'y': 0, 'z': 0, 'G_x': 0, 'G_y': 0, 'G_z': 1}], + 'grid_settings': {'regular_grid_resolution': [50, 50, 50], 'regular_grid_extent': [0, 1, 0, 1, 0, 1]}, + 'series': [ + {'name': 'Series1', 'surfaces': ['Surface1'], 'colors': ['#015482']}, # Valid color + {'name': 'Series2', 'surfaces': ['Surface2']} # No colors specified + ] + } + JsonIO._validate_json_schema(data) + assert data['series'][0]['colors'] == ['#015482'] # Color should be preserved + assert 'colors' not in data['series'][1] # No colors should be added if not specified + + +@pytest.fixture +def sample_model_with_series(): + """Create a sample model with series and fault relations.""" + # Create a simple model with two series and a fault + from gempy.core.data.structural_frame import StructuralFrame + from gempy.core.data import StructuralElement, StructuralGroup + from gempy_engine.core.data.stack_relation_type import StackRelationType + + # Create structural frame first + structural_frame = StructuralFrame.initialize_default_structure() + + # Create model with structural frame + model = gp.create_geomodel( + project_name="test_model", + extent=[0, 10, 0, 10, 0, 10], + resolution=[10, 10, 10], + structural_frame=structural_frame + ) + + # Create structural elements + fault_element = StructuralElement( + name="fault", + color=next(model.structural_frame.color_generator), + surface_points=gp.data.SurfacePointsTable.initialize_empty(), + orientations=gp.data.OrientationsTable.initialize_empty() + ) + + rock1_element = StructuralElement( + name="rock1", + color=next(model.structural_frame.color_generator), + surface_points=gp.data.SurfacePointsTable.initialize_empty(), + orientations=gp.data.OrientationsTable.initialize_empty() + ) + + rock2_element = StructuralElement( + name="rock2", + color=next(model.structural_frame.color_generator), + surface_points=gp.data.SurfacePointsTable.initialize_empty(), + orientations=gp.data.OrientationsTable.initialize_empty() + ) + + # Create structural groups + fault_group = StructuralGroup( + name="fault_series", + elements=[fault_element], + structural_relation=StackRelationType.FAULT + ) + + strat_group = StructuralGroup( + name="strat_series", + elements=[rock1_element, rock2_element], + structural_relation=StackRelationType.ERODE + ) + + # Set up the structural frame + model.structural_frame.structural_groups = [fault_group, strat_group] + + # Set fault relations (2x2 matrix: fault affects strat_series) + model.structural_frame.fault_relations = np.array([ + [0, 1], # fault_series affects strat_series + [0, 0] # strat_series doesn't affect any series + ]) + + # Add surface points + gp.add_surface_points( + geo_model=model, + x=[1, 2, 3, 4, 5], + y=[1, 2, 3, 4, 5], + z=[1, 2, 3, 4, 5], + elements_names=["fault", "rock1", "rock1", "rock2", "rock2"] + ) + + # Add orientations + gp.add_orientations( + geo_model=model, + x=[1.5, 2.5, 3.5, 4.5], + y=[1.5, 2.5, 3.5, 4.5], + z=[1.5, 2.5, 3.5, 4.5], + elements_names=["fault", "rock1", "rock1", "rock2"], + pole_vector=[[0, 0, 1], [0, 0, 1], [0, 0, 1], [0, 0, 1]] + ) + + return model + + +def test_metadata_handling(tmp_path, sample_model_with_series): + """Test metadata handling in JSON I/O.""" + model = sample_model_with_series + + # Add custom metadata + model.meta.owner = "test_value" + + # Save model + file_path = tmp_path / "test_metadata.json" + JsonIO.save_model_to_json(model, str(file_path)) + + # Load model + loaded_model = JsonIO.load_model_from_json(str(file_path)) + + # Verify metadata + assert loaded_model.meta.owner == "test_value" + assert loaded_model.meta.creation_date is not None + assert loaded_model.meta.last_modification_date is not None + + +def test_grid_settings_handling(tmp_path, sample_model_with_series): + """Test grid settings handling in JSON I/O.""" + model = sample_model_with_series + + # Modify grid settings + model.grid.regular_grid.resolution = [20, 20, 20] + model.grid.regular_grid.extent = [0, 20, 0, 20, 0, 20] + + # Save model + file_path = tmp_path / "test_grid_settings.json" + JsonIO.save_model_to_json(model, str(file_path)) + + # Load model + loaded_model = JsonIO.load_model_from_json(str(file_path)) + + # Verify grid settings + np.testing.assert_array_equal( + model.grid.regular_grid.resolution, + loaded_model.grid.regular_grid.resolution + ) + np.testing.assert_array_equal( + model.grid.regular_grid.extent, + loaded_model.grid.regular_grid.extent + ) + + +def test_interpolation_options_handling(tmp_path, sample_model_with_series): + """Test interpolation options handling in JSON I/O.""" + model = sample_model_with_series + + # Set custom interpolation options + model.interpolation_options.kernel_options.range = 10.0 + model.interpolation_options.kernel_options.c_o = 15.0 + model.interpolation_options.mesh_extraction = False + model.interpolation_options.number_octree_levels = 2 + + # Save model + file_path = tmp_path / "test_interpolation_options.json" + JsonIO.save_model_to_json(model, str(file_path)) + + # Load model + loaded_model = JsonIO.load_model_from_json(str(file_path)) + + # Verify interpolation options + assert loaded_model.interpolation_options.kernel_options.range == 10.0 + assert loaded_model.interpolation_options.kernel_options.c_o == 15.0 + assert loaded_model.interpolation_options.mesh_extraction is False + assert loaded_model.interpolation_options.number_octree_levels == 2 + + +def test_fault_relationships(tmp_path): + """Test saving and loading models with fault relationships.""" + # Create temporary CSV files for surface points and orientations + surface_points_data = pd.DataFrame({ + 'X': [500, 500, 500], + 'Y': [500, 500, 500], + 'Z': [800, 500, 200], + 'surface': ['fault', 'rock2', 'rock1'] + }) + orientations_data = pd.DataFrame({ + 'X': [500, 500, 500], + 'Y': [500, 500, 500], + 'Z': [800, 500, 200], + 'dip': [90, 0, 0], + 'azimuth': [90, 90, 90], + 'polarity': [1, 1, 1], + 'surface': ['fault', 'rock2', 'rock1'] + }) + + surface_points_path = tmp_path / "surface_points.csv" + orientations_path = tmp_path / "orientations.csv" + surface_points_data.to_csv(surface_points_path, index=False) + orientations_data.to_csv(orientations_path, index=False) + + # Create a model with fault relationships + importer_helper = ImporterHelper( + path_to_surface_points=str(surface_points_path), + path_to_orientations=str(orientations_path) + ) + model = gp.create_geomodel( + project_name='fault_test', + extent=[0, 1000, 0, 1000, 0, 1000], + resolution=[10, 10, 10], + importer_helper=importer_helper + ) + + # Map stack to surfaces + gp.map_stack_to_surfaces( + gempy_model=model, + mapping_object={ + "Fault_Series": ['fault'], + "Strat_Series1": ['rock2'], + "Strat_Series2": ['rock1'] + } + ) + + # Set structural relations + model.structural_frame.structural_groups[0].structural_relation = StackRelationType.FAULT + model.structural_frame.structural_groups[0].fault_relations = FaultsRelationSpecialCase.OFFSET_ALL + + # Save the model + file_path = tmp_path / "fault_model.json" + JsonIO.save_model_to_json(model, str(file_path)) + + # Load the model and verify + loaded_model = JsonIO.load_model_from_json(str(file_path)) + + # Check structural groups + assert len(loaded_model.structural_frame.structural_groups) == 3 + assert loaded_model.structural_frame.structural_groups[0].name == "Fault_Series" + assert loaded_model.structural_frame.structural_groups[0].structural_relation == StackRelationType.FAULT + assert loaded_model.structural_frame.structural_groups[0].fault_relations == FaultsRelationSpecialCase.OFFSET_ALL + + # Check fault relations matrix + np.testing.assert_array_equal( + loaded_model.structural_frame.fault_relations, + np.array([[0, 1, 1], [0, 0, 0], [0, 0, 0]]) + ) + + +def test_multiple_series_relationships(tmp_path): + """Test saving and loading models with multiple series and different structural relationships.""" + # Create temporary CSV files for surface points and orientations + surface_points_data = pd.DataFrame({ + 'X': [500, 500, 500, 500], + 'Y': [500, 500, 500, 500], + 'Z': [800, 600, 400, 200], + 'surface': ['fault', 'rock3', 'rock2', 'rock1'] + }) + orientations_data = pd.DataFrame({ + 'X': [500, 500, 500, 500], + 'Y': [500, 500, 500, 500], + 'Z': [800, 600, 400, 200], + 'dip': [90, 0, 0, 0], + 'azimuth': [90, 90, 90, 90], + 'polarity': [1, 1, 1, 1], + 'surface': ['fault', 'rock3', 'rock2', 'rock1'] + }) + + surface_points_path = tmp_path / "surface_points.csv" + orientations_path = tmp_path / "orientations.csv" + surface_points_data.to_csv(surface_points_path, index=False) + orientations_data.to_csv(orientations_path, index=False) + + # Create a model with multiple series relationships + importer_helper = ImporterHelper( + path_to_surface_points=str(surface_points_path), + path_to_orientations=str(orientations_path) + ) + model = gp.create_geomodel( + project_name='multiple_series_test', + extent=[0, 1000, 0, 1000, 0, 1000], + resolution=[10, 10, 10], + importer_helper=importer_helper + ) + + # Map stack to surfaces + gp.map_stack_to_surfaces( + gempy_model=model, + mapping_object={ + "Fault_Series": ['fault'], + "Series3": ['rock3'], + "Series2": ['rock2'], + "Series1": ['rock1'] + } + ) + + # Set structural relations + model.structural_frame.structural_groups[0].structural_relation = StackRelationType.FAULT + model.structural_frame.structural_groups[0].fault_relations = FaultsRelationSpecialCase.OFFSET_ALL + model.structural_frame.structural_groups[1].structural_relation = StackRelationType.ONLAP + model.structural_frame.structural_groups[2].structural_relation = StackRelationType.ERODE + model.structural_frame.structural_groups[3].structural_relation = StackRelationType.ERODE + + # Set colors for verification + model.structural_frame.get_element_by_name("fault").color = '#527682' # Fault color + model.structural_frame.get_element_by_name("rock3").color = '#9f0052' # rock3 color + model.structural_frame.get_element_by_name("rock2").color = '#015482' # rock2 color + model.structural_frame.get_element_by_name("rock1").color = '#728f02' # rock1 color + + # Save the model + file_path = tmp_path / "multiple_series_model.json" + JsonIO.save_model_to_json(model, str(file_path)) + + # Load the model and verify + loaded_model = JsonIO.load_model_from_json(str(file_path)) + + # Check structural groups + assert len(loaded_model.structural_frame.structural_groups) == 4 + assert loaded_model.structural_frame.structural_groups[0].structural_relation == StackRelationType.FAULT + assert loaded_model.structural_frame.structural_groups[1].structural_relation == StackRelationType.ONLAP + assert loaded_model.structural_frame.structural_groups[2].structural_relation == StackRelationType.ERODE + + # Check colors are preserved + assert loaded_model.structural_frame.get_element_by_name("fault").color == '#527682' + assert loaded_model.structural_frame.get_element_by_name("rock3").color == '#9f0052' + assert loaded_model.structural_frame.get_element_by_name("rock2").color == '#015482' + assert loaded_model.structural_frame.get_element_by_name("rock1").color == '#728f02' + + +def test_combination_model(tmp_path): + """Test saving and loading a complex model that combines faults and unconformities.""" + # Create temporary CSV files for surface points and orientations + surface_points_data = pd.DataFrame({ + 'X': [500, 500, 500], + 'Y': [500, 500, 500], + 'Z': [800, 500, 200], + 'surface': ['fault', 'rock2', 'rock1'] + }) + orientations_data = pd.DataFrame({ + 'X': [500, 500, 500], + 'Y': [500, 500, 500], + 'Z': [800, 500, 200], + 'dip': [90, 0, 0], + 'azimuth': [90, 90, 90], + 'polarity': [1, 1, 1], + 'surface': ['fault', 'rock2', 'rock1'] + }) + + surface_points_path = tmp_path / "surface_points.csv" + orientations_path = tmp_path / "orientations.csv" + surface_points_data.to_csv(surface_points_path, index=False) + orientations_data.to_csv(orientations_path, index=False) + + # Create a combination model + importer_helper = ImporterHelper( + path_to_surface_points=str(surface_points_path), + path_to_orientations=str(orientations_path) + ) + model = gp.create_geomodel( + project_name='Combination Model', + extent=[0, 1000, 0, 1000, 0, 1000], + resolution=[10, 10, 10], + importer_helper=importer_helper + ) + + # Map stack to surfaces + gp.map_stack_to_surfaces( + gempy_model=model, + mapping_object={ + "Fault_Series": ['fault'], + "Strat_Series1": ['rock2'], + "Strat_Series2": ['rock1'] + } + ) + + # Set structural relations + model.structural_frame.structural_groups[0].structural_relation = StackRelationType.FAULT + model.structural_frame.structural_groups[0].fault_relations = FaultsRelationSpecialCase.OFFSET_ALL + model.structural_frame.structural_groups[1].structural_relation = StackRelationType.ERODE + model.structural_frame.structural_groups[2].structural_relation = StackRelationType.ONLAP + + # Set metadata + model.meta.creation_date = "2024-03-24" + model.meta.last_modification_date = "2024-03-24" + + # Save the model + file_path = tmp_path / "combination_model.json" + JsonIO.save_model_to_json(model, str(file_path)) + + # Load the model and verify + loaded_model = JsonIO.load_model_from_json(str(file_path)) + + # Check structural setup + assert len(loaded_model.structural_frame.structural_groups) == 3 + assert loaded_model.structural_frame.structural_groups[0].structural_relation == StackRelationType.FAULT + assert loaded_model.structural_frame.structural_groups[1].structural_relation == StackRelationType.ERODE + assert loaded_model.structural_frame.structural_groups[2].structural_relation == StackRelationType.ONLAP + + # Check fault relations + np.testing.assert_array_equal( + loaded_model.structural_frame.fault_relations, + np.array([[0, 1, 1], [0, 0, 0], [0, 0, 0]]) + ) + + # Check metadata + assert loaded_model.meta.creation_date == "2024-03-24" + assert loaded_model.meta.last_modification_date == "2024-03-24" \ No newline at end of file