## Where communities thrive

• Join over 1.5M+ people
• Join over 100K+ communities
• Free without limits
##### Activity
Júlio Hoffimann
@juliohm
Thank you @immaryw the issue is in the print method here: https://github.com/JuliaEarth/GeoStatsBase.jl/blob/92e3ad70d5e24a0e1bed825faa9bd20baa4bde2b/src/ensembles.jl#L64 I will try to find the time to fix it. Could you please open an issue on the GeoStats.jl repo so that I don't forget? You shouldn't worry about bugs related to this order, the entire solution is constructed using a dictionary so we always assign the data to the correct variable regardless of order. We only need to fix the printing methods to preserve the order that was specified in the solver.
Jing Wang
@immaryw
@juliohm Gotya! Thanks for your clarification! I opened the issue (JuliaEarth/GeoStatsBase.jl#86).
Eduardo Lenz
@CodeLenz
Dear all. First, thanks for this amazing package. I am just starting to use it and I have some (silly, for sure) doubts.
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.
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.