## Where communities thrive

• Join over 1.5M+ people
• Join over 100K+ communities
• Free without limits
##### Activity
Eduardo Lenz
@CodeLenz
I started by generating a Latin Hipercube using LatinHypercubeSampling and than creating an estimation for a simple function (x[1]-2)^2, just to test it.
The solution is as expected. My question is: how can I evaluate the approximation in a given point not in my CartesianGrid?
Júlio Hoffimann
@juliohm
Hi @CodeLenz I am not sure I understand your question. Can you clarify what you are trying to achieve? latin hypercube sampling is just a method for sampling a space
Eduardo Lenz
@CodeLenz
for example, my solution is in a variable of type GeoData{CartesianGrid{1, Float64}, .....), such that I believe I can only obtain numerical solutions in the grid previously created.
Júlio Hoffimann
@juliohm
You are trying to define an estimation problem on a point set?
Eduardo Lenz
@CodeLenz
ops...indeed @juliohm . I was just trying to explain my approach.
Júlio Hoffimann
@juliohm
If you want to perform estimation on a point set you just need to replace the CartesianGrid by the PointSet of interest
Eduardo Lenz
@CodeLenz
But...after using the LHC to set an initial support points and the associated function values, I created the grid.
OK...so I can use the previous solution on a PointSet?
I am trying to understand the proper usage.
Júlio Hoffimann
@juliohm
You can define EstimationProblem(georef(table, points), newpoints, :variable). Perhaps if you share some snippet of code it will be easier to help.
Eduardo Lenz
@CodeLenz
I just did this
props = (Y=vals1,);
coord = ptos;
D = georef(props,coord);
prob = EstimationProblem(D,grid,:Y);
a = 2.0
solver = Kriging(
:Y => (variogram=GaussianVariogram(range=a),)
);
sol = solve(prob, solver);
Júlio Hoffimann
@juliohm
Yes, you can replace the grid by a pset = PointSet(ptos2)
and everything should work fine. The estimation will be performed only on the pointset instead of in a full grid.
Eduardo Lenz
@CodeLenz
I see! Thank you very much! And thank you for this amazing package!
Júlio Hoffimann
@juliohm
You are welcome! Glad it is useful! :)
Júlio Hoffimann
@juliohm
Our first workshop at a conference was really nice :heart:
Looking forward to the next one
@hameye
Hi everyone, I am trying to use GSLIB convention, does anyone how use it when having strike dip and pitch as angles ?
Jing Wang
@immaryw

Hello everyone, I have a question about exporting/writing spatial data. After doing simulations with GeoStats.jl, I got data with the data structure "GeoData{CartesianGrid{2,Float64},TypedTables.Table{NamedTuple{(:X1,),Tuple{Float64}},1,NamedTuple{(:X1,),Tuple{Array{Float64,1}}}}}". How to export it as raster data?

I tried NCDatasets.jl, but seems like they don't provide a method to convert GeoData to netCDF? Do you happen to know the best package to export the simulation results of GeoStats.jl? Thanks!

Júlio Hoffimann
@juliohm
Hi @immaryw , GeoData is simply a table of values together with a domain. In this case the domain is a grid. You can get access to the table of values with values(geodata) or you can interpret the data as a Julia array with asarray(geodata, variable). You can then use any package to write the array or table to disk.
Hi @hameye , the GSLIB convention is super strange. First, if I recall correctly the rotations are in Z then X then Y. With the Z rotation in clockwise orientation and the others in counter-clockwise orientation. Did you manage to match the convention to your angles? Also, we should probably add a plot recipe for ellipsoids to help visualize these conventions.
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