## Where communities thrive

• Join over 1.5M+ people
• Join over 100K+ communities
• Free without limits
##### Activity
Jing Wang
@immaryw
I see. Thanks a lot :)
@hameye
@juliohm It turns out angles convention of GSLIB is more far-fetched (for a geologist) than I thought, I have not spent enought time to figure it out yet
Júlio Hoffimann
@juliohm
@hameye did you have any luck with the other implemented conventions? Angle conventions are a nightmare and I really hope we can do a better job than commercial software out there which simply pick a convention and ignore the rest of the world. The plan is to write angles conventions that are easy to convert back and forth. We have most of them implemented already but can´t convert yet.
anton-degterev
@anton-degterev
Hi everyone! The question arose: how in the actual version of GeoStats to get the coordinates and values from the cell centers of the resulting spatial solution? Around February, I used for this the command "coordinates()", but now this approach does not work, no information can be found on this command and what it was replaced with. What should it be replaced with?
Júlio Hoffimann
@juliohm
Hi @anton-degterev , if you have an ensemble as the solution to a simulation problem, you can index the ensemble to get specific realizations. For example, r = ensemble[1] will give you the first realization of the process. You can then retrieve the domain of this realization with d = domain(r). Let's assume that d is a CartesianGrid{2}. This domain has finite elements, in this case quadrangles. For example, d[1] is the first quadrangle. You can get the centroid of all quadrangles with centroid.(d) and the result will be a vector of Point. Finally, you can get the coordinates of these points in a global reference system using coordinates.(centroid.(d)).
anton-degterev
@anton-degterev
@juliohm, thank you, it works. But as a result of all these manipulations, I only get a set of coordinates, but not the values in these coordinates. How is it proposed to get the values at these points?
Júlio Hoffimann
@juliohm
@anton-degterev you can get the values(r) as a table for example.
anton-degterev
@anton-degterev
But values(r) are in the CartesianGrid domain. How is it proposed to extract them in order to map them to the coordinates?
Júlio Hoffimann
@juliohm
The values are associated to the centroids of the quadrangles of the domain. I'm not sure I followed the question.
anton-degterev
@anton-degterev
Data type of values(r) is EstimationProblem{GeoData{PointSet{2, Float64}, TypedTables.Table{NamedTuple{(:X,), Tuple{Float64}}, 1, NamedTuple{(:X,), Tuple{Vector{Float64}}}}}, CartesianGrid{2, Float64}, 1}. Data type of coordinates.(centroid.(domain(r))) is Vector{StaticArrays.SVector{2, Float64}} (alias for Array{StaticArrays.SArray{Tuple{2}, Float64, 1, 2}, 1}). How can we get the values in the given coordinates by having values on a regular grid and an array of coordinate pairs?
Júlio Hoffimann
@juliohm
There is something wrong in your example. Can you please share a minimum working example?
anton-degterev
@anton-degterev
data = rand(Uniform(-1, 1), 100, 3)
sdata = georef((X=data[:,3],), collect(zip(data[:,1],data[:,2])))
P   = EstimationProblem(sdata, CartesianGrid((101,101),(-1.,-1.),(0.02,0.02)), :X)
g = ExponentialVariogram(sill=0.9, range=1.0, nugget=0.0)
LU = Kriging(:X => (variogram=g,))
sol = solve(P, LU)
typeof(coordinates.(centroid.(domain(sol))))
typeof(values(P))
Júlio Hoffimann
@juliohm
@anton-degterev you are taking the values of a problem object
In your script P is a problem
anton-degterev
@anton-degterev
@juliohm Thank you, everything works now. Now the shortest entry to get the full XYa array is hcat(getindex.(coordinates.(centroid.(domain(sol))), [1 2]),values(sol).X) This is a bit more cumbersome compared to what it was before, but also acceptable. I understand correctly, the increase in the complexity of the procedure for obtaining coordinates is a consequence of the transition from RegularGrid to CartesianGrid?
Júlio Hoffimann
@juliohm
Hi @anton-degterev , notice that every geospatial data object is already a table. If you are using Pluto as your "IDE" of development it should already be showing you the solution as a table. If you are using the normal Julia REPL, you can pipe the solution to a DataFrame sol |> DataFrame. You will see that you have almost everything you need already, including a special column named geometry with the quadrangles. So the only missing ingredient is a conversion from quadrangle to coordinates of centroid of quadrangle. There are many ways to do this conversion. You can operate with the DataFrame directly for example. I myself like Query.jl to transform the solution into other formats sol |> @mutate(geometry=centroid(_.geometry)) |> DataFrame. Perhaps if you share your end goal, we can provide a helper function.
You are correct that the increase in complexity is related to the more sophisticated domains from Meshes.jl We are working on more advanced methods that require this complexity, but trying to hide it from the end user whenever possible. So if you have a use case that justifies a helper function, we can add that for future use.
@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.
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