Where communities thrive


  • Join over 1.5M+ people
  • Join over 100K+ communities
  • Free without limits
  • Create your own community
People
Activity
    Eduardo Lenz
    @CodeLenz
    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
    geostats-workshop.jpg
    Our first workshop at a conference was really nice :heart:
    Looking forward to the next one
    Hadrien Meyer
    @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 :)
    Hadrien Meyer
    @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.
    Hadrien Meyer
    @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.
    beethoven.png
    Júlio Hoffimann
    @juliohm
    New paper is out (PDF version is better): https://www.frontiersin.org/articles/10.3389/fams.2021.689393/full
    Hadrien Meyer
    @hameye
    Thanks @juliohm ! Great reminder that GeoStatistics are not statistics applied to geospatial objects !!
    Unfortunately it has been forgotten the past years when trying do to regression/classification ...
    Júlio Hoffimann
    @juliohm
    Exactly @hameye !
    Júlio Hoffimann
    @juliohm
    JuliaCon2021 :tada: https://youtu.be/75A6zyn5pIE
    Denis Laprise
    @nside
    Hi I'm curious if anyone tried or is aware of any attempt to compile/run GeoStats.jl on webassembly?
    Júlio Hoffimann
    @juliohm
    @nside GeoStats.jl is fully written in Julia so the question you may want to ask is if anyone tried Julia on webassembly. You will have better luck asking in the Julia Discourse forum.
    Júlio Hoffimann
    @juliohm
    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/
    Markus Zechner
    @MarkusZechner
    Screen Shot 2021-08-13 at 4.09.02 PM.png

    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

    Júlio Hoffimann
    @juliohm
    Hi @MarkusZechner , both work in 3D but FFTGS is currently unconditional only The SGS is being refactored currently to handle change of support. It should be out soon.
    7 replies
    Júlio Hoffimann
    @juliohm
    Talk at the Geostats2021 conference: https://youtu.be/eytHg5jU7r4
    Markus Zechner
    @MarkusZechner
    image.png