Good, I’ll add some tests and check the fitting function and then I’ll open up a PR
2 Likes
@stepanoslejsek we merged your PR yesterday. We are working on some minor internal changes before releasing the new feature. Thank you 
2 Likes
@stepanoslejsek we fixed some typos in your PR (swapped optimization parameters) and released a patch. Please let us know if you need additional help.
1 Like
Thank you, I’m glad I could help. Now I’m at the final stage - interpolating the elevation along the tram track (I have my fitted variogram and data ready to perform the interpolation).
When I follow the 11 Simple interpolation – Geospatial Data Science with Julia with my data (i.e.DEM in close neighborhood of the tram track with altitude values given by series of points and tram track given by series of ropes), my interpolation results in NaN values by calling (projected into WebMercator)
interp = dem_railway |> InterpolateNeighbors(pointify(r.geometry), Kriging(γ)) or
interp = dem_railway |> InterpolateNeighbors(r.geometry, Kriging(γ)). What can cause the interpolation to return NaN?
I understand that the second call will return the interpolated value at the center of the rope, which is not the desired result, therefore I call pointify on r.geometry.
Can you try with other variogram models to see if this is specific to PowerVariogram? Did you pick an exponent < 2 as explained in the paper?
Yeah, the exponent is within the boundaries and I tried GaussianVariogram with the same result.
Appreciate if you can share a MWE. We can investigate and get back to you with further info.
That is correct. By default the interpolation results refer to the centroid of the geometries in the domain. If you need to interpolate at other geometries you have some options.
You can sample the domain homogeneously based on the measure of the geometries with
points = sample(r.geometry, HomogeneousSampling(N))
or sample points with a minimum distance criteria:
points = sample(r.geometry, MinDistanceSampling(...))
The function returns an iterator of points. You can collect and send the vector of points to the InterpolateNeighbors:
dem |> InterpolateNeighbors(collect(points), Kriging(gamma))
Can you please confirm the crs(tracks.geometry) and the crs(dem.geometry)?
@stepanoslejsek In case you are interested. I use the file that you sent me and simulated this exercise with SRTM01, since I don’t have your LIDAR data. This is what I got. With your data the profile should be smoother.
using GMT
# Read the rail file
rail = gmtread("railway_tram7.geojson");
# Compute the 150 m buffer zone
bf = buffergeo(rail[1], width=150);
# Simulate your LIDAR data with SRTM01
xyz = grd2xyz("@earth_relief_01s", region=bf);
# Retain only the points inside the buffer zone (not really needed)
xyz_rail = gmtselect(xyz, polygon=bf);
# Compute a grid at ~30 m resolution with minimum curvature
G = gridit(xyz_rail, method=:surface, inc=0.0003);
# Extract the track from the grid (returns a lon lat z)
track = grdtrack(G, rail[1]);
# Compute the accumulated distance along the track
rail_track = mapproject(track, track_distances=(accumulated=true,), o="3,2");
# Viz it.
viz(rail_track, title="Rail track bumped by SRTM 1 arc sec", name="rail_track")
1 Like
@juliohm Crs of both is
WebMercator{WGS84Latest, CoordRefSystems.Shift{Quantity{Float64, NoDims, Unitful.FreeUnits{(°,), NoDims, nothing}}, Quantity{Float64, 𝐋 , Unitful.FreeUnits{(m,), 𝐋 , nothing}}}(0.0°, 0.0 m, 0.0 m), Quantity{Float64, 𝐋 , Unitful.FreeUnits{(m,), 𝐋 , nothing}}}
The DEM geometry looks like this
julia> dem.geometry
2500998 view(::PointSet, [1168390, 1011504, 1168334, 1011445, ..., 234287, 182630, 182956, 293135])
├─ Point(x: 2.0283759246077256e6 m, y: 6.4071342973804595e6 m)
├─ Point(x: 2.0283833377949062e6 m, y: 6.407126130970184e6 m)
├─ Point(x: 2.0283830182996714e6 m, y: 6.407134193308169e6 m)
├─ Point(x: 2.0283900836383237e6 m, y: 6.407129096925305e6 m)
├─ Point(x: 2.028389762812237e6 m, y: 6.4071331582979355e6 m)
⋮
├─ Point(x: 2.0215452846179784e6 m, y: 6.415756048713601e6 m)
├─ Point(x: 2.0215759514491577e6 m, y: 6.415701598773762e6 m)
├─ Point(x: 2.0215760135708791e6 m, y: 6.415709699428169e6 m)
├─ Point(x: 2.021576859022117e6 m, y: 6.415717285723413e6 m)
└─ Point(x: 2.021583630664356e6 m, y: 6.4157023605668e6 m)
These points are used to construct the variogram (for each point I know the altitude).
I’d like to know the elevation along the track which looks like this
julia> tracks.geometry
34 GeometrySet
├─ Rope((x: 2.02857e6 m, y: 6.40726e6 m), (x: 2.02857e6 m, y: 6.40727e6 m))
├─ Rope((x: 2.02857e6 m, y: 6.40727e6 m), ..., (x: 2.02875e6 m, y: 6.40896e6 m))
├─ Rope((x: 2.02875e6 m, y: 6.40896e6 m), (x: 2.02881e6 m, y: 6.40918e6 m))
├─ Rope((x: 2.02881e6 m, y: 6.40918e6 m), ..., (x: 2.02889e6 m, y: 6.40946e6 m))
├─ Rope((x: 2.02889e6 m, y: 6.40946e6 m), ..., (x: 2.02965e6 m, y: 6.41079e6 m))
⋮
├─ Rope((x: 2.02338e6 m, y: 6.41731e6 m), (x: 2.02337e6 m, y: 6.41732e6 m))
├─ Rope((x: 2.02337e6 m, y: 6.41732e6 m), ..., (x: 2.02328e6 m, y: 6.41735e6 m))
├─ Rope((x: 2.02328e6 m, y: 6.41735e6 m), ..., (x: 2.02258e6 m, y: 6.41763e6 m))
├─ Rope((x: 2.02258e6 m, y: 6.41763e6 m), ..., (x: 2.02151e6 m, y: 6.41554e6 m))
└─ Rope((x: 2.02151e6 m, y: 6.41554e6 m), ..., (x: 2.02137e6 m, y: 6.41562e6 m))
As mentioned above, using pointify(tracks.geometry) will be used in InterpolateNeighbors function.
Wow, this looks very promising! Looking at the elevation from noisy GPS measurements, the altitude profile looks similar to yours. I’ll try to incorporate the DEM.
Perhaps you can debug with a few ropes along the track:
ropes = tracks.geometry[1:3]
dem |> InterpolateNeighbors(ropes, model)
and try other interpolation models like NN or IDW to nail down the issue.
1 Like
Both NN and IDW works fine. The kriging model does not (even for 3 ropes)
Maybe you can execute the interpolation in debug mode to identify where the NaNs are coming from? You can do the following on VSCode:
interp = InterpolateNeighbors(ropes, model)
@enter interp(dem)
The menu will have the instructions on how to enter into the function calls line by line. You can watch the values of intermediate variables to spot the issue.
@stepanoslejsek it may be the case that you have repeated points in the DEM table? What happens if you send the table to UniqueCoords() before the InterpolateNeighbors?
Covariance matrices with repeated rows or columns can be problematic in the Kriging model. If that solves the issue we can think of ways to improve the usability.
1 Like
Yes, you hit the spot. I removed duplicate coordinates before constructing an EmpiricalVariogram using pipe. But I haven’t done same thing for the inference of the kriging model. Maybe adding same warning as in EmpiricalVariogram would be helpful. Thanks for pointing this out.
1 Like
Glad it worked. We will think about ways to improve this issue. Either by calling UniqueCoords internally in InterpolateNeighbors and Interpolate or by emitting a warning as you suggested.
1 Like
As far as I know, kriging interpolation can provide a variance of the estimate. Is this feature built in GeoStats right now? I am just curious about it.
Yes, use prob=true in InterpolateNeighbors for probabilistic interpolation. It will return a Distributions.jl object at every point. You can extract any statistic from it afterwards.
1 Like
That works perfectly!
1 Like