import geosoft.gxpy.gx as gx
import geosoft.gxpy.view as gxview
import geosoft.gxpy.group as gxgroup
import geosoft.gxpy.vox as gxvox
import geosoft.gxpy.vox_display as gxvoxd
import geosoft.gxpy.viewer as gxviewer
import geosoft.gxpy.utility as gxu
import geosoft.gxpy.map as gxmap
import geosoft.gxpy.surface as gxsurf
from IPython.display import Image
import numpy as np

gxc = gx.GXpy()

url = 'https://github.com/GeosoftInc/gxpy/raw/9.3.1/examples/data/'
vox_file = 'rjsmith_voxi_density.geosoft_voxel'
gxu.url_retrieve(url + vox_file)
gxu.url_retrieve(url + vox_file + '.xml')
vectorvox_file = 'mvi.geosoft_vectorvoxel'
gxu.url_retrieve(url + vectorvox_file)
gxu.url_retrieve(url + vectorvox_file + '.xml')

"""
Explanation: Copyright (c) 2018 Geosoft Inc.
https://github.com/GeosoftInc/gxpy
BSD 2-clause License
Geosoft Voxels
An important geoscience data structure is a voxel (voxset), which is a 3D cell-based model of the subsurface volume of part of the Earth. Voxel models are commonly used to represent some property of the subsurface in which the entire volume of interest is subdivided into cells, each with a value that represents a property of the cell volume. A voxel property can also be a vector quantity, in which case the property of a voxel cell will have a vector quantity defined by (vx, vy, vz).
We work with voxels using the Vox class and the Vox_display class.  The Vox class can be used to open a Geosoft voxel (or vector voxel) or create a new Geosoft voxel, and the Vox_display class can create a visualization of a voxel that can be placed in a 3D view.
Lessons
<!--- # Run this from a code cell to create TOC markdown: -->
<!--- import geosoft.gxpy.utility; print(geosoft.gxpy.utility.jupyter_markdown_toc('geosoft voxel')) -->

Python imports, GX context and get exercise data from GitHub
Voxel properties
Referenced by elevation, or by depth
View the voxel as a solid
View the voxel as isosurfaces
New voxel from a numpy data array
Locate a new voxel spatially
Define a variable cell size
Vector Voxel
Display a vector voxel as vector cones

Python imports, GX context and get exercise data from GitHub
End of explanation
"""


# open the vox
with gxvox.Vox.open(vox_file) as vox:
    
    # show some of the vox properties
    print(vox.name)
    print('dimensions (nx, ny, nz):', vox.nx, vox.ny, vox.nz)
    print('coordinate system:', vox.coordinate_system)
    print('surface extent:', vox.extent_xy)
    print('volume extent:', vox.extent)
    print('origin location:', (vox.origin_x, vox.origin_y, vox.origin_z))
    print('unit of measure:', vox.unit_of_measure)
    print('depth cell sizes in', vox.coordinate_system.unit_of_measure, ', from the bottom:')
    for cell in vox.cells_z:
        print('\t', cell)

"""
Explanation: Voxel properties
End of explanation
"""


# open the vox
with gxvox.Vox.open(vox_file, depth=True) as vox:
    
    # show some of the vox properties
    print(vox.name, 'properties by depth...')
    print('volume extent:', vox.extent)
    print('origin location:', (vox.origin_x, vox.origin_y, vox.origin_z))
    print('depth cell sizes in', vox.coordinate_system.unit_of_measure, ', from the top:')
    for cell in vox.cells_z:
        print('\t', cell)

"""
Explanation: Referenced by elevation, or by depth
The vertical dimension of a voxel (the 'z' dimension), can be referenced by elevation or by depth.  If referenced by elevation, which is the default, the positive z direction is up and values can be assumed to represent elevations relative to the vertical reference datum of the coordinate system. Also, when referenced by elevation, the voxel origin is on the bottom slice of the voxel and z arrays of voxel cells sizes or locations start at the bottom and sequence up, in the positive direction.
You can also choose to work with a voxel by depth, either by setting the depth=True parameter when a new voxel is created, or an existing voxel is openned, or by changing the is_depth property to True.
For example, here we change the reference to depth and print depth-based properties. Note the z property differences from the example above.
End of explanation
"""


# open the vox
with gxvox.Vox.open(vox_file) as vox:
    
    # create a display instance
    with gxvoxd.VoxDisplay.solid(vox) as voxd:
        
        # limit the display to densities above 0.005 g/cc
        voxd.shell_limits = (0.005, None)
        
        # Create a title and include the coordinate system.
        title = "VOXI Density Model - RJ Smith Test Range" + '\n' + str(vox.coordinate_system)
        
        # create a figure map for the vox. This creates a 3D view that is placed on a map and annotated 
        # to create a default presentation figure. By not specifying a name for the map a temporary map
        # is created only for the life of this Python session.
        vox_map_file = voxd.figure_map(title=title, legend_label='density\ng/cc').file_name

# create a PNG from the figure map and display in Jupyter
Image(gxmap.Map.open(vox_map_file).image_file(pix_width=800))

"""
Explanation: View the voxel as a solid
A Vox instance can be displayed withion a 3D view either as a solid with each cell assigned a color, or as a set of iso-surfaces. Vector voxels can be displayed as a solid, as vectors or as iso-surfaces.
Here we create a solid with starting minimum value of 0.005, though this can be changed dynamically in the Geosoft 3D viewer. The figure_map() method will render the solid into a group of a 3d view in a 3d map.
End of explanation
"""


# create an isosurface dataset, shell at 0.005
isosurface_file = gxsurf.SurfaceDataset.vox_surface(gxvox.Vox.open(vox_file), 0.005, overwrite=True).file_name

# create a temporary 3d view and place the isosurface dataset in the view
with gxview.View_3d.new() as v3d:
    v3d_file = v3d.file_name
    gxgroup.surface_group_from_file(v3d, isosurface_file)
    
# display here
Image(gxmap.Map.open(v3d_file).image_file(pix_width=800))

"""
Explanation: View the voxel as isosurfaces
To view the voxel as an isosurface we first create a SurfaceDataset from the voxel, which we then add to a 3d view.
End of explanation
"""


# create an isosurface dataset, shells at 0.002, 0.01 and 0.05
isosurface_file = gxsurf.SurfaceDataset.vox_surface(gxvox.Vox.open(vox_file), 
                                                    surfaces=(0.002, 0.01, 0.05),
                                                    color=(gxgroup.C_YELLOW, gxgroup.C_GREEN, gxgroup.C_RED),
                                                    overwrite=True).file_name

# create a temporary 3d view and place the isosurface dataset in the view
with gxview.View_3d.new() as v3d:
    v3d_file = v3d.file_name
    gxgroup.surface_group_from_file(v3d, isosurface_file)
    
# display here
Image(gxmap.Map.open(v3d_file).image_file(pix_width=800))

"""
Explanation: Here we create three isosurfaces with specified colors.
End of explanation
"""


# create a 3D data array in which each z slice is set to the z index
data = np.empty((20, 65, 100), dtype=np.float64)
for i in range(data.shape[0]):
    data[i, :, :] = i
    
# create a new voxel, and display
with gxvox.Vox.new(name='my vox', data=data, overwrite=True) as vox:
    vox_map_file = gxvoxd.VoxDisplay.solid(vox).figure_map().file_name

# create a PNG from the figure map and display in Jupyter
Image(gxmap.Map.open(vox_map_file).image_file(pix_width=800))

"""
Explanation: New voxel from a numpy data array
A simple way to create a new voxel is to collect the voxel data into a 3D numpy array in which each dimension represents the cells in the X, Y and Z spatial directions. Data in an array is organized by X, then by Y, and lastly by Z, which results in data indexing being in the reverse order.  For example, a voxel that has spatial dimensions (nx, ny, nz) of (100, 65, 20) will be stored in a 3d array with shape (20, 65, 100).
Note also that the Z direction is positive up, which means that a voxel origin, which is the center of the first cell (cell[0, 0, 0]) is the at the center of the deepest, lower-left cell in the voxel.
In this example we will start with a simple indexed data array in which the value of each cell is the Z index.
End of explanation
"""


# create spatially located voxel
with gxvox.Vox.new(name='my vox', 
                   data=data, 
                   overwrite=True, 
                   origin=(450025, 6800025, -975), # the origin is the center of the bottom cell[0, 0, 0]
                   cell_size=50, 
                   coordinate_system='NAD83 / UTM zone 19N') as vox:
    
    # display in a map
    vox_map_file = gxvoxd.VoxDisplay.solid(vox).figure_map(title=str(vox.coordinate_system)).file_name

# display the map here
Image(gxmap.Map.open(vox_map_file).image_file(pix_width=800))

"""
Explanation: Locate a new voxel spatially
In the last example we created a voxel without true spatial dimensions. In this example we will scale the same voxel to have true spatial dimensions and be located at a known location on Earth. For this example we will assume that the voxels are 50 m square cells and the model is located on the NAD83 datum using projection UTM zone 19N.
End of explanation
"""


# create a z cell dimension increasing 10% as a function of depth (positive down)
nz = data.shape[0]
cell = 50
z_cell = [cell] * nz
for i in range(1, nz):
    z_cell[i] = cell * (1. + i * 0.1)
    
# create voxel, and display
with gxvox.Vox.new(name='my vox', 
                   data=data, 
                   overwrite=True,
                   depth=True, # work using z as depth
                   origin=(450025, 6800025, 25), # the origin is the center of the top cell when working from depth 
                   cell_size=(cell, cell, z_cell), # the z_cell size array defines variable z cell sizes
                   coordinate_system='NAD83 / UTM zone 19N') as vox:
    vox_map_file = gxvoxd.VoxDisplay.solid(vox).figure_map(title=str(vox.coordinate_system)).file_name

# display here
Image(gxmap.Map.open(vox_map_file).image_file(pix_width=800))

"""
Explanation: Define a variable cell size
Voxel cell sizes may vary along any dimension. For example, when defining a voxel model for geophysical modelling one might increase the vertical size of voxel as a function of depth.  To do this you can provide the cell dimensions as an array of dimensions along any, or all axis. In this example we will increase the cell size as a function of depth, with each deeper cell being 10% larger. We will referencing depth rather than elevation as it makes the cell depth calculation simpler.
End of explanation
"""


# open the vector voxel
with gxvox.Vox.open(vectorvox_file) as vox:
    
    # create a solid display
    with gxvoxd.VoxDisplay.solid(vox) as voxd:
        
        # limit the display to MVI susceptibility above 0.01 SI
        voxd.shell_limits = (0.04, None)
        
        # Create a title and include the coordinate system.
        title = "VOXI MVI Susceptibility" + '\n' + str(vox.coordinate_system)
        
        # create a figure map for the vox.
        vox_map_file = voxd.figure_map(title=title, legend_label='MVI SI').file_name

# create a PNG from the figure map and display in Jupyter
Image(gxmap.Map.open(vox_map_file).image_file(pix_width=800))

"""
Explanation: Vector Voxel
Voxel cells may also store a vector value defined by a (vx, vy, vz) for each cell. For example, Geosoft MVI models define model magnetization as a vector quantity. Vector voxel files have extension .geosoft_vectorvoxel.
The is_vectorvox property of a Vox instance that holds a vector voxel will be True, and the Numpy array returned by the Vox.np() method will have 4 dimensions with the last dimension being 3 which contains the (vx, vy, vz) vector for each voxel cell.
Displaying a vector voxel as a solid will display the Vox as scalar colour based on vector amplitude, as follows.
End of explanation
"""


# open the vector voxel
with gxvox.Vox.open(vectorvox_file) as vox:
    
    # create a vector display
    with gxvoxd.VoxDisplay.vector(vox) as voxd:
        
        # limit the display to MVI susceptibility above 0.04 SI
        voxd.shell_limits = (0.04, None)
        
        # set the vector specifications
        voxd.vector_cone_specs = (3, None, None, 1000)
        
        # Create a title and include the coordinate system.
        title = "VOXI MVI Susceptibility"
        
        # create a figure map for the vox.
        vox_map_file = voxd.figure_map(title=title, legend_label='MVI SI').file_name

# create a PNG from the figure map and display in Jupyter
Image(gxmap.Map.open(vox_map_file).image_file(pix_width=800))

"""
Explanation: Display a vector voxel as vector cones
Vector voxels can also be displayed as cones that provide a visualization of the vector direction.  In this case the cone indicates the direction of magnetization in the MVI model.
End of explanation
"""


