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.
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!
Bit of hacks, but I think this works:
view(𝒟.domain, LinearIndices(size(A))[ismissing.(A)])
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.
GeoStats.jl v0.24 is out :tada: More details on Discourse: https://discourse.julialang.org/t/ann-geostats-jl-v0-24/58585
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
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
3 replies
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!
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
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.
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.
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.
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?
OK...so I can use the previous solution on a PointSet?
I am trying to understand the proper usage.
props = (Y=vals1,);
coord = ptos;
D = georef(props,coord);
prob = EstimationProblem(D,grid,:Y);
a = 2.0
solver = Kriging(
:Y => (variogram=GaussianVariogram(range=a),)
);
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);
and everything should work fine. The estimation will be performed only on the pointset instead of in a full grid.
Our first workshop at a conference was really nice :heart:
Looking forward to the next one