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?
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))
In your script P is a problem
@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?
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.
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
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.
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.
New paper is out (PDF version is better): https://www.frontiersin.org/articles/10.3389/fams.2021.689393/full
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/
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
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.
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)
1 reply
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
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.
Can I intject something with fault intersection check via neighborhood parameter, or I must write cusom Kriging Estimator?
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.
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
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.
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