JuliaEarth/GeoStats.jl

using GeoStats

# random fault as a Segment
function randomfault()
  A = rand(Point2)
  B = rand(Point2)
  Segment(A, B)
end

# midpoint of a fault
function midpoint(fault)
  A = coordinates(fault(0))
  B = coordinates(fault(1))
  Point((A + B) / 2)
end

# normal direction to a fault
function normal(fault)
  A = fault(0)
  B = fault(1)
  d = B - A
  # 90 degree rotation
  Vec(d[2], -d[1])
end

# simulate 10 random faults
faults  = [randomfault() for i in 1:10]

# define bisectors (partition methods)
points    = [midpoint(fault) for fault in faults]
normals   = [normal(fault) for fault in faults]
pairs     = zip(normals, points)
bisectors = [BisectPointPartition(normal, point) for (normal, point) in pairs]

# the product of bisectors creates compartments
compartments = reduce(*, bisectors)

# random data for Kriging
obs = georef((Z=rand(100),), rand(2,100))

# partition data to use in different Kriging calls
partition(obs, compartments)