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
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!
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.
@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.
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?
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)).