## Where communities thrive

• Join over 1.5M+ people
• Join over 100K+ communities
• Free without limits
##### Activity
@hameye
@juliohm I created two different conventions : one is working, but the other is not .. I manage to create one, meaning I understand how to do it, but I cannot make the second work ... I have not spent more time on it since I managed to do one.
anton-degterev
@anton-degterev
Hi @juliohm , The problem I am solving is the following. I have an analytically defined function and I want to calculate the error of its prediction in each of the model cells. It seemed to me that the most efficient way to get such residuals is to get the coordinates of the cell centers and calculate the values of the original function at these points
Júlio Hoffimann
@juliohm
Hi @anton-degterev , in that case I think you don´t need to construct a full table. I would just do something like @. sol.z - f(centroid(sol.geometry))where the function f is defined for any point, for example f(p::Point) = sum(coordinates(p)). The syntax sol.z returns the column zof the solution as a vector. And the syntax sol.geometry returns the quadrangles as a vector. You may need to update to the latest Meshes.jl to get this syntax.
Júlio Hoffimann
@juliohm
What if we could learn functions directly on 3D meshes, which are discretizations of geospatial domains? In the example below, I am using GeoStats.jl to (1) simulate two correlated Gaussian processes on the statue of Beethoven and to (2) learn a function directly on this mesh with geostatistical learning (Hoffimann et al 2021). I will be sharing more details in my JuliaCon2021 talk next month.
Júlio Hoffimann
@juliohm
@hameye
Thanks @juliohm ! Great reminder that GeoStatistics are not statistics applied to geospatial objects !!
Unfortunately it has been forgotten the past years when trying do to regression/classification ...
Júlio Hoffimann
@juliohm
Exactly @hameye !
Júlio Hoffimann
@juliohm
Denis Laprise
@nside
Hi I'm curious if anyone tried or is aware of any attempt to compile/run GeoStats.jl on webassembly?
Júlio Hoffimann
@juliohm
@nside GeoStats.jl is fully written in Julia so the question you may want to ask is if anyone tried Julia on webassembly. You will have better luck asking in the Julia Discourse forum.
Júlio Hoffimann
@juliohm
FOSS4G is the largest international conference on geospatial open source software. I invite everyone to participate in the GeoStats.jl workshop therein. We will cover the basics of the Julia programming language and the powerful stack we are building for advanced geostatistics: https://callforpapers.2021.foss4g.org/foss4g-2021-workshop/talk/RHA3TH/
Markus Zechner
@MarkusZechner

Hey!
I hope you guys do well!
I have a question: Is SGS or FFTGS working in 3D? I would need conditional as well as unconditional.

I tried to construct a 3D grid with ‘CartesianGrid’ and fed it into a Simulation problem - that worked but when I tried to run SGS or FFTGS I am getting an error:

Can you maybe point me at a working example?

Thanks a lot!

Markus

Júlio Hoffimann
@juliohm
Hi @MarkusZechner , both work in 3D but FFTGS is currently unconditional only The SGS is being refactored currently to handle change of support. It should be out soon.
7 replies
Júlio Hoffimann
@juliohm
Talk at the Geostats2021 conference: https://youtu.be/eytHg5jU7r4
Markus Zechner
@MarkusZechner
Júlio Hoffimann
@juliohm
Can you try update to v0.26? We recently changed some internals to handle change of support and may have introduced a bug in v0.25 series. Try to Pkg.update()
Markus Zechner
@MarkusZechner
Perfect! Now it works. Thanks!
yizhengw
@yizhengw
Hi @juliohm! Amazing work. Thank you so much! I am wondering if you are planning to expand FFTGS to conditional as well?
Júlio Hoffimann
@juliohm
@yizhengw thanks. We could provide a helper option to do the Kriging trick on the unconditional realizations to transform them to conditional realizations. You can do this manually today with the Kriging solver with some auxiliary script
yizhengw
@yizhengw
Thanks!! I'll try to implement in the next days - would be nice if i could show you for verification.
Júlio Hoffimann
@juliohm
Of course. Feel free to share here or on Zulip.
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
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: