sol = solve(prob, solver);
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!
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))
.
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?
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))
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?
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.
@. 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 z
of 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.