Discontinuity

A discontinuity is a plane which is typically used to mark a change in the physical or chemical properties of a rock mass. They are represented in the Maptek Python SDK as a set of coplanar points which are turned into a surface using facets similar to a surface. Unlike a surface, discontinuities do not support any kind of point or facet property. The example below demonstrates how to use a set of coplanar points to define a new discontinuity.

from mapteksdk.project import Project
from mapteksdk.data import Discontinuity
 
points = [[0, 0, -0.5], [1, 0, 0], [0, 1, -2], [1, 1, -1.5]]
 
project = Project()
 
with project.new("geotechnical/-x+3y+2z+1=0", Discontinuity) as plane:
  plane.planar_points = points

The plane of the discontinuity is defined using the first three points. Any points beyond the first are projected to be on the plane defined by the first three points. This is not necessary for the discontinuity created in the above example because all four points used lie on the plane -x + 3y + 2z + 1 = 0. The discontinuity is shown below.

Alternatively, you can create a discontinuity by defining the dip, dip direction and location of the discontinuity where:

  • Dip is the angle the discontinuity is rotated by in the dip direction.

  • Dip direction is the angle about the z axis which the discontinuity is rotated in.

  • Location is the coordinate of the centre of the discontinuity.

When you use this method to create a discontinuity, the planar points are automatically generated. However, because no points were provided when creating the discontinuity, querying the length or area of the discontinuity results in NaN (“Not a Number”). The example below shows creating a discontinuity in this way.

import math
 
from mapteksdk.project import Project
from mapteksdk.data import Discontinuity
 
project = Project()
 
with project.new("geotechnical/angled_discontinuity", Discontinuity) as plane:
  plane.dip_direction = math.radians(32.5)
  plane.dip = math.radians(16.35)
  plane.location = [-12, 42.8, 32.7]

The image below shows the discontinuity created by the above example. Note that the annotations match with the angles in degrees used to create the discontinuity.

The most robust method for creating a discontinuity requires specifying both the points and the dip, dip direction and location. If you use this method to create a discontinuity, the points are projected onto the plane defined by the dip and dip direction and translated to be centred around the specified location. This allows for the shape of the discontinuity and its properties to be specified independently.

Here is a an example of how to specify both the points and the dip, dip direction and location using code.

import math
 
from mapteksdk.project import Project
from mapteksdk.data import Discontinuity
 
project = Project()
 
# These are the points of a unit hexagon.
sqrt3_2 = math.sqrt(3) / 2
points = [[0, 0, 0], [1, 0, 0], [0.5, sqrt3_2, 0], [-0.5, sqrt3_2, 0],
          [-1, 0, 0], [-0.5, -sqrt3_2, 0], [0.5, -sqrt3_2, 0]]
 
with project.new("geotechnical/both_discontinuity", Discontinuity) as plane:
  plane.planar_points = points
  plane.dip_direction = math.radians(32.5)
  plane.dip = math.radians(16.35)
  plane.location = [-12, 42.8, 32.7]

If you run the above code, you should produce a discontinuity which looks like the image below.

Discontinuities possess other properties aside from those shown in the previous examples. The following example queries many of these properties for the discontinuities created in the previous examples. Make sure to run all previous examples before running the example below.

import math
 
from mapteksdk.project import Project
 
project = Project()
 
discontinuities = ("-x+3y+2z+1=0", "angled_discontinuity", "both_discontinuity")
 
for discontinuity in discontinuities:
  path = f"geotechnical/{discontinuity}"
  print(f"Properties for: {discontinuity}")
  with project.read(path) as plane:
    print("Dip:", math.degrees(plane.dip))
    print("Dip direction:", math.degrees(plane.dip_direction))
    print("Strike:", math.degrees(plane.strike))
    print("Plunge", math.degrees(plane.plunge))
    print("Trend", math.degrees(plane.trend))
    print("Location:", plane.location)
    print("Area:", plane.area)
    print("Length:", plane.length)
    print("Colour", plane.planar_colour)
  print("-" * 25)