juliaearth / geostats.jl Goto Github PK

In the example AnisotropicModels (accessible through GeoStats.examples()), the following example does not work:

using Plots; gr()

dim, nobs = 2, 50

X = 100*rand(dim, nobs)
z = rand(nobs)

scatter(X[1,:], X[2,:], c=z, label="data")

The error is the following:

MethodError: no method matching plot_color(::Float64)
Closest candidates are:
  plot_color(::Any, ::Void) at C:\Users\tomasi\.julia\v0.6\PlotUtils\src\colors.jl:15
  plot_color(::Any, ::Number) at C:\Users\tomasi\.julia\v0.6\PlotUtils\src\colors.jl:18
  plot_color(::Symbol) at C:\Users\tomasi\.julia\v0.6\PlotUtils\src\colors.jl:6

Removing the keyword c=z makes the error disappear, although the plot might look different from what was intended:

scatter(X[1,:], X[2,:], label="data")

Version information:

julia> versioninfo()                                                
Julia Version 0.6.2                                                 
Commit d386e40c17* (2017-12-13 18:08 UTC)                           
Platform Info:                                                      
  OS: Windows (x86_64-w64-mingw32)                                  
  CPU: Intel(R) Core(TM) i7-4702MQ CPU @ 2.20GHz                    
  WORD_SIZE: 64                                                     
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)  
  LAPACK: libopenblas64_                                            
  LIBM: libopenlibm                                                 
  LLVM: libLLVM-3.9.1 (ORCJIT, haswell)                             

julia> Pkg.status("GeoStats")                                       
 - GeoStats                      0.6.0              a18650da (dirty)

I have global (as in, spanning the globe) rock data that is sparse, given as 3 vecrtors (longitudes x, latitudes y, and values z), and I want to try to use some Kriging to interpolate/extrapolate these onto a global grid. Here is my current code:

using GeoStats, Plots, Distances
x = [81.45,81.42,83.35,85.24,89.51,91.01,93.05,93.07,96.04,94.09,95.03,310.0,31.0,324.0,301.0,285.0,12.0,359.0,90.0,67.0,236.0,284.0,271.0,292.0,123.0,225.0,107.0,237.0,126.0,126.0,127.5,129.0,129.0,12.56,13.23,10.23,303.0,307.0,314.0,301.0,305.0,312.0,322.0,194.0,192.0,191.0,177.0,352.99,349.75,343.61,342.84,341.85,342.83,341.43,338.93,346.0,341.7,335.5,342.5,9.0,358.55,355.75,351.45,347.79,346.92,344.42,342.73,343.7,342.0,347.65,155.39,153.29,154.08,153.45,152.0,152.2,139.17,310.44,310.52,312.0,309.93,310.47,352.41,314.85,299.72,12.36,255.64,214.41,121.32,130.58,254.7,278.76,103.07,269.0,92.0,110.0,119.0,354.0,342.0,350.0,110.0,228.3,276.1,264.4,161.4,240.82,302.0,130.13,107.11,106.43,88.57,353.0,341.5,337.92,344.94,246.4,286.27,300.36,299.65,300.63,300.42,311.6,332.0,338.0,353.0,332.0,5.25,352.0,6.8,346.01,326.0,333.0,284.0,288.0,106.4,49.2,182.13,182.7,183.18,277.0,73.3,60.0,28.0,23.0,136.0,43.0,124.2,309.5,286.45,286.6,289.76,358.0,132.65,115.63,115.4,119.0,145.05,275.0,265.0,262.66,149.3,44.0,49.0,150.0,290.0,288.0,168.0,241.0,226.0,34.5,301.0,289.7,121.9,288.0,326.0,332.0,333.0,160.0,280.0,85.0,33.0,127.0,140.0,170.0,150.0,106.0,122.0,298.0,69.0,180.0,147.0,43.0,40.0,140.0,145.0,83.0,171.0,115.0,25.0,323.0,55.5,356.8,352.2,351.3,349.5,346.25,351.5,35.0,5.0,12.0,36.0,36.2,41.0,42.5,48.0,43.15,314.58,304.68,324.95,319.5,287.88,300.0,8.0,132.0,132.0,128.0,146.15,172.4,289.1,38.85,32.6,32.87,18.45,13.28,11.3,2.9,280.0,240.0,170.0,120.0,60.0,20.0,320.0,37.0]
y = [9.39,11.32,14.42,16.38,18.28,18.21,17.11,14.3,12.46,11.45,10.12,0.0,32.0,-11.0,-34.5,11.2,-6.0,53.0,21.0,24.0,46.3,40.3,30.0,49.0,31.0,70.0,9.0,37.5,37.0,35.0,34.5,35.0,37.0,-29.42,-25.36,-14.53,64.0,62.0,60.0,59.0,55.0,59.5,60.0,54.0,54.0,54.0,53.0,34.2,32.24,25.02,23.44,21.2,19.29,13.5,4.55,28.0,17.0,12.5,8.0,42.0,34.41,31.58,28.59,22.44,20.32,18.09,14.38,14.08,14.0,11.43,50.51,48.59,49.56,48.23,47.0,-31.5,-35.1,65.07,65.1,7.5,65.1,65.08,45.04,43.29,-52.55,-12.12,18.13,54.37,2.42,36.39,19.06,-7.35,-7.21,13.0,5.0,7.0,23.0,3.0,14.0,36.0,35.0,46.3,-3.6,4.1,46.3,47.09,-35.0,-9.06,-8.25,-8.7,-6.02,62.0,63.4,60.57,46.3,27.0,-16.92,14.72,13.03,12.55,11.03,3.4,69.0,65.0,62.0,78.0,60.5,71.0,65.0,72.9,67.0,72.3,-46.2,-38.0,-7.9,-67.37,-18.57,-18.5,-20.13,3.0,0.4,-10.0,-32.0,-33.0,-36.0,15.0,-8.3,13.25,-49.75,-52.33,-54.95,58.0,-1.2,-33.32,-34.42,-20.0,-4.08,10.0,13.0,19.83,-5.0,-12.0,-19.0,-5.0,-65.0,-66.0,-15.0,37.62,53.0,-20.0,-63.5,-23.85,20.25,-42.0,65.0,70.5,71.0,-10.0,13.0,-59.0,-27.0,-8.0,35.0,54.0,10.0,-6.0,-10.0,16.0,-49.0,-80.0,-43.0,12.0,20.0,27.0,15.0,-65.0,-45.0,25.0,67.0,66.0,-21.0,36.1,36.8,35.8,45.4,56.5,47.5,26.0,43.0,45.0,37.0,37.0,17.5,11.6,12.0,11.6,-60.95,-66.05,-74.4,-77.15,69.0,66.5,68.0,-4.0,-7.0,-9.0,-17.3,-41.0,-69.5,-23.37,-31.5,-29.63,-35.78,-30.58,-25.5,-50.58,-70.0,-70.0,-70.0,-65.0,-65.0,-70.0,-75.0,25.0]
z = [-16.2,-16.1,-13.9,-15.1,-12.6,-12.8,-8.6,-9.3,-11.3,-9.2,-10.0,-9.2,-3.3,-12.9,-10.3,-8.3,-16.1,-13.5,-15.7,-12.2,-6.6,-11.3,-10.9,-5.3,-10.7,-14.3,-9.5,-3.6,-16.9,-13.2,-18.4,-13.4,-23.9,-11.0,-12.2,-18.8,-27.5,-24.1,-16.5,-17.0,-27.7,-8.8,-7.9,7.4,7.6,7.2,8.5,-11.7,-13.0,-15.3,-14.4,-14.8,-14.3,-15.7,-18.6,-12.8,-14.6,-12.0,-14.4,-13.0,-11.8,-13.6,-17.1,-17.8,-17.9,-17.1,-16.2,-14.3,-14.5,-12.1,7.6,7.3,7.7,9.6,10.1,-1.3,-6.0,-33.0,-24.2,-10.2,-34.3,-26.1,-9.9,-24.5,-4.4,-24.9,4.2,1.6,4.2,-9.3,1.3,-5.0,-13.8,2.0,-13.8,-10.6,-9.6,-9.1,-13.6,-8.5,-9.7,-6.1,-4.3,-3.4,-0.5,-6.0,-9.0,-9.3,-3.1,-2.8,-7.7,-6.5,7.6,3.4,-11.6,7.6,-3.2,-9.1,-9.6,-11.6,-13.3,-10.7,4.0,8.0,7.0,-13.0,-14.0,5.0,-14.0,-13.7,-38.0,-29.8,5.6,2.5,-2.6,-38.2,9.2,8.8,7.5,9.8,6.4,7.3,-9.2,-9.5,-21.1,4.1,-1.9,-12.0,-1.4,1.1,9.8,-10.6,4.1,1.7,-4.6,-15.1,7.1,7.3,5.9,1.5,8.3,3.4,-22.1,8.3,1.1,-3.2,7.4,-0.5,9.6,-14.5,5.7,5.5,-2.8,2.1,-38.0,4.0,-29.8,1.8,4.1,-0.5,-30.0,-12.4,7.7,7.4,7.4,0.0,-3.3,5.5,0.0,-20.8,-5.2,5.0,5.0,2.3,6.7,-8.0,-5.3,-10.2,-19.4,-26.0,4.0,-10.1,-8.4,-11.8,-11.6,-11.1,-11.6,-9.1,-9.7,-10.8,-6.2,-6.3,-4.2,6.2,10.0,6.3,-6.2,-5.7,-7.7,-10.4,-23.1,-25.8,-15.1,-8.7,-8.1,-9.5,-6.3,1.2,5.5,-15.5,-13.4,-15.3,-11.2,-11.6,-11.6,-4.8,-4.1,-3.7,-6.8,-14.3,-18.8,-14.0,-8.0,-11.2]
# list of properties with coordinates
props = OrderedDict(:prop => z)
coord = [(lon,lat) for (lon,lat) in zip(x,y)]
# define spatial data
sdata = PointSetData(props, coord)
# the solution domain (basically a grid covering the entire globe)
sdomain = begin
    dims = (180, 91)
    start = (1.0, -89.01098901098901)
    finish = (359.0, 89.01098901098901)
    RegularGrid(start, finish, dims=dims)
end
# define estimation problem for any data column(s) (e.g. :precipitation)
problem = EstimationProblem(sdata, sdomain, :prop)
# choose a solver from the list of solvers
solver = Kriging(
    :prop => (variogram=GaussianVariogram(range=35.),
              distance=Haversine(1.0)) # use the Haversine distance because the earth is round
)
# solve the problem
solution = GeoStats.solve(problem, solver)
# plot the solution
plt2 = contourf(solution)
plt1 = scatter(sdata)
plot(plt1, plt2, layout=(2,1))

I was expecting the solution (prop mean) to look like the data (prop), but I'm probably making a mistake along the way. (First time trying your package.) I was put off by issue #61 as well (origin and spacing not used in plotting), but I don't think that's my only issue here. Could you help me figure this out?

Hello, I tried the GeoStats.examples() command from the REPL and got the following error:

julia> GeoStats.examples()
ERROR: UndefVarError: examples not defined

julia> Pkg.status("GeoStats")
 - GeoStats                      0.3.5+             master

julia> versioninfo()
Julia Version 0.6.0
Commit 903644385b* (2017-06-19 13:05 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: Intel(R) Core(TM) i5-4570 CPU @ 3.20GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, haswell)

edit - The function is perhaps not exported?

Examples of usage would be great for users without background in Julia.

Tests for package run successfully, but several functions appear to not be defined.

Hi

At the moment, the covariance models do not allow nugget variances. Since it is an often desired feature, I think it should be added. @juliohm : If you agree to it, please check my pull request for an example addition of nugget field in the GaussianCovariance composite type definition. Thanks.

Missing values are now well-supported by the language. The project already takes missing values into account, but current implementations can certainly be improved.

• Deprecate uses of the valid function throughout the GeoStats.jl stack, including solver packages
• Deprecate the isvalid function and assume that the only missing value is missing

Thanks for the package!

I'm in the process of adding interpolation methods to https://github.com/evetion/GeoRasters.jl using GeoStats (GeoStatsDevTools). I hope to just provide an interpolate(::GeoArray, ::Solver) method.

However, I run into trouble with Array{Union{Missing, Float64}} support. The following code works flawlessly:

z = Array{Float64}(rand(10, 10))
z[2,2] = NaN
z[3,3] = NaN

problemdata = RegularGridData(Dict(:z=>z), (0.,0.), (1.,1.))
problemdomain = PointSet([1.5 2.5; 1.5 2.5])  # center coords of two missing values
problem = EstimationProblem(problemdata, problemdomain, :z)
solve(problem, Kriging())

While the same code, using Missing support:

z = Array{Union{Missing, Float64}}(rand(10, 10))
z[2,2] = missing
z[3,3] = missing

problemdata = RegularGridData(Dict(:z=>z), (0.,0.), (1.,1.))
problemdomain = PointSet([1.5 2.5; 1.5 2.5])  # center coords of two missing values
problem = EstimationProblem(problemdata, problemdomain, :z)

fails with:

ERROR: MethodError: Cannot convert an object of type Type{Union{Missing, Float64}} to an object of type DataType

Convenient syntax for learning tasks is missing. Below is one idea:

@regress (x1,x2) \to y
@classify (x1,x2) \to y
@cluster (x1,x2) \to y

These macros could be easily added to the corresponding files for the learning tasks.

The Makie.jl ecosystem is mature enough to support plot recipes. Migrating to this new plotting ecosystem will enable high-performance 3D visualisations on the GPU.

• SimpleMesh
• CartesianGrid
• PointSet
• GeometrySet
• Views
• Partitions
• GeoStatsBase.jl
• Variography.jl

Any idea what's going on here - it seems like something about the grid is causing this error.

mi2m = 1609.
grid_spacing = 1*mi2m
grid_domain = RegularGrid((x_extent, y_extent), (ceil(min_x), ceil(min_y)),
              (1., 1.))

I was doing a mile spacing, but nothing works

Here's the full error:

InexactError: Int64(0.23931789643928889) Type at float.jl:703 [inlined] convert at number.jl:7 [inlined] setindex!(::Array{Int64,1}, ::Float64, ::Int64) at array.jl:767 solve_globally(::EstimationProblem{GeoDataFrame{Float64,2,DataFrame},RegularGrid{Float64,2},SimpleMapper}, ::Symbol, ::Dict{Symbol,NamedTuple}) at kriging.jl:252 solve(::EstimationProblem{GeoDataFrame{Float64,2,DataFrame},RegularGrid{Float64,2},SimpleMapper}, ::Kriging) at kriging.jl:151 top-level scope at none:0

Currently the Kriging system is assembled for every single estimate at a new location, which is unacceptable in terms of performance. At some point I would like to have KrigingEstimator objects that are initialized with the data configuration once and evaluated at new locations with LU factorization or similar.

@JuliaRegistrator register()

Hi there, it could be nice to add the slicing functionality to a grid when working with 3d estimation result.

I was testing Kriging my data with different variograms and found that Gaussian variogram seems to be constantly producing problems. When I used Matern (and also other types such as Spherical) variograms the results I got looked normal (see MaternExample.png). However, when I used Gaussian variogram, the results looked awful with patches of very extreme values (see GaussianExample.png ). The specific variograms I used for the two examples are:
γtheo = MaternVariogram(sill=1.35, range=175.0, nugget=0.15, distance=Haversine(6371))
γtheo = GaussianVariogram(sill=1.2, range=300.0, nugget=0.15, distance=Haversine(6371))
I am also attaching my dataset in (lon,lat,val) format. Thanks a lot!

StnRes_label.txt

Julia Version 0.6.2
GeoStats 0.5.5

The EstimationProblem(geodata, domain, :poro) does not work if the geodata's coordinates and domain data's coordinates do have the same data type.

To transform the dataset in a Gaussian space, we can define an empirical distribution

d = EmpiricalDistribution(values) and then transform to Gaussian with transform!(values, d, Normal(0,1))

It is possible to pass a column of a DataFrame instead of an Array as values?

Hi,

I installed the package in Julia (Pkg.add("GeoStats.jl)
Later, using GeoStats

And then I was running the runtests.jl, where I encountered the problem:

Covariance models: Error During Test
Got an exception of type UndefVarError outside of a @test
UndefVarError: GaussianCovariance not defined

It seems covmodels.jl is not automatically added with the package installation. I downloaded the file and ran it in the Julia terminal, and then the tests ran smoothly.

Shouldn't the definitions for covariance models be added in the package, for subsequent use in spatial predictions using kriging functions?

Please excuse if this is a silly comment.

The Tables.jl API is finally settling down. Now it is a good time to start adhering to the interface.

It seems that origin and spacing are not taken into account when plotting things onto grids.

(not-so-M)WE:

using Plots
origin = (-100.0, -200.0)
spacing = (5.0, 5.0)
X = origin[1]:100
Y = origin[2]:200
Z = [10sin(x/10) + y for x in X, y in Y]
plt1 = heatmap(X,Y,Z')
using GeoStats
Ω = RegularGridData{Float64}(OrderedDict(:Z=>Z))
RegularGridData(OrderedDict(:Z=>Z), origin .+ 1000, spacing)
plt2 = plot(Ω)
plot(plt1, plt2, layout=(2,1))

I was expecting the bottom plot to show the shifted origin and (5.0,5.0) spacing, but it does not?

With the unupdated packages like InverseDistanceWeighting.jl, the quick example of the documentation throws at

julia> # solve the problem
       solution = solve(problem, solver)
ERROR: MethodError: no method matching LinearPath()
Closest candidates are:
  LinearPath(::D) where D<:AbstractDomain at /Users/benoitpasquier/.julia/packages/GeoStatsBase/gtoIM/src/paths/linear.jl:11
Stacktrace:
[1] preprocess(problem::EstimationProblem{PointSetData{Float64,2},RegularGrid{Float64,2},NearestMapper}, solver::Kriging)
  @ KrigingEstimators ~/.julia/packages/KrigingEstimators/RONev/src/solvers/kriging.jl:124

[2] solve(problem::EstimationProblem{PointSetData{Float64,2},RegularGrid{Float64,2},NearestMapper}, solver::Kriging)
  @ KrigingEstimators ~/.julia/packages/KrigingEstimators/RONev/src/solvers/kriging.jl:141

[3] top-level scope
  @ REPL[9]:2

I think this is because the changes in GeoStatsBase.jl are necessary for the latest version of KrigingEstimators.jl. So maybe the compat needs to strictly impose GeoStatsBase v0.10 instead of currently allowing older versions?

A very nice feature I would love to see is the ability to simulate variables on spatial domains that are much smaller than the domain of the input data. For example, I am often only interested in values along a fixed cross-section while my data is spread over 3 dimensions.

This works nicely with estimation problems, but I'm not sure if it is feasible for conditional simulation problems.

After some careful considerations, I've concluded that we don't want to depend on MLJBase.jl. Some design decisions behind the MLJ.jl project are sub-optimal for our research, and we are better off switching to the more lightweight MLJModelInterface.jl. Hopefully, this lightweight interface will be replaced by a nicer trait-based interface in the future as well.

When I tried to generate an empirical variogram and plot it using the same command as in the example, an error message was given:

'c' argument looks like a single numeric RGB or RGBA sequence, which should be avoided as value-mapping will have precedence in case its length matches with 'x' & 'y'. Please use a 2-D array with a single row if you really want to specify the same RGB or RGBA value for all points. WARNING: Base.γ is deprecated, use MathConstants.γ instead. likely near /Users/tianze/.julia/packages/IJulia/DL02A/src/kernel.jl:41

and the figure shown (see attached figure) is always the same no matter what nlags I choose. I am using Julia 0.7.0 and the version of Jupyter is 4.4.0. Thanks a lot!

Currently, spatial objects such as spatial data, solutions and statistics have their own plot recipes. This approach works fine, but is hard to maintain. I am opening this issue to discuss possible fixes.

A first proposal would be to define plot recipes for spatial data types only, and write plot recipes for other types in terms of this recipes. That is, to plot a spatial solution, we would find out the corresponding domain onto which it is defined and call the plot recipe for that domain. A solution on a PointSet for example could be plotted by converting the solution into a PointSetData and plotting it.

@riyadm what is your take on this? Remember we discussed it when you introduced the spatial statistics? Let me know if you have other suggestions.

Thank you @cormullion for your incredible work. I am excited to try ideas for the logo and learn from you.

First, let me describe some of the visual elements in our current logo to give you some background in our research field. The logo shown below is made of a "scatter" of colored points in space representing measurements of a physical property at a few locations in a map. For example, a group of geologists doing field work can report rock properties at a few sites. We always attach to each measurement its coordinates in the map.

Given these scarce measurements, the job of geostatistics is to fill in the map by using the spatial correlation between the sampled locations. There are many methods for achieving this goal, each of them relying on some prior information.

The first classical approach is based on simple 2-point correlation as the prior. Instead of using correlation as a function of distance h, we use another related concept that is more general known as the variogram:

which is basically "variogram(h) = 1 - correlation(h)". The variogram has a well-known symbol in the literature, it is denoted universally in textbooks with lowercase gamma: https://en.wikipedia.org/wiki/Gamma It may be useful as a visual reminder for practioners.

Another classical approach is based on training images. Instead of building a covariance model with the variogram function, we provide an image of reference that contains the patterns that we want to reproduce in the domain while filling the gaps. For many examples of such images, we have this repo: https://github.com/juliohm/GeoStatsImages.jl

In all geostatistical approaches, we can either estimate things locally or estimate the entire map globally. In the local approach, we usually define neighborhoods of points like the ellipsoid I drew in greyish color in the logo. The idea is that we fit a model only considering points inside of a neighborhood, and move this neighborhood around in order to fill the entire map.

I was thinking on a smaller, minimalist logo. Something that could reflect some of these visual elements I described above, and at the same time be easily remembered. Currently, the logo is too big making it impossible to fit in the documentation page tab in the top left. The Augmentor.jl package solved this logo size issue nicely by creating two variants of the logo, one of the top left tab, and a bigger one with the name written explicitly: https://evizero.github.io/Augmentor.jl/

I am open to any ideas, feel free to innovate = )

We are gradually moving to more complex stochastic processes, and the current solver infrastructure is not coping with the complexity planned ahead. This issue needs to be addressed in a major release after careful brainstorming.

Among other things, we need better API to express:

• co-estimation/co-simulation
• dependence relationships
• computation DAGs

There would be nice to have the BlockPartitioner a generalized to allow blocks of different side sizes. It should generalize the BlockCrossValidation method as well.

dear sir:
when I type such commands in Julia(Windows 7 system), there was an error prompted,
using GeoStats
using Distances
d = Ellipsoildal([2., 1.],[0.])
ERROR: UndefVarError: Ellipsoidal not defined.
maybe I need to install another Package?
Thanks,
jiamingtao

The sequential Gaussian simulation solver is currently not type stable. There should be lots of room for performance improvements. If someone can spot the places were types are not inferred correctly, that would be super helpful. Otherwise, I will tackle this when I find time.

Hello,

I'm looking to a solution to the following problem. With climate model data, the grids are often "curvilinear" (e.g. rotated poles, lambert conformal conic , etc...). I'd like to interpolate from one grid to another. I have access to the longitude and latitude grids (i.e. 2D array for longitude "LON" and 2D array for latitude "LAT") for both the source and destination grid. Hence, the data values coordinates [i, j] are specified by LON[i, j] and LAT[i, j].

My aim is to simply do an inverse weight distance for each points of the destination grid, using the 2D array of longitude and latitude. However, when I look at the codes, the "x" and "y" are vectors and not array. I was wondering if it's possible to use GeoStats with such 2D coordinates arrays. Perhaps it's simply a matter of unfolding the arrays to 1D vectors?

Thanks for any help!

A complete rewrite is needed before the next minor release of the project.

Julia Version 0.6.2
GeoStats 0.5.5

The plot for the solutions is current only available to regular grids, while cannot be used for point collections.

Currently the variogram plane functionality lives inside a plot recipe. We need to refactor the code to create a first-class VariogramPlane type. This should enable additional features and facilitate maintenance.

The current implementation is simulating the secondary variables over the entire domain and discarding the locations according to the primary variable. Now we should be able to rewrite it in terms of domain views.

The Kriging system is quite ill-conditioned when data locations are too close from each other. At some point, I have to fix this by adding regularization or any other preconditioning approach.

References:

1. https://arxiv.org/pdf/1602.00853.pdf
2. https://link.springer.com/article/10.1007/BF02065878

Hi @visr, @masiddique,

First of all I would like to thank you for reporting bugs with Float32, these should all be fixed by now.

Could you please share how single precision is important in your work? I want to make sure GeoStats.jl functions seamlessly with any Real types.

All the estimators, with the exception of Universal Kriging, have a stable numerical solution now. I will try to replace the current LU factorization of the Universal Kriging system by block elimination. In theory, this should improve the accuracy and decrease the differences between estimates on 32 vs. 64 bits data.

The current heuristic is dispatching the full-search algorithm even with large datasets. Perhaps it makes sense to dispatch the ball-search algorithm when the number of points exceeds 10000 for example.

Hi,

Sorry I'm here to bother you again. I'm noticing some anomalous, or say, undesired behavior of the unikrig() function when changing the float precision for the inputs. I first noticed it when working with my project data; unfortunately I cannot share it with you due to restrictions from the data provider. However, I managed to recreate the situation by slightly modifying the example code you've put in the "Documentation" section.

# create some data
import Distributions: Gamma
srand(0);rand(Gamma(2.0))
dim, nobs = 3, 1000
X1 = rand(dim, nobs); z = rand(nobs)
X2 = convert(Matrix{Float32},X1)
X3 = convert(Matrix{Float64},X2)

# target location
x₀ = rand(dim)

# define a covariance model
cov = GaussianCovariance(1.,1.) # sill and range

# estimation
μ1, σ² = unikrig(x₀, X1, z, cov=cov)
μ2, σ² = unikrig(x₀, X2, z, cov=cov)
μ3, σ² = unikrig(x₀, X3, z, cov=cov)

(μ1, μ2, μ3)

(2.198385878270528,1.843946672568211,-2.185270204569795)

Note that each of the three predicted values are different.

Regarding μ1 and μ2:
X1 is of type Matrix{Float64}. Truncating the precision to Float32 for X2 implies X1 != X2, and therefore, we can expect that μ1 != μ2. However, the user should generally be aware of such a possibility.

Regarding μ2 and μ3:
The difference between these two is problematic. Increasing the precision of an existing float from 32 to 64 bits doesn't change the value (X3 should be equal to X2, except for additional 0's at the far decimal places). Therefore, the difference between μ2 and μ3 seems like a bug!

I did some preliminary troubleshooting. I think this problem has to do with the following statement (and probably others similar to it):

g = Float64[γ(norm(X[:,j]-x₀)) for j=1:n] # line 100 (kriging.jl)

The code γ(norm(X[:,j]-x₀)) for j=1:n may not yield the same values when X is float32 or float64. Please look into this issue when you've some time. Thanks.

Add a section to the documentation about the GslibIO.jl package. It may help people coming from a GSLIB background to get going with the project.

Hi!

I was thinking it'd be useful to have a set of simple methods for solution postprocessing, e.g.

mean(solution)
var(solution)
quantile(solution, p)

Some of the other possibilities: mean above/below a threshold, probability (frequency) above/below a threshold.

These are simple, yet provide a nice and quick way of looking at a summary statistic of a solution. Some of these methods may be extended from Statistics module.

Would also be great to hear feedback from users on which methods would be useful. Thanks!

Sequential Gaussian simulation is not reproducing the imposed variogram structure. Below is a comparison with direct Gaussian simulation:

using GeoStats
using DirectGaussianSimulation
using Plots; pyplot(size=(1000,600))

problem = SimulationProblem(RegularGrid{Float64}(20, 20), :z => Float64, 3)

solver1 = DirectGaussSim(:z => (variogram=GaussianVariogram(range=10., nugget=0.001),))
solver2 = SeqGaussSim(:z => (variogram=GaussianVariogram(range=10., nugget=0.001),))

compare([solver1,solver2], problem, VisualComparison())

The plot shown in the Quick Example section of the documentation is misleading: if the user runs the code, it will not produce the same plot, as no Gaussian variogram is used. This should be explicitly stated in the documentation, or the plot changed.

Update outdated URLs due to migration to JuliaEarth organization.

Spatial data and domain types are implemented in terms of generic coordinate types. More tests are necessary to ensure that Unitful.jl types propagate correctly in algorithms.

The current implementation of paths as iterators is sub-optimal. It makes it difficult to construct parallel paths over multiple spatial objects. We should refactor the spatial path interface, and consider transducer-based implementations.

The Kriging-based solvers implemented so far assume that the measurements have a fixed volume. This assumption is quite strong, and we need to generalize the implementations to achieve block Kriging.

First we need to generalize the spatial object types to include volumetric information. If anyone has suggestions on efficient data structures for representing spatial domains with volumes, please feel free to share. I am tempted to continue with types a la VTK. After that, we can start working on the Kriging-based solvers to incorporate these modifications.

The implementation is almost complete. There are two items missing:

• Gaussian transform (waiting for JuliaStats/Distributions.jl#661)
• Consistent interface for reporting singular Kriging matrices (waiting for Julia v0.7)

Hi,

I'm new to julia and this package, though I was wondering is there any gpu implementation available yet?
From what I understand implementation would need specific libraries
https://juliacomputing.com/domains/gpus.html

I'm trying to apply ordinary kriging on a problem that has ~ 100.00 observations (and ~200/1000 observations in range) and was looking for best scaling implementation (besides sampling or reducing the observations), if you have any advice on this since you mention it in your documentation?

Thanks!

GeoStats v 0.4.3
Julia v 0.6.0

DomainError is raised for certain point pairs x, y passed to HaversineDistance.

Error happens here: GeoStats/src/distances.jl:138
More specific: DomainError() in sqrt at base/math.jl:425

136  # haversine formula
137  a = sin(Δφ/2)^2 + cos(φ₁)*cos(φ₂)*sin(Δλ/2)^2
138  c = 2atan2(√a, √(1-a))

a gets larger than 1 in some cases, which makes √(1-a) complain

a is very close to 1 (1.0000000000000002), I think it is a rounding error.

MWE:

using GeoStats
dist = HaversineDistance()
dist([1., -15.625], [-179, 15.625])

Some other point pairs that cause error (the pattern is obvious):
[7.125, -15.625]
[-172.875, 15.625]

[-62.375, -79.875]
[117.625, 79.875]

[-54.125, -82.875]
[125.875, 82.875]

[-124.375, -84.125]
[55.625, 84.125]