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?
sol |> DataFrame. You will see that you have almost everything you need already, including a special column named
geometrywith 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.
@. sol.z - f(centroid(sol.geometry))where the function
fis defined for any point, for example
f(p::Point) = sum(coordinates(p)). The syntax
sol.zreturns the column
zof the solution as a vector. And the syntax
sol.geometryreturns the quadrangles as a vector. You may need to update to the latest Meshes.jl to get this syntax.
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!
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, -d) 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)
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
domain(obsthat are consistent. Let me know if you have questions.
s1, s2, ..., skas
Pointand 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
neighborhood=Ellipsoid(...)for example: https://juliaearth.github.io/GeoStats.jl/stable/solvers/builtin.html#GeoEstimation.Kriging