## Where communities thrive

• Join over 1.5M+ people
• Join over 100K+ communities
• Free without limits
##### Activity
yizhengw
@yizhengw
Thanks!!
Júlio Hoffimann
@juliohm
GeoStats.jl v0.26 has been released. New algorithms include Geostatistical Hierarchical Clustering (GHC), Block Kriging and simulation equivalents.
Ruslan Vorobev
@ruslan_vorobev:matrix.org
[m]
Hi, everyone. I am looking for a 2D kriging estimation that taking into account core faults. Faults is defined as polylines. In theory it not takes neighbors located across (on the other side) of the fault. Is there way to do something like this in GeoStats? I didn't find anything similar in the standard examples, maybe i have missed something. Can you point me at some info or example how to do this? I have no experience on GeoStats.
Júlio Hoffimann
@juliohm
Hi @ruslan_vorobev:matrix.org , you won't find a general geostats package out there with this functionality out of the box because it is too specific to your application. The good news is that it should be pretty easy to code something like that in GeoStats.jl because of its sophisticated set of geometric operations. Suppose you have a collection of faults represented as line segments, here is a code you can use to partition your data into comparments:
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)
You can partition both the data and the domain(obs) if necessary. The result of the partition is a collection of views as explained in the docs and you can use these views to defined various subproblems to solve with Kriging for example. You can also get the indices(p) of the partition to construct views of obs and domain(obs that are consistent. Let me know if you have questions.
Gitter doesn't support Julia syntax unfortunately. If you want to discuss with more code, let's move the conversation to our Zulip channel: https://julialang.zulipchat.com/#narrow/stream/276201-geostats.2Ejl
Ruslan Vorobev
@ruslan_vorobev:matrix.org
[m]
Júlio Hoffimann
@juliohm
Hi @ruslan_vorobev:matrix.org you should be able to cook a solution using our low-level Kriging estimator objects instead of the Kriging solver. Suppose you have a set of locations with data s1, s2, ..., sk as Point and a new location where you want to perform estimation s0. You can create segments si = Segment(si, s0) and check if they intersect with any fault segment using qi = si ∩ fault. You can then check if isnothing(qi) and include the point in the list of valid points for that location s0. After you have a list of points of interest, you can create the kriging estimator object just for the subset: https://juliaearth.github.io/GeoStats.jl/stable/kriging.html
If you are working in 2D I think you are all set. For 3D faults represented as planes, we still need to implement intersection algorithms between segments and planes in Meshes.jl. Contributions are very welcome.
Ruslan Vorobev
@ruslan_vorobev:matrix.org
[m]
Can I intject something with fault intersection check via neighborhood parameter, or I must write cusom Kriging Estimator?
Júlio Hoffimann
@juliohm
I think writing a custom function for your purposes will be more flexible. Although faults don't fit very well in the concept of neighborhoods, maybe we could think of a new neighoborhood type that takes geometric objects into account.
Ruslan Vorobev
@ruslan_vorobev:matrix.org
[m]
As I understand I need to create custom solver KrigingWithFaults, rewrite function preprocess(problem::EstimationProblem, solver::KrigingWithFaults), copy other from kriging solver. And use it in solve(problem, solver). Thank you for help.
Júlio Hoffimann
@juliohm
You don't need to do all that just to solve your problem. I'd write a simple function to check intersection and solve the problem with the low level estimators. The Kriging solver is only used if you want to write a custom neighborhood
If you van define a neighbohorhood object that satisfies the interface of other neighborhoods then you can use the Kriging solver
The estimators that I was talking about use simple fit and predict calls
Ruslan Vorobev
@ruslan_vorobev:matrix.org
[m]
I have not yet figured out how to call another neighbor search function by setting a different neighborhood structure. The problem is a certain intermediate structure (BallSearch) with algorithm parameters is always created. In general the plan is clear. If I will fail to replace the neighbors search algorithm through parameters, I can implement it through a custom solver or a hardcode of its function. Thanks alot for help in understanding how it works.
Júlio Hoffimann
@juliohm
You can specify neighborhood=Ellipsoid(...) for example: https://juliaearth.github.io/GeoStats.jl/stable/solvers/builtin.html#GeoEstimation.Kriging
These neighborhood objects are defined in Meshes.jl here: https://github.com/JuliaGeometry/Meshes.jl/tree/master/src/neighborhoods
So if you follow the same interface with your custom neighborhood that considers faults (collection of segments), you can pass it directly to Kriging.
It uses the neighborhood objects to perform the search.
Ruslan Vorobev
@ruslan_vorobev:matrix.org
[m]
Thanks for links. The problem is I pass neighborhood object (NormBall for example) to parameters, but neighbor search data object (BallSearch for example) is created inside function preprocess(problem::EstimationProblem, solver::Kriging) in https://github.com/JuliaEarth/GeoEstimation.jl/blob/master/src/kriging.jl
And highly likely I have to override preprocess function
Júlio Hoffimann
@juliohm
I see your point, the search method is a function of the neighborhood and just passing a new neighborhood won't work. We should probably think of ways to generalize these concepts with multiple-dispatch. Maybe after you have your own custom solution :+1:
@hameye
Hi everyone, I'm trying to use the ImageQuilting.
@juliohm I saw in your paper you used Flumy to generate your 3D example, did you just plug the output from Flumy as a 3D array with value of facies ?
Júlio Hoffimann
@juliohm
Hi @hameye , yes, I plugged the output of Flumy directly with NaN values for inactive cells. Back in the days ImageQuilting.jl had GPU support for large 3D cases, but Julia changed dramatically the way GPU code is implemented and I didn't have time yet to reintroduce the feature in recent versions of the language. Other recent updates in the LightGraphs.jl stack and in the language may have introduced issues with the cut algorithm. If you have time to double check the solver, that would be a great contribution. For the GPU version to work, we need to update this filtering function: https://github.com/JuliaEarth/ImageQuilting.jl/blob/master/src/imfilter_gpu.jl Do you think you can work on a GPU filtering in latest Julia?
If we have a working imfilter_gpu we can plug it back as an option to handle these 3D cases. The code is currently commented out, sitting in the source tree until someone finds time to fix it.
Do you have experience with the CUDA stack for GPU programming?
@hameye
@juliohm For the time being, I just need a functional example of 3D image quilting, even your example might work. I need it to generate a pattern of meanders on a pointset.
I've never done GPU computation before, but this might be a good reason to have a look at it !
Júlio Hoffimann
@juliohm
The 3D example should be exactly equal to the 2D examples in the docs. You just need to inform an additional tile size for the Z dimension.
Let me know if you encounter problems. The image quilting needs a review of the cut algorithm, I remember the last time I ran it in Julia v1.6 the overlap area was patchier than usual.
Júlio Hoffimann
@juliohm
Dear @/all , we are deprecating the Gitter channel in favor of the Zulip channel, which is much richer in terms of Julia syntax highlight, LaTeX rendering, etc. There we can share images of our use cases more easily, exchange snippets of code with each other more quickly, and setup video calls to meet each other: https://julialang.zulipchat.com/#narrow/stream/276201-geostats.2Ejl I will keep this channel open while we slowly migrate to Zulip. As soon as you join the Zulip channel, press Shift+? to learn about all keyboard shortcuts. They are very useful and may help you understand the Zulip chat model as well. Thank you for being part of this community :heart:
jaquetdiamantino
@jaquetdiamantino
Hi, I'm new to Julia and I'm trying to create a package following the steps in the developer guide in the GeoStats.jl package. When I try to use the package I'm developing, the following error appears. Does anyone know where I'm going wrong? julia> Pkg.add("MyTest")
ERROR: The following package names could not be resolved: