Where communities thrive


  • Join over 1.5M+ people
  • Join over 100K+ communities
  • Free without limits
  • Create your own community
People
Activity
    Maarten Pronk
    @evetion:matrix.org
    [m]
    using GeoStatsBase
    using GeoEstimation
    
    A = Array{Union{Missing,Float64}}(rand(100,100))
    A[2,2] = missing  # set only 1 value to missing
    
    𝒟 = georef((z=A,), origin=(0.,0.), spacing=(1.,1.))
    𝒫 = EstimationProblem(𝒟, 𝒟.domain, :z)
    @time r = solve(𝒫, IDW())
    A_interpolated = reshape(r[:z], (100,100))
    But this ofcourse interpolates everything and I want it to just interpolate the 1 missing value. I've used CopyMapper() before. How do we approach this? Should I convert everything to pointsets?
    Júlio Hoffimann
    @juliohm
    Hi @evetion:matrix.org , I think you want the second argument of the problem to only contain the missing locations. Would that be too problematic? You can pass in a view(domain(𝒟), inds) as your second argument. This will be a solution over this view.
    Maybe this use case you have of filling gaps in a single domain deserves a new API.
    Something like 𝒟new = interpolate(𝒟, inds)
    The EstimationProblem is a more general problem where the domain of estimation doesn't necessarily match the domain of the data. Think of applications in geosciences where the data is only available along drillholes and the domain is a 3D grid model.
    Maarten Pronk
    @evetion:matrix.org
    [m]
    I can do that. I understand the estimation, this was mainly the question of how to use the estimation in very specific places, but still making use of the grid. I might move to pointsets altogether for rotated grids and actually warping (a la gdalwarp) in the future.
    Thanks for the very quick reply!
    Maarten Pronk
    @evetion:matrix.org
    [m]
    Bit of hacks, but I think this works:
    view(𝒟.domain, LinearIndices(size(A))[ismissing.(A)])
    Júlio Hoffimann
    @juliohm
    A quick comment is that 𝒟.domain is undocumented internal behavior. You should use the public API domain(𝒟) to avoid issues with views of views.
    Also, remember that in the latest release 𝒟 is a Tables.jl table. So you could find the missing rows with some table filtering operation instead of using the LinearIndices hack above.
    What we are trying to achieve is an interface where the user can work with these objects as if they were simple tables in a statistical/machine learning workflow.
    Júlio Hoffimann
    @juliohm
    GeoStats.jl v0.24 is out :tada: More details on Discourse: https://discourse.julialang.org/t/ann-geostats-jl-v0-24/58585
    Hadrien Meyer
    @hameye
    Hi @rmcaixeta, did you get the chance to work on this issue (JuliaEarth/GeoStats.jl#158) since February ? I'd be interested in the requests you mentioned
    Rafael Caixeta
    @rmcaixeta
    Hey @hameye . Nice to know. I didn't worked with that since then. There were some major refactoring and the octants/sectors needed some more adjusts. I can start the pull request again this weekend with what I had so far adapted to the new version. If you could help me with it after that, it'd be great
    Hadrien Meyer
    @hameye
    @rmcaixeta Sure, I'm glad to help !
    ammar1510
    @ammar1510
    Hi, I'm Ammar Shaikh, a sophomore from IIT ISM Dhanbad with a major in Computer Science and Engineering. I'd like to be a contributor to geostats.jl at GSoC.
    Júlio Hoffimann
    @juliohm
    Hi @ammar1510 can you please confirm that you have good programming experience in Julia?
    3 replies
    Jing Wang
    @immaryw
    Hi, I'm Jing. I just started to use GeoStats and have a question about how to set the order of each realization. I generated three processes (1 realization for each) in the order of "Z1,Z2,Z3" and I wanted them plot in the same order. But after SpecGaussSim function, the order changed and I can only plot them as "Z2,Z3,Z1". Is there any way to control the order for SpecGaussSim or change the order of the subplot? Thanks!
    image.png
    Júlio Hoffimann
    @juliohm
    Hi @immaryw , can you update to the latest version of the project? GeoStats.jl v0.24
    The SpecGaussSim solver has been replaced by the built-in FFTGS solver
    Same options as before. The ecosystem has changed dramatically recently and so if you are starting now you should try the most recent version for better support.
    In the most recent version I think we preserve the order of processes
    Jing Wang
    @immaryw
    @juliohm Thanks for your reply! I updated it to 0.24.1 and used FFTGS instead. But the order is still messy...
    image.png
    Júlio Hoffimann
    @juliohm
    Thank you @immaryw. The issue was fixed in JuliaEarth/GeoStatsBase.jl@12a29db and a patch release has been triggered with the fix. You should be able to update after this PR is merged by the bot: JuliaRegistries/General#34943 (approximately 15min)
    You can type ] update in your environment or ] update GeoStatsBase to only update the package.
    Jing Wang
    @immaryw
    Thank you very much! I updated GeoStatsBase to 0.21.2 and copied your test/solvers.jl to re-run. I also added a line "sol = solve(prob, solv)" after your test code. After running "solve(prob,solv)", it turns out "variables: Z2, Z3 and Z1". Looks like the order is still messy. But it changes back to the correct order when I plot(sol). I was wondering why this happened and hoping that it may not cause any bug in the future.
    image.png
    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.
    image.png
    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!