# These two lines are necessary only if gempy is not installed
import sys, os
sys.path.append("../")

# Importing gempy
import gempy as gp

# Embedding matplotlib figures into the notebooks
%matplotlib inline

# Aux imports
import numpy as np

"""
Explanation: Chapter 2: A real example. Importing data and setting series
Data Management
In this example we will show how we can import data from a csv and generate a model with several depositional series.
End of explanation
"""


# Importing the data from csv files and settign extent and resolution
geo_data = gp.create_data([696000,747000,6863000,6950000,-20000, 200],[50, 50, 50],
                         path_f = os.pardir+"/input_data/a_Foliations.csv",
                         path_i = os.pardir+"/input_data/a_Points.csv")

gp.get_raw_data(geo_data, 'interfaces').head()

"""
Explanation: In this case instead loading a geo_data object directly, we will create one. The main atributes we need to pass are:
- Extent: X min, X max, Y min, Y max, Z min, Z max
- Resolution: X,Y,Z
Additionaly we can pass the address to csv files (GeoModeller3D format) with the data.
End of explanation
"""


gp.visualize(geo_data)

"""
Explanation: You can visualize the points in 3D (work in progress)
End of explanation
"""


gp.plot_data(geo_data, direction='z')

"""
Explanation: Or a projection in 2D:
End of explanation
"""


gp.get_series(geo_data)

"""
Explanation: This model consist in 3 different depositional series. This mean that only data in the same depositional series affect the interpolation. To select with formations belong to witch series we will use the set_data_series function which takes a python dictionary as input.
We can see the unique formations with:
End of explanation
"""


# Assigning series to formations as well as their order (timewise)
gp.set_data_series(geo_data, {"EarlyGranite_Series": 'EarlyGranite', 
                              "BIF_Series":('SimpleMafic2', 'SimpleBIF'),
                              "SimpleMafic_Series":'SimpleMafic1'}, 
                      order_series = ["EarlyGranite_Series",
                                      "BIF_Series",
                                      "SimpleMafic_Series"], verbose=1)

"""
Explanation: Setting the series we also give the specific order of the series. In python 3.6 and above the dictionaries conserve the key order so it is not necessary to give explicitly the order of the series. 
Notice as well that the order of the formations within each series is not relevant for the result but in case of being wrong can lead to confusing color coding (work in progress). 
In the representation given by get_series the elements get repeated but is only how Pandas print tables.
End of explanation
"""


interp_data = gp.InterpolatorInput(geo_data)

sol = gp.compute_model(interp_data)

"""
Explanation: Computing the model
Now as in the previous chapter we just need to create the interpolator object and compute the model.
End of explanation
"""


import matplotlib.pyplot as plt
gp.plot_section(geo_data, sol[0,0,:], 11)
plt.show()
gp.plot_section(geo_data, sol[1,0,:], 11)
plt.show()
gp.plot_section(geo_data, sol[2,0,:], 11)
plt.show()

"""
Explanation: Now if we analyse the results we have a 3D array where the axis 0 represent the superposition of the series (potential fields). The color coding is working process yet.
End of explanation
"""


gp.plot_potential_field(geo_data, sol[0,1,:], 11, cmap='inferno_r')
plt.show()
gp.plot_potential_field(geo_data, sol[1,1,:], 11, cmap='inferno_r')
plt.show()
gp.plot_potential_field(geo_data, sol[2,1,:], 11, cmap='inferno_r')
plt.show()

"""
Explanation: The axis 1 keeps the potential field:
End of explanation
"""


gp.export_vtk_rectilinear(geo_data, sol[-1, 0, :], path=None)

"""
Explanation: And the axis 2 keeps the faults network that in this model since there is not faults does not represent anything.
Additionally with can export the blocks to vtk in order to visualize them in Paraview. We are working in visualization in place as well.
End of explanation
"""


