grail.ag Module

Version: 16.2

Table of Contents

Overview

A collection of routines for performing analysis on spatial and geometric data. The term "ag" comes from "Analytical Geometry".

Note

This module will not provide the full suite of standard matrix functions. Matrix work is done much, much more efficiently with the Numeric [1] library that is distributed alongside MineSight® Grail.

Standard Definitions

Some standards for working with this module,

Co-ordinates:
All co-ordinates are assumed to be in world co-ordinates, such that Z is up, X is to the left and Y is away from you.
Vectors:
A list of co-ordinates, defined as: [i, j, k]. Note that a vector and point follow the same simple definition.
Points:
A list of co-ordinates, defined as: [x, y, z].
Point Lists:
A list of points, defined as: [[x0, y0, z0], ..., [xN, yN, zN]].
Face Lists:

All face lists are described such that the first index indicates how many faces are in the face. For example,

[3, 1, 2, 3, 3, 4, 5, 6]

States that there are 3 vertexes in the first face, and those vertexes are 1, 2, 3. Then the next face contains 3 vertexes, which are 4, 5, 6. This may seem a bit confusing, but using this structure buys us some definite efficiencies.

Note that all faces in MineSight® products have 3 vertexes.

Plane Equation:

The plane equation is defined as,

a + b + c + d = 0

Whenever you are asked for a, b, c, and d, these match to the above plane equation.

Exceptions Handling

In general most functions in this module will either raise an ArithmeticError or one of its sub-classes, FloatingPointError [2].

This exceptions are defined by the Python Library itself.

Functions

The functions in this module can be divided into four distinct groups. These are functions that work with vectors, functions that are used in unit conversions, functions used in spatial/geometric computations, and functions used to compute equality.

Conversion Functions

cm2inch(centimeters)

Converts centimeters into inches.

feet2meter(feet)

Converts feet into meters.

inch2cm(inches)

Converts inches into centimeters.

meter2feet(meters)

Converts meters into feet.

round_sig(x, digits)

Rounds the value of x to the number of digits given. The standard is the value given for ag.STD_RND_PRECISION.

A value of digits that is less than 1 will generate a ValueError.

todegrees(radians)

Converts radians to degrees.

toradians(degrees)

Converts degrees to radians.

Equality Function

isequal(a, b[, epsilon])

Determines if two real numbers are approximately equal.

Since mathematical reals have infinite precision, and since a machine can not represent infinite precision, there is no precise equality on a machine.

Instead this function looks at the absolute value of the difference and compares it against epsilon, if the difference is less than epsilon then the numbers are considered equal.

The comparison can be written as,

isequal = (| a - b | < epsilon)

Note that specifying your epsilon value leads to an approximate 30% decrease in efficiency.

The default value for epsilon is ag.EPSILON.

Returns:
1 if the absolute difference is less than epsilon; 0 otherwise.

Geometric/Spatial Functions

shiftPlane(a, b, c, d, vector)

Returns a plane based on the argument plane (a, b, c, d) and the vector.

Arguments:
a, b, c, d : floats
Plane (see Standard Definitions).
vector : list
Based on standard point definition
Returns:
The Plane shifted by supplied vector : list of plane coefficients Invalid planes, such as (0.0, 0.0, 0.0, 0.0) will simply be returned without being shifted.

Azimuth and Dip

get2dazimuth(x1, y1, x2, y2)

Returns azimuth from point (x1, y1) to point (x2,y2) on a plane.

Arguments:
x1, y1, x2, y2 : floats
Points (see Standard Definitions).
Returns:
The azimuth of the 2D points.
getdip(p1, p2)

Returns the dip from one point to another.

Arguments:
p1, p2 : lists
Points (see Standard Definitions).
Returns:
Dip of a vector from p1 to p2: float

Points

acuteAngle2D(p1, p2, p3)

Returns the 2D angle between line segments p1p2 and p2p3 (where p2 is the elbow). This will only consider the X,Y components of the points (i.e indices 0 and 1).

Arguments:
p1, p2, p3 : list
Points (see Standard Definitions).
Returns:
The angle in degrees
findEndPoint(point, dip, azimuth, length)

Finds the endpoint of a line given a starting point, 3d length, dip, and azimuth. New in version 4.70.

Arguments:
point: list
Point (see Standard Definitions).
dip: float
dip in degrees
azimuth: float
azimiuth in degrees
length: float
3d length of line
Returns:
The end point (see Standard Definitions).
getBlockCentroid(pcfPath, lev, row, col)

Finds the centroid of a block model cell. Support for regular models and horizontally-rotated models. New in version 5.10.

Arguments:
pcfPath: string
Full path to the Project Control File.
lev: integer
Level index.
row: integer
Row index.
col: integer
Column index.
Returns:
The centroid point (see Standard Definitions).
intersect2D(p1, p2, p3, p4)

Determines if line segments p1p2 and p3p4 intersect. This will only consider the X,Y components of the points (i.e indices 0 and 1).

Arguments:
p1, p2, p3, p4 : list
Points (see Standard Definitions).

Returns 1 if the segments intersect, 0 otherwise.

ispointcoincident(p1, p2 [, epsilon])

Indicates if two points are nearly the same (within a given epsilon).

Arguments:

p1, p2 : list
Points (see Standard Definitions).
epsilon : float
The tolerance that decides the level of equality (default = ag.EPSILON).
Returns:
1 if the two points are coincident (nearly the same) and 0 otherwise.
nearestpointonline(point, lp, lv)

Determines the nearest point on the line to the given point.

Given a 3-D unbounded line, this computes the nearest point on that line to the point given.

See Standard Definitions for more information about how a point is defined.

Arguments:
point : list
The search point.
lp : list
One point on the 3-D line.
lv : list
Another point on the 3-D line.
Returns:
A point that lies on the line defined by lp, lv.
point_inside_pointlist2d(point, pointlist)

Determines if the point is within the given pointlist. This will only consider the X,Y components of point and pointlist (i.e indices 0 and 1).

Note

If you supply a pointlist that is not a closed polygon (i.e. the end point is not the same as the start point), then this algorithm will effectively "close" the pointlist.

This may not be the desired behavior, and you can check your if your pointlist is a polygon using the polygon_check function.

pointdistance(p1, p2)

Returns the distance between two points.

Arguments:
p1, p2 : list
Points (see Standard Definitions).
Returns:
The distance between the two points.

pointdistancefromplane(p, a, b, c, d)

Returns the distance from a given plane.

Arguments:
p : list
A point (see Standard Definitions).
a, b, c, d : floats
Coefficients for the plane equation (see Standard Definitions).
Returns:
The distance from the plane.

Point Lists

ispointlistaline(pointlist)

Determine if the point list defines a line.

Line is defined by taking the vector between two points and determining if there is any variation from one segment to the next. Therefore the following two point lists would be considered a line,

points = [
    
[1.0, 1.0, 1.0],
    
[2.0, 2.0, 2.0],
    
[3.0, 3.0, 3.0],
    
]

points = [
    
[2.0, 2.0, 2.0],
    
[1.0, 1.0, 1.0],
    
[3.0, 3.0, 3.0],
    
]

But this would not be a valid line,

points = [
    
[0.0, 0.0, 12121.0],
    
[-1., 0.0, 0.0],
    
[0.0, 1.0, 0.0],
    
]
Argument:
pointlist : list
A list of points (see Standard Definitions).

Returns 1 if the points are define a line, 0 otherwise.

ispointlistonplane(pointlist, a, b, c, d [,epsilon])

Indicates if a group of points exist on the specified plane.

Arguments:
pointlist : list
A list of points (see Standard Definitions).
a, b, c, d : floats
Coefficients for the plane equation (see Standard Definitions).
epsilon : float
The tolerance that decides the level of equality (default = ag.EPSILON).
Returns:
1 if all the points exist on the plane within the given epsilon; otherwise return 0.
pointlist_equal(pointlist1, pointlist2[, epsilon])
Compares pointlist1 to pointlist2 and returns true (1) if they both have the same number of points, and each point is coincident to within the given epsilon. The default value for epsilon is ag.EPSILON.
pointlist_shift(pointlist, shiftVector)

Shifts the contents of the pointlist by the given shiftVector. Useful for shifting points quickly based on some given offset.

Argument:
pointlist : list
A list of points (see Standard Definitions).
Returns:
A new shifted pointlist.
pointlist2dlength(pointlist)

Returns the length of a point list in 2-D.

This just sums along from one point to the next only considering the XY dimensions.

Argument:
pointlist : list
A list of points (see Standard Definitions).
Returns:
The length of the lines.
pointlistarea(pointlist)

Determines the area of a list of points.

To get the correct area all points should lie on the same plane. Furthermore, it is implicit in the calculation of the area that the points in the point list are supposed to connect completely (for example, a polygon).

Argument:
pointlist : list
A list of points (see Standard Definitions).
Returns:
The area of the pointlist. ArithmeticError [2] is raised if the area can not be computed (i.e. the pointlist is not something that defines an area, or there is a duplicate end point).
pointlistboundingvolume(pointlist)

Returns the minimum and maximum points in the point list.

The minimum and maximum points describe the bounding volume for the list of all points. They should describe the extreme corners of a box that is the bounding volume.

Notice that this calculation can apply to a polyline or to just a cloud of points.

Argument:
pointlist : list
A list of points (see Standard Definitions).
Returns
A tuple containing two points, and defined as, ([min_x, min_y, min_z], [max_x, max_y, max_z]).
pointlistcenter(pointlist)

Computes the center of mass for a point list.

If there is only one point in the list, then that point is returned.

If all the points on the list are co-incident, then the first point is returned.

Note

If your data is polygon (end-point repeated) this returns the center of the polygon.

Usage:

Consider the simple case where we have a square polygon defined by the following points,

from grail import ag
points = [[2.0,  0.0,  0.0],
          
[2.0,  2.0,  0.0],
          
[0.0,  2.0,  0.0],
          
[0.0,  0.0,  0.0],
          
[2.0,  0.0,  0.0]]
center = ag.pointlistcenter(points)

At this point center will be equal to [1., 1., 0.].

The above call could have also been made as,

x, y, z = ag.pointlistcenter(points)

With this call, you would have x=1.0, y=1.0, and z=0.0.

Argument:
pointlist : list
A list of points (see Standard Definitions).
Returns:
The center point (see Standard Definitions).
pointlistdirection(pointlist)

Returns the direction of the points.

The points are read in sequence to determine if they are leading in a clockwise or counter-clockwise direction.

Argument:
pointlist : list
A list of points (see Standard Definitions).
Returns:
Either ag.COUNTER_CLOCKWISE or ag.CLOCKWISE. Both values are defined in this module (see Constants).
pointlistlength(pointlist)

Returns the length of segments along the point list.

Traverses the points list and returns the total length. Just like walking along the line with a pedometer in 3D space.

Argument:
pointlist : list
A list of points (see Standard Definitions).
Returns:
The length of the lines.
pointlistplane(pointlist)

Returns the plane equation for the pointlist.

Arguments:
pointlist : list
A list of points (see Standard Definitions).
Returns:
a, b, c, d : floats
Coefficients for the plane equation (see Standard Definitions).
rotateX(pointlist, angle)

Rotate each point in pointlist about the X-axis.

Arguments:
pointlist : list
A list of points (see Standard Definitions).
angle : float
Angle of rotation in degrees (positive=CW, negative=CCW)
Returns:
The rotated pointlist
rotateY(pointlist, angle)

Rotate each point in pointlist about the Y-axis.

Arguments:
pointlist : list
A list of points (see Standard Definitions).
angle : float
Angle of rotation in degrees (positive=CCW, negative=CW)
Returns:
The rotated pointlist
rotateZ(pointlist, angle)

Rotate each point in pointlist about the Z-axis.

Arguments:
pointlist : list
A list of points (see Standard Definitions).
angle : float
Angle of rotation in degrees (positive=CW, negative=CCW)
Returns:
The rotated pointlist

Planes

lineIntersectsPlane(p1, p2, a, b, c, d)

Determines if the line segment represented by endpoints p1 and p2 intersects the plane. New in version 4.70.

Arguments:
p1, p2 : list
Points (see Standard Definitions).
a, b, c, d : floats
Coefficients for the plane equation (see Standard Definitions).
Returns:
The point of intersection (if a single intersection point is found); 1 if line is coincident to plane (infinite solutions); otherwise returns 0.
nearestpointonplane(point, a, b, c, d)

Determines the nearest point on the plane to the given point.

Given a 3-D plane equation co-efficients and a point, this returns the nearest point on the plane. Can also be thought of as projecting a point directly onto that plane.

Arguments:
pointlist : list
A list of points (see Standard Definitions).
a, b, c, d : floats
Coefficients for the plane equation (see Standard Definitions).
Returns:
A point that is nearest the given point, and lies on the defined plane.
planenormal(a, b, c, d)

Computes the normal to the given plane.

Arguments:
a, b, c, d : floats
Coefficients for the plane equation (see Standard Definitions).
Returns:
A vector that defines the normal to the inputed plane.
plane_reorient_normal(a, b, c, d)

Reorients the plane defined by a, b, c, and d into the standard used by MS3D. This is useful for when you have the following plane equation:

0.00000     1.00000     0.00000    1895.00000

and:

0.00000    -1.00000     0.00000   -1895.00000

both are the same, but the latter one is pointing downwad. The big problem comes from the fact that MS3D will report two different planes one at 1895.00 and one at -1895. By reorienting your planes in a consistent manner, you guarantee that they are stored correctly.

An example of calling this function is,

        nrm_pln = (0.00000, 1.00000, 0.00000, 1895.00000,)
        
dwn_pln = (0.00000, -1.00000, 0.00000, -1895.00000,)

        
reorient_pln = ag.plane_reorient_normal(*dwn_pln)

after calling plane_reorient_normal, the reorient_pln will be equal to the nrm_pln. Notice the trick using the * to pass in a tuple with 4 elements as a, b, c, and d.

If you supply an invalid plane equation, such as 0, 0, 0, 0 or 0, 0, 0, N, where N is any non-zero number, this equation will not generate an error. Instead it will merrily attempt to reorient the plane.

Polygons

polygon_check(pointlist [,epsilon])

Returns True if the pointlist is a polygon (i.e. it has an end point that is equivalent to it's start point).

The epsilon is used to help determine equality similar to ispointcoincident. It has a default value defined as EPSILON.

isConvex(pointlist)

Determine if the polygon is convex or concave.

Arguments:
pointlist : list
A list of points (see Standard Definitions).

Returns 1 if the polygon is convex, 0 if the polygon is concave.

isSelfIntersecting(pointlist)

Detects whether a polygon intersects itself.

Arguments:
pointlist : list
A list of points (see Standard Definitions).

Returns 1 if the polygon is self-intersecting, 0 otherwise.

boundingPolygon(point_list)

Returns a polygon containing the vertices of a minimal bounding hull of the point list.

Arguments:
point_list : list
Points (see Standard Definitions).
Returns:
A polygon containing all points in the point_list.

Smoothing

smoothaverage(pointlist, averagecount, pin)

Returns the new smoothed pointlist.

Creates a new smoothed pointlist by computing a running average of the points for all the nodes within the pointlist.

The averagecount indicates how many nodes to average across.

Arguments:
pointlist : list
A list of points (see Standard Definitions).
averagecount : integer

An integer equal to the number of nodes to average when creating the new smoothed pointlist.

The minimum value for averagecount is 2.

pin : integer

An integer indicating whether to preserve end points.

SMOOTH_PIN : preserve end points

SMOOTH_NO_PIN : do not preserve end points

Pinning applies to polylines only, which are defined as a pointlist where the first point is not equal to the last point.

Returns:
Smoothed pointlist or None if an error encountered.

Usage:

This function applies to polyline and polygon pointlists.

For this function to work correctly you will require at least three points to perform an average across. If you have less than three points this function will return the original pointlist.

For speed, you may want to screen the pointlist to see if it actually can be smoothed before calling this function. If the pointlist can not be smoothed, the original pointlist is returned. Examples of pointlists that can not be smoothed are: 1. a polyline with two points; 2. a polygon or polyline having a total number of points less than or equal to averagecount.

smoothbezier(pointlist, addnodes, smoothfactor)

Returns the new smoothed pointlist. Function also known as "Node Preserving Spline."

Creates a new smoothed pointlist by applying the Bezier algorithm to the polyline/polygon for all nodes in the pointlist. The Bezier smoothing algorithm will intersect nodes while it smooths the curve.

The addnodes specifies how many nodes to add during the smoothing process, and the smoothfactor indicates how smooth to make the curve. The larger the smoothfactor the greater the curve will deviate from the original.

Arguments:
pointlist : list
A list of points (see Standard Definitions).
addnodes : integer
An integer equal to the number of nodes to add between points in the pointlist during the smoothing process.
smoothfactor : float
An integer indicating how smooth to make the curve. 'smoothfactor' must be greater than zero.
Returns:
Smoothed pointlist or None if an error encountered.

Usage:

This function applies to polyline and polygon pointlists.

For this function to work correctly you will require at least three points to perform the smooth. If you have less than three points this function will return the original pointlist.

For speed, you may want to screen the pointlist to see if it actually can be smoothed before calling this function.

smoothspline(pointlist, addnodes, pin)

Returns the new smoothed pointlist.

Creates a new smoothed pointlist by applying the Spline smoothing algorithm to the polyline/polygon for all nodes within the pointlist. The Spline smoothing algorithm does not respect nodes while it smooths the curve.

The addnodes specifies how many nodes to add during the smoothing process.

Arguments:
pointlist : list
A list of points (see Standard Definitions).
addnodes : integer
An integer equal to the number of nodes to add between points in the pointlist during the smoothing process.
pin : integer

An integer indicating whether to preserve end points.

SMOOTH_PIN : preserve end points

SMOOTH_NO_PIN : do not preserve end points

Pinning applies to polylines only, which are defined as a pointlist where the first point is not equal to the last point.

Returns:
Smoothed pointlist or None if an error encountered.

Usage:

This function applies to polyline and polygon pointlists.

For this function to work correctly you will require at least three points to perform the smooth. If you have less than three points this function will return the original pointlist.

For speed, you may want to screen the pointlist to see if it actually can be smoothed before calling this function.

surfacearea(pointlist, facelist)

Computes the area for a surface.

Arguments:
pointlist : list
A list of points (see Standard Definitions).
:a:facelist : list
A list of faces using the 'hoops' convention (see Standard Definitions).
Returns:
The area of a surface.

Vector Functions

angle2d(x, y)

Returns the angle between point at x,y and the origin.

Positive angle computed in a counter-clockwise fashion. When x is 0 and y is 0, the return angle will be 0.

Arguments:

x : float
A point along the x-axis.
y : float
A point along the y-axis.
Returns:
The angle in radians between the origin and the x,y point.
crossproduct(v1, v2)

Computes the cross product between two vectors.

Arguments:

v1, v2 : list
Vectors (see Standard Definitions).
Returns:
Another vector which is the crossproduct between v1 crossed with v2.
dotproduct(v1, v2)

Computes the dot product between two vectors.

Arguments:

v1, v2 : list
Vectors (see Standard Definitions).
Returns:

The dot product between v1 and v2, in other words,

v1.i*v2.i+v1.j*v2.j+v1.k*v2.k
vectorlength(vector)

Returns the length of a vector.

Argument:
vector : list
A vector (see Standard Definitions).
Returns:
The length of the vector.
vectorscale(vector, scale)

Scales a vector.

Arguments:
vector : list
A vector (see Standard Definitions).
scale : float
A scaling factor to apply to each element of a vector.
Returns:

A vector [i, j, k] that is scaled by the scaling factor.

Equivalent to,

[el*scale for el in vector]

Matrix Functions

changebasisfromunitvectors(x, y, z)

Returns change basis matrix as a 4x4 list of floats.

Arguments:

x : list
Vector (see Standard Definitions).
y : list
Vector (see Standard Definitions).
z : list
Vector (see Standard Definitions).

example: [1,0,0], [0,1,0], [0,0,1]

Returns:
The change basis matrix as a 4x4 list of floats or None if error.
Usage:

This will create a change basis matrix going from cartesian coordinates to the given cartesian coordinates in the x,y,z vectors out in space.

To do the opposite of going from the x,y,z vectors out in space to 0,0,0 cartesian then transpose the matrix.

The x, y, z points in represent the coordinate system out in space. These should be at right angles to each other. In other words, the dot product of each vector with any of the other vectors will be 0.

matrixfromposition(position, origin)

Returns modeling matrix for the desired point.

Arguments:

position : a list of three vectors representing the point, up vector, right vector

example: [[59661.46,58083.77,0.00],[59661.46,58083.77,1.00],[59662.46,58083.77,0.00]]

origin : a point representing the origin of the geometry object

example: [55300.0, 52300.0, -400.0]
Returns:
A modeling matrix as a 4x4 list of floats or None if error.
positionfrommatrix(matrix, origin)

Returns position from the modeling matrix and origin of the geometry object.

Arguments:

matrix : a 4x4 list of floats representing the modeling matrix

example: [[1.0, 0.0, 0.0, 4361.],[0.0, 1.0, 0.0, 5783.],[0.0, 0.0, -1.0, 100.0],[0.0, 0.0, 0.0, 1.0]]

origin : a point representing the origin of the geometry object

example: (55300.0, 52300.0, -400.0)
Returns:
The position as a list of three vectors representing the point, up vector, right vector
transposematrix(matrix)

Returns transposed 4x4 matrix of floats.

Arguments:

matrix : List of 4 lists (or tuples) containing four floats each.
Returns:
The transposed 4x4 matrix (list of 4 lists of 4 floats), or None if error.

Constants

CLOCKWISE COUNTER_CLOCKWISE

Constants indicate direction of polyline/polygon.

EPSILON

Default constant generally used to determine closeness between two floating point numbers. The value is 0.001.

MAX_REAL

The maximum size of a real number on the operating system.

MAX_REAL_DIGITS

The maximum number of digits in a real number on the operating system.

SMOOTH_PIN SMOOTH_NO_PIN

Constants indicate whether to preserve original nodes of polyline/polygon.

STD_RND_PRECISION

Number of digits to use when doing a signficant figure reduction.

References

[1]"Homepage". Numerical Python. 11 May 2004.
[2](1, 2) "Built-in Exceptions" Python Library Reference. 19 December 2001. 31 May 2004. <http://www.python.org/doc/2.2/lib/module-exceptions.html>