# Recommended workflow and packages for computing a slope along a tramway track

**URL:** <https://discourse.julialang.org/t/recommended-workflow-and-packages-for-computing-a-slope-along-a-tramway-track/119975>\
**Category:** Geo\
**Tags:** slope\
**Created:** [September 27, 2024, 8:09pm UTC](https://discourse.julialang.org/t/recommended-workflow-and-packages-for-computing-a-slope-along-a-tramway-track/119975 "2024-09-27T20:09:54Z")\
**Posts on this page:** 20\
**Page:** 2

<div class="post-metadata">

**Author:** ![stepanoslejsek](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stepanoslejsek/32/218476_2.png) [@stepanoslejsek](https://discourse.julialang.org/u/stepanoslejsek)\
**Post date:** [October 22, 2024, 12:33pm UTC](https://discourse.julialang.org/t/recommended-workflow-and-packages-for-computing-a-slope-along-a-tramway-track/119975/21 "2024-10-22T12:33:04Z")

</div>

Good, I’ll add some tests and check the fitting function and then I’ll open up a PR

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [October 23, 2024, 11:08am UTC](https://discourse.julialang.org/t/recommended-workflow-and-packages-for-computing-a-slope-along-a-tramway-track/119975/22 "2024-10-23T11:08:41Z")

</div>

@stepanoslejsek we merged your PR yesterday. We are working on some minor internal changes before releasing the new feature. Thank you ❤

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [October 23, 2024, 3:00pm UTC](https://discourse.julialang.org/t/recommended-workflow-and-packages-for-computing-a-slope-along-a-tramway-track/119975/23 "2024-10-23T15:00:45Z")

</div>

@stepanoslejsek we fixed some typos in your PR (swapped optimization parameters) and released a patch. Please let us know if you need additional help.

---

<div class="post-metadata">

**Author:** ![stepanoslejsek](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stepanoslejsek/32/218476_2.png) [@stepanoslejsek](https://discourse.julialang.org/u/stepanoslejsek)\
**Post date:** [October 23, 2024, 4:12pm UTC](https://discourse.julialang.org/t/recommended-workflow-and-packages-for-computing-a-slope-along-a-tramway-track/119975/24 "2024-10-23T16:12:48Z")

</div>

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&nbsp; Simple interpolation – Geospatial Data Science with Julia](https://juliaearth.github.io/geospatial-data-science-with-julia/11-interpolation.html) 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`.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [October 23, 2024, 5:06pm UTC](https://discourse.julialang.org/t/recommended-workflow-and-packages-for-computing-a-slope-along-a-tramway-track/119975/25 "2024-10-23T17:06:50Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![stepanoslejsek](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stepanoslejsek/32/218476_2.png) [@stepanoslejsek](https://discourse.julialang.org/u/stepanoslejsek)\
**Post date:** [October 23, 2024, 5:09pm UTC](https://discourse.julialang.org/t/recommended-workflow-and-packages-for-computing-a-slope-along-a-tramway-track/119975/26 "2024-10-23T17:09:00Z")

</div>

Yeah, the exponent is within the boundaries and I tried `GaussianVariogram` with the same result.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [October 23, 2024, 5:15pm UTC](https://discourse.julialang.org/t/recommended-workflow-and-packages-for-computing-a-slope-along-a-tramway-track/119975/27 "2024-10-23T17:15:56Z")

</div>

Appreciate if you can share a MWE. We can investigate and get back to you with further info.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [October 23, 2024, 5:28pm UTC](https://discourse.julialang.org/t/recommended-workflow-and-packages-for-computing-a-slope-along-a-tramway-track/119975/28 "2024-10-23T17:28:05Z")

</div>

> [@stepanoslejsek](#):
>
> 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`.

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

```julia
points = sample(r.geometry, HomogeneousSampling(N))

```

or sample points with a minimum distance criteria:

```julia
points = sample(r.geometry, MinDistanceSampling(...))

```

The function returns an iterator of points. You can collect and send the vector of points to the InterpolateNeighbors:

```julia
dem |> InterpolateNeighbors(collect(points), Kriging(gamma))

```

Can you please confirm the `crs(tracks.geometry)` and the `crs(dem.geometry)`?

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [October 23, 2024, 6:10pm UTC](https://discourse.julialang.org/t/recommended-workflow-and-packages-for-computing-a-slope-along-a-tramway-track/119975/29 "2024-10-23T18:10:18Z")

</div>

@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.

```julia
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")

```

 ![rail_track](https://global.discourse-cdn.com/julialang/original/3X/a/f/afe75610380e24a09087779d18d9f4353f65b036.png)

---

<div class="post-metadata">

**Author:** ![stepanoslejsek](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stepanoslejsek/32/218476_2.png) [@stepanoslejsek](https://discourse.julialang.org/u/stepanoslejsek)\
**Post date:** [October 23, 2024, 6:20pm UTC](https://discourse.julialang.org/t/recommended-workflow-and-packages-for-computing-a-slope-along-a-tramway-track/119975/30 "2024-10-23T18:20:46Z")

</div>

@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
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
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.

---

<div class="post-metadata">

**Author:** ![stepanoslejsek](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stepanoslejsek/32/218476_2.png) [@stepanoslejsek](https://discourse.julialang.org/u/stepanoslejsek)\
**Post date:** [October 23, 2024, 6:23pm UTC](https://discourse.julialang.org/t/recommended-workflow-and-packages-for-computing-a-slope-along-a-tramway-track/119975/31 "2024-10-23T18:23:04Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [October 23, 2024, 6:28pm UTC](https://discourse.julialang.org/t/recommended-workflow-and-packages-for-computing-a-slope-along-a-tramway-track/119975/32 "2024-10-23T18:28:33Z")

</div>

> [@stepanoslejsek](#):
>
> I’d like to know the elevation along the track which looks like this

Perhaps you can debug with a few ropes along the track:

```julia
ropes = tracks.geometry[1:3]

dem |> InterpolateNeighbors(ropes, model)

```

and try other interpolation models like `NN` or `IDW` to nail down the issue.

---

<div class="post-metadata">

**Author:** ![stepanoslejsek](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stepanoslejsek/32/218476_2.png) [@stepanoslejsek](https://discourse.julialang.org/u/stepanoslejsek)\
**Post date:** [October 23, 2024, 6:31pm UTC](https://discourse.julialang.org/t/recommended-workflow-and-packages-for-computing-a-slope-along-a-tramway-track/119975/33 "2024-10-23T18:31:46Z")

</div>

Both `NN` and `IDW` works fine. The kriging model does not (even for 3 ropes)

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [October 23, 2024, 6:40pm UTC](https://discourse.julialang.org/t/recommended-workflow-and-packages-for-computing-a-slope-along-a-tramway-track/119975/34 "2024-10-23T18:40:16Z")

</div>

Maybe you can execute the interpolation in debug mode to identify where the NaNs are coming from? You can do the following on VSCode:

```julia
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.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [October 23, 2024, 8:29pm UTC](https://discourse.julialang.org/t/recommended-workflow-and-packages-for-computing-a-slope-along-a-tramway-track/119975/35 "2024-10-23T20:29:57Z")

</div>

@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.

---

<div class="post-metadata">

**Author:** ![stepanoslejsek](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stepanoslejsek/32/218476_2.png) [@stepanoslejsek](https://discourse.julialang.org/u/stepanoslejsek)\
**Post date:** [October 23, 2024, 9:33pm UTC](https://discourse.julialang.org/t/recommended-workflow-and-packages-for-computing-a-slope-along-a-tramway-track/119975/36 "2024-10-23T21:33:41Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [October 23, 2024, 9:42pm UTC](https://discourse.julialang.org/t/recommended-workflow-and-packages-for-computing-a-slope-along-a-tramway-track/119975/37 "2024-10-23T21:42:15Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![stepanoslejsek](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stepanoslejsek/32/218476_2.png) [@stepanoslejsek](https://discourse.julialang.org/u/stepanoslejsek)\
**Post date:** [October 23, 2024, 9:56pm UTC](https://discourse.julialang.org/t/recommended-workflow-and-packages-for-computing-a-slope-along-a-tramway-track/119975/38 "2024-10-23T21:56:44Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [October 23, 2024, 11:10pm UTC](https://discourse.julialang.org/t/recommended-workflow-and-packages-for-computing-a-slope-along-a-tramway-track/119975/39 "2024-10-23T23:10:06Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![stepanoslejsek](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stepanoslejsek/32/218476_2.png) [@stepanoslejsek](https://discourse.julialang.org/u/stepanoslejsek)\
**Post date:** [October 24, 2024, 9:00am UTC](https://discourse.julialang.org/t/recommended-workflow-and-packages-for-computing-a-slope-along-a-tramway-track/119975/40 "2024-10-24T09:00:37Z")

</div>

That works perfectly!

[Previous page](https://discourse.julialang.org/t/recommended-workflow-and-packages-for-computing-a-slope-along-a-tramway-track/119975.md?page=1)