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?
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))
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 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.
@. 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 z
of 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.
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
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)
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.
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
neighborhood=Ellipsoid(...)
for example: https://juliaearth.github.io/GeoStats.jl/stable/solvers/builtin.html#GeoEstimation.Kriging