In [1]:
# load some modules
import pandas as pd
import matplotlib.pylab as plt
import numpy as np
import pygslib
In [2]:
# see block model help
help(pygslib.blockmodel)
Help on module pygslib.blockmodel in pygslib:
NAME
pygslib.blockmodel - PyGSLIB Blockmodel, Module to handle blockmodel data.
FILE
c:\og_python\pygslib\pygslib\blockmodel.pyd
DESCRIPTION
Copyright (C) 2015 Adrian Martinez Vargas
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 3 of the License, or
any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>
CLASSES
__builtin__.object
Blockmodel
class Blockmodel(__builtin__.object)
| Blockmodel object.
|
|
|
| Parameters
| ----------
| nx,ny,nz : integer
| Number of rows columns and levels
| xorg,yorg,zorg : float
| Coordinates of the lower left corner (not centroid)
| dx,dy,dz : float
| Size of the parent block
|
| Attributes
| ----------
| nx, ny, nz : int
| number of rows, columns, levels
|
| xorg, yorg, zorg : float
| coordinates of the left lower corner of the first block
|
| dx, dy, dz : float
| sizes of the blocks
|
| bmtable : Pandas DataFrame
| Table with blocks
|
| gridtable : Pandas Dataframe
|
| Methods defined here:
|
| block2point(...)
| block2point(np.ndarray [double, ndim=1] x, np.ndarray [double, ndim=1] y, np.ndarray [double, ndim=1] z, str prop_name)
|
| Assigns a block property to the points inside a block
|
|
| Parameters
| ----------
| x,y,z : float 1D numpy arrays
| coordinates of the points. Array shapes may be equal
| prop_name : str
| name of the property at block model
|
|
| Example
| -------
| >>>
| >>> Au_at_point= point2block( x, y, z, prop_name= 'Au')
| >>>
|
|
| Note
| ----
| Points not laying in an existing block will be returned with
| value numpy.nan
|
| blockinsurface(...)
| blockinsurface(object surface, str field, double azm=0, double dip =90, int test=1, bint overwrite=False)
|
| Creates blocks given a VTK surface (polydata) depending on
| test criteria:
|
| Parameters
| ----------
| surface : VTK polydata
| this may work with any 3D object..., no necessarily triangles
| field : str
| Name of the new field with selection results
| azm, dip: float, default 0, 90
| rotation defining the direction we will use to test the points
| azm 0 will point north and dip positive means downward direction
| (like surface drillholes)
| test : integer, default 1
| ``test == 1`` test inside closed surface. Here we use
| vtkOBBTree::InsideOrOutside. Closed surface are required
| ``test == 2`` test 'above' surface
| ``test == 3`` test 'below' surface
| ``test == 4`` test 'inside' surface (the surface can be open)
| overwrite : boolean
| overwrite flag, if true and field exist in bmtable the
| values will be overwritten
|
| Example
| -------
| >>>
| >>> fillblocks(surface, azm, dip, test, overwrite=False)
| >>>
|
| Note
| ----
| This function calls vtktools.pointquering for all the points
| existing in the block model. The function is not optimized
| for large model.
|
| A new set of blocks will be created if there is not block
| defined. If there are blocks in the model and overwrite==True
| the blocks will be removed first.
|
| blocks2vtkRectilinearGrid(...)
| blocks2vtkRectilinearGrid(str path)
|
| Exports blocks of a full grid to a vtkRectilinearGrid file.
|
| Parameters
| ----------
| path : string
| file name and path, without extension. The file extension
| (*.vtr) will be added automatically.
|
| Examples
| --------
| >>>
| >>> blocks2vtkRectilinearGrid('myfile')
| >>>
|
| Note
| ----
| This will only work for full grid, in other words, if all the
| nx*ny*nz are defined.
|
| All the fields defined in the block model will be exported
|
| blocks2vtkUnstructuredGrid(...)
| blocks2vtkUnstructuredGrid(str path, str varname = None)
|
| Exports blocks of a partial grid to a vtkRectilinearGrid file.
|
| Parameters
| ----------
| path : string
| file name and path, without extension. The file extension
| (*.vtr) will be added automatically.
|
| varname : string
| If None export all columns, otherwise export only
| the column varname
|
|
| Examples
| --------
| >>>
| >>> blocks2vtkUnstructuredGrid('myfile', 'AU')
| >>>
|
| Note
| ----
| Require XC, YC and ZC.
|
| Only numeric variables can be exported
| TODO: export all
|
| This call pygslib.vtktools.partialgrid2vtkfile
|
| calc_ijk(...)
| calc_ijk(bint overwrite=False)
|
| Calculates the IJK field from IX, IY, IZ indices
| If IJK exist and overwrite=True the existing values
| will be overwritten.
|
|
| Parameters
| ----------
| overwrite : Boolean, default False
| If True overwrites any existing IJK field
|
| Example
| -------
| >>>
| >>> myblockmodel.calc_ijk()
| >>> myblockmodel.calc_ijk(overwrite=True)
| >>>
|
| calc_ixyz_fromijk(...)
| calc_ixyz_fromijk(bint overwrite=False)
|
| Calculates the IX, IY, IZ fields from IJK index
| If IX, IY, IZ exist and overwrite=True the existing values
| will be overwritten.
|
|
| Parameters
| ----------
| overwrite : Boolean, default False
|
| Example
| -------
| >>>
| >>> myblockmodel.calc_ixyz_fromijk()
| >>> myblockmodel.calc_ixyz_fromijk(overwrite=True)
| >>>
|
| calc_ixyz_fromxyz(...)
| calc_ixyz_fromxyz(bint overwrite=False)
|
| Calculates the IX, IY, IZ fields from XC, YC, ZC coordinates
| If IX, IY, IZ exist and overwrite=True the existing values
| will be overwritten.
|
|
| Parameters
| ----------
| overwrite : Boolean, default False
|
| Example
| -------
| >>>
| >>> myblockmodel.calc_ixyz_fromxyz()
| >>> myblockmodel.calc_ixyz_fromxyz(overwrite=True)
| >>>
|
| calc_xyz_fromixyz(...)
| calc_xyz_fromixyz(bint overwrite=False)
|
| Calculates the XC, YC, ZC coordinates from IX, IY, IZ indices
| If XC, YC, ZC exist and overwrite=True the existing values
| will be overwritten.
|
|
| Parameters
| ----------
| overwrite : Boolean, default False
|
| Examples
| --------
| >>>
| >>> myblockmodel.calc_xyz_fromixyz()
| >>> myblockmodel.calc_xyz_fromixyz(overwrite=True)
| >>>
|
| create_3Dgrid(...)
| create_3Dgrid(bint overwrite=False)
|
| Creates a full block 3D grid/block model
|
|
| Parameters
| ----------
| overwrite : boolean
| overwrite flag, if true the entire grid will be
| overwritten
|
| Example
| -------
| >>>
| >>> mymodel.create_3Dgrid(overwrite=False)
| >>>
|
| create_IJK(...)
| create_IJK(bint overwrite=False)
|
| Creates a new block model table with IJK indices.
|
| Examples
| --------
| >>>
| >>> create_IJK(overwrite=True)
| >>>
|
| Note
| ----
| A new set of blocks will be created if there is not block
| defined. If there are blocks in the model and overwrite==True
| the blocks will be removed first.
|
| create_gridxy(...)
| create_gridxy(bint overwrite=False)
|
| Creates a full block 2D grid in defined in the XY direction.
|
| The grid definition is taken from the block model definition in
| the X, and Y directions.
|
|
| Parameters
| ----------
| overwrite : boolean
| overwrite flag, if true the entire grid will be
| overwritten
|
| Example
| -------
| >>>
| >>> mymodel.create_gridxy(overwrite=False)
| >>>
|
| create_gridxz(...)
| create_gridxz(bint overwrite=False)
|
| Creates a full block 2D grid in defined in the XZ direction.
|
| The grid definition is taken from the block model definition in
| the X, and Z directions.
|
|
| Parameters
| ----------
| overwrite : boolean
| overwrite flag, if true the entire grid will be
| overwritten
|
| Example
| -------
| >>>
| >>> mymodel.create_gridxz(overwrite=False)
| >>>
|
| create_gridyz(...)
| create_gridyz(bint overwrite=False)
|
| Creates a full block 2D grid in defined in the YZ direction.
|
| The grid definition is taken from the block model definition in
| the Y, and Z directions.
|
|
| Parameters
| ----------
| overwrite : boolean
| overwrite flag, if true the entire grid will be
| overwritten
|
| Example
| -------
| >>>
| >>> mymodel.create_gridyz(overwrite=False)
| >>>
|
| delete_blocks(...)
| delete_blocks()
|
| Deletes the block model table
|
| Examples
| --------
| >>>
| >>> myblockmodel.delete_blocks()
| >>> print myblockmodel.bmtable
| None
| >>>
|
| Note
| ----
| This functions makes Blockmodel.bmtable=None. The data will
| be preserved in any external instance of bmtable.
|
| fillwireframe(...)
| fillwireframe(object surface, float toll=0, bint overwrite=False)
|
| Creates a full block model given a VTK surface (polydata) using
| vtkPolyDataToImageStencil. The results consist of blocks with
| parameter ``in`` with value between 0 and 1.
|
| ``in = 0`` represent blocks completely outside the wireframe
| ``in = 1`` represent blocks completely inside the wireframe
| ``1 >in > 0`` represent blocks cut by the wireframe
|
| The parameter in is calculated as the sum of corners of the
| block, each corner is equal to 1 if is inside the block or equal
| to 0 if is outside.
|
| The parameter ``in`` can be used as percentage of inclusion in
| the block with an error proportional to the size of the block
| and the resolution/complexity of the wireframe.
|
| The method may fail if the wireframe is narrower than a block,
| e.j. larg blocks in narrow veins. In this case you can set
| a tolerance value 1>=toll>=0. The output may not be used as
| volume, only as selection of blocks within a wireframe
| (with ``in>0``).
|
|
| Parameters
| ----------
| surface : VTK polydata
| this may work with any 3D object..., no necessarily triangles
| tol : float, default 0.0
| number in interval [0. , 1.]
| you may use tol>0 if the blocks are large and the
| wireframe is narrower than the block size.
| overwrite : boolean
| overwrite flag, if true the entire block model will be
| overwritten
|
| Example
| -------
| >>>
| >>> mymodel.fillwireframe(surface, toll=0, overwrite=False)
| >>>
|
| Note
| ----
| The tolerance only works on x, y plane, not in z.
|
| This function works well with closed surfaces
|
| This function works well with large and complicated wireframes,
| for example, Leapfrog grade shells.
|
| point2block(...)
| point2block(np.ndarray [double, ndim=1] x, np.ndarray [double, ndim=1] y, np.ndarray [double, ndim=1] z,object prop, str prop_name, bint average = False, bint overwrite = False)
|
| Assigns data in a point to a block
|
| Parameters
| ----------
| x,y,z : Float 1D numpy arrays
| coordinates of the points. Array shapes may be equal
| prop : array like
| array with property shape[0] == x.shape[0]
| note that prop can be an array of array, for example, an
| array of indicators arrays
| prop_name: str
| name of the property at block model
| average : bool
| If there are more than one point at block
| average the points. If true the average will be created,
| this may fail if prop is not numeric. If False the last
| point in the block will be used to assign the value.
| overwrite : bool
| If False an existing property will not be overwritten
|
| Example
| -------
| >>>
| >>> point2block( x, y, z, prop, prop_name= 'Au', average = False, overwrite = False)
| >>>
|
| TODO: pass array of properties np.ndarray [double, ndim=2] prop and [prop_names]
|
| Note
| ----
| Points not laying in an existing block will be ignored.
|
| reblock(...)
| reblock(double dx, double dy, double dz)
|
| Reblocks a block model
|
| The model is reblocked using IJK calculated from XC, YC, ZC
| coordinates and averaging all points with same IJK.
|
| Note that this is like a point to block conversion and subblock
| overlapping are not taken into account.
|
| Parameters
| ----------
| dx,dy,dz: floats
| New block dimensions. Note that new dimentions may be
| greater or equal than old dimensions in all directions.
|
|
| Returns
| -------
| block: model instance with model reblocked
|
| err: list of variables excluded (ej no numeric)
|
| Warning
| -------
| The existing model will be overwritten
|
| set_block_size(...)
| set_block_size(float dx,float dy, float dz)
|
| Set block sizes
|
| Parameters
| ----------
| dx,dy,dz : float
| sizes of the blocks
|
| Examples
| --------
| >>>
| >>> myblockmodel.set_block_size(dx=10,dy=10,dz=5)
| >>>
|
| set_blocks(...)
| set_blocks(object bmtable)
|
| Assigns an external block model stored in a Pandas DataFrame
| table.
|
| One or many conditions apply
|
| - the block model has IJK field
| - or/and has IX, IY and IZ fields
| - or/and has XC, YC and ZC fields
|
|
| Parameters
| ----------
| bmtable : Pandas DataFrame object
| block model with special coordinate or index fields
|
|
| Examples
| --------
| >>>
| >>> bmtable=pd.DataFrame({'IJK':np.array([1,2,3,4], dtype=np.int32)})
| >>> myblockmodel.set_blocks(bmtable)
| >>>
|
| Note
| ----
| Make sure IJK, IX,IY and IZ have dtype int32. XC,YC and ZC may
| have dtype float32
|
| One of the above conditions is required. In case that two or
| more of the special fields mentioned above are defined the
| validity of the fields is not verified, for example the
| valid IJK representation of fields IX, IY and IZ is not tested.
|
| set_origin(...)
| set_origin(float xorg,float yorg, float zorg)
|
| Set the block model origin of coordinate. This is the lower
| left corner of the lower left block (not-the centroid).
|
| Parameters
| ----------
| xorg,yorg,zorg : float
| coordinates of the left lower corner of the first block
|
| Examples
| --------
| >>>
| >>> myblockmodel.set_origin(xorg=-5,yorg=-5, zorg=0)
| >>>
|
| set_rcl(...)
| set_rcl(int nx,int ny, int nz)
|
| Set number of blocks at row, call, level
|
| Two conditions apply
|
| - nx,ny and nz may be greater than zero
| - nx*ny*nz may be equal to the actual number of blocks
|
| Examples
| --------
| >>>
| >>> myblockmodel.set_rcl(nx=10,ny=10,nz=10)
| >>> myblockmodel.set_rcl(nx=10,ny=100,nz=1)
| >>>
|
| Note
| ----
| This function basically reorders the blocks along rows, columns
| and levels without changing the total number of blocks.
|
| ----------------------------------------------------------------------
| Data descriptors defined here:
|
| bmtable
|
| dx
|
| dy
|
| dz
|
| grid2DXY
|
| grid2DXZ
|
| grid2DYZ
|
| ndx
|
| ndy
|
| ndz
|
| nx
|
| ny
|
| nz
|
| xorg
|
| yorg
|
| zorg
|
| ----------------------------------------------------------------------
| Data and other attributes defined here:
|
| __new__ = <built-in method __new__ of type object>
| T.__new__(S, ...) -> a new object with type S, a subtype of T
|
| __pyx_vtable__ = <capsule object NULL>
FUNCTIONS
ijk2ind(...)
ijk2ind(np.ndarray [long, ndim=1] ijk, unsigned int nx, unsigned int ny, unsigned int nz)
Calculates the raw, column, level indices ``ix, iy, iz`` from IJK.
The IJK block index is a unique identifier of each block position.
This is equivalent to the position of a block in a `gslib` grid file.
All the points within a block will have the same IJK number.
From IJK you can calculate ix, iy and iz as::
iz = ijk / nx*ny
iy = (ijk-iz*nx*ny)/nx
ix = ijk-iz*nx*ny - iy*nx
Parameters
----------
ijk : 1D array of integers
arbitrary raw, level and column indices
nx,ny,nz : integers
number of blocks per row, column, level
Returns
-------
ix,iy,iz : 1D array of integers
The raw, column and level indices
Example
-------
>>>
>>> ijk= np.array([0, 1, 2, 3, 4, 5, 6, 7])
>>> ix,iy,iz = ijk2ind(ijk,2,2,2)
>>> print ix
>>> print iy
>>> print iz
[0 1 0 1 0 1 0 1]
[0 0 1 1 0 0 1 1]
[0 0 0 0 1 1 1 1]
>>>
See Also
--------
x2ix, ijk2ind
Note
----
The indices ix, iy and iz start at zero
ind2ijk(...)
ind2ijk(np.ndarray [long, ndim=1] ix, np.ndarray [long, ndim=1] iy, np.ndarray [long, ndim=1] iz, unsigned int nx, unsigned int ny, unsigned int nz)
Calculates the IJK block index
The IJK block index is a unique identifier of each block position.
This is equivalent to the position of a block in a `gslib` grid file.
All the points within a block will have the same IJK number.
IJK is calculated as follow:
``ijk = iz*nx*ny+ iy*nx + ix ``
Parameters
----------
ix,iy,iz : 1D array of integers
arbitrary raw, level and column indices
nx,ny,nz : integers
number of blocks per row, column, level
Returns
-------
ijk : 1D array of integers
Unique identifier of block location
Example
--------
>>>
>>> # a 2D grid with 2x2x1 cells
>>> ix=np.array([0,1,0,1])
>>> iy=np.array([1,1,0,0])
>>> iz=np.array([0,0,0,0])
>>> ind2ijk(ix,iy,iz,2,2,1)
array([2, 3, 0, 1])
>>>
See Also
--------
x2ix, ijk2ind
Note
-----
The index IJK start at zero (first block) and ends at nx*ny*nz
ix2x(...)
ix2x(np.ndarray [long, ndim=1] ix, double xorg, double dx)
Calculates the block coordinate (x, y or z)
Return array of coordinates x calculated from indices ix.
Parameters
----------
ix : 1D array of positive integers
arbitrary index ix (or iy, or iz)
xorg : float
origin of coordinate (lower left corner of the block)
dx : float
size of the block
Returns
-------
x : 1D array of floats
Coordinates x corresponding to index ix
See Also
--------
x2ix
Notes
-----
The index starts at zero (first block)
If there is a point i with coordinate x[i] < xorg and error will be
raised
x2ix(...)
x2ix(np.ndarray [double, ndim=1] x, double xorg, double dx)
Calculates the block index (ix, iy or iz)
Return array of ix indices representing the raw, column or level
in which a point with coordinate x is located. To get the three
indexes ix, iy, iz you may run this functions three times, changing
coordinates to x, y and z, respectively.
You can use this function to find the block which contain a points
with arbitrary coordinates.
Parameters
----------
x : 1D array of floats
arbitrary coordinates x (or y, or z)
xorg : float
origin of coordinate (lower left corner of the block)
dx : float
size of the block
Returns
-------
ix : 1D array of integers
Index of the blocks where points with coordinates x are located
Example
-------
>>>
>>> a=np.array([2.1,7.8,15.9,20,30,31])
>>> x2ix(a,xorg=2., dx=1)
array([ 0, 5, 13, 18, 28, 29])
>>>
See Also
--------
ind2ijk, ijk2ind
Note
----
The index starts at zero (first block)
If there is a point i with coordinate x[i] < xorg and error will be
raised
DATA
__test__ = {u'Blockmodel.block2point (line 946)': u"block2point(np.nda...
In [3]:
# Create an empty block model
mymodel=pygslib.blockmodel.Blockmodel(nx=5,ny=5,nz=5,xorg=-6,yorg=-6,zorg=-6,dx=3,dy=3,dz=3)
In [4]:
# there is no data, it is empty
print mymodel.bmtable
None
In [5]:
# Create a model from external DataFrame.
# In this case IJK or IX, IY, IZ or XC, YC, ZC are required
bmtable=pd.DataFrame({'IJK':np.array([1,2,3,4], dtype=np.int32)})
mymodel.set_blocks(bmtable)
print mymodel.bmtable
IJK
0 1
1 2
2 3
3 4
In [6]:
bmtable.dtypes
Out[6]:
IJK int32
dtype: object
In [7]:
# Calculate row, column and level indices from IJK
mymodel.calc_ixyz_fromijk(overwrite=False)
mymodel.calc_ixyz_fromijk(overwrite=True) # if IX, IY, IZ exist we need to overwrite
print mymodel.bmtable
IJK IX IY IZ
0 1 1 0 0
1 2 2 0 0
2 3 3 0 0
3 4 4 0 0
In [8]:
# Removing data from block model
mymodel.delete_blocks()
print mymodel.bmtable
None
In [9]:
# Creating a full block model with IJK
mymodel.create_IJK(overwrite=True)
print mymodel.bmtable.tail()
IJK
120 120
121 121
122 122
123 123
124 124
In [10]:
# Calculate row, column and level indices from IJK
mymodel.calc_ixyz_fromijk(overwrite=True)
print mymodel.bmtable.tail()
IJK IX IY IZ
120 120 0 4 4
121 121 1 4 4
122 122 2 4 4
123 123 3 4 4
124 124 4 4 4
In [11]:
# Calculate coordinates from indices
mymodel.calc_xyz_fromixyz(overwrite=True)
print mymodel.bmtable.tail()
IJK IX IY IZ XC YC ZC
120 120 0 4 4 -4.5 7.5 7.5
121 121 1 4 4 -1.5 7.5 7.5
122 122 2 4 4 1.5 7.5 7.5
123 123 3 4 4 4.5 7.5 7.5
124 124 4 4 4 7.5 7.5 7.5
In [12]:
mymodel.blocks2vtkRectilinearGrid('blocks')
The results can be ploted in an external viewer, for example mayavi or paraview:
In [ ]:
mymodel.
Content source: opengeostat/pygslib
Similar notebooks: