I would like to ask for a recommendation on a workflow for computing a slope along a tramway track using Julia’s geo ecosystem (organizations JuliaEarth, JuliaGeo, …).
The tramway tracks are obtained from OpenStreetMaps as sequences of points (in WGS84). Using GeoIO package, we can import them into a GeoTable with its geometry as a sequence of, possibly unordered, ropes.
The elevation data is obtained from a Digital Terrain Model (DTM), in ETRS89, in triangular network (in GML format), possibly formatted also as a grid (raster, TIF format) if needed.
How can these two data sets be combined? The DTM is huge, but ultimately we are only interested in what is going on along the tramway tracks.
We do not need just the elevation but we ultimately need the slope. Should the computation of the slope be performed before the elevation data is restricted to the tramway track?
Which packages are relevant here? In particular the DTM->slope and projection to the tramway track part. What can be used for this computation? I have found that slope computation is mentioned in Geomorphometry package. Anything else?
Geomorphometry will do rasterized slope. If you have very high resolution that may be fine.
If you load your data as a Rasters.jl Raster Geomorphometry should return a Raster.
Then you can just rasterize you track into the slope raster with mask(dtrm_rast; with=tramway_linestrings) where tramway_linestring might be a shape file, geo json or whatever.
But this may not be the best method. I would instead simply convert your lines to points. Maybe points = collect(GeoInterface.getpoint(linestring)) is the easiest and most generic way to do that - it should work on any geometries in the ecosystem.
Then use Rasters.extract to extract the elevation of each node. extract(dtm_raster, points) |> DataFrame will give you a DataFrame of you points and their elevations. Then you can use the GeometryOps.distance to find the distance between each point, and with that you can use the elevation value to calculate slope.
This may work if your dtm is very high resolution and your tram lines lower or the slope is very gradual. Otherwise you may way to interpolate points into the dtm, and for that you could try Interpolations.jl or maybe something in GeoStats.jl.
If you get this written up, or close to it feel free to post some code and we can workshop it.
Hi @zdenek_hurak! It would be nice to support GML in GeoIO.jl to avoid unnecessary conversions to raster formats. In the meantime, you can try saving the data as geotiff before trying GeoIO.jl on it.
Given a geotable tracks with the tramtracks and another geotable dtm with the terrain model, you could try the syntax dtm[tracks.geometry[1], :] to extract all variables from the terrain that intersect with a geometry of the tracks.
If all your slopes are constrained to a small region where geodesics don’t play an important role, you can rely on the existing algorithms in Meshes.jl to compute slopes along the tracks. Geodesic geometry is not available yet, but is coming soon.
In theory, you should be able to target your computation near the tracks and avoid raster approximations, which may compromise location accuracy. You will likely encounter rough edges with any approach. If you need help with GeoIO.jl or any other package in the GeoStats.jl stack, feel free to reach out in a separate post.
Once interpolation is implemented in Rasters.jl (maybe once Raf finishes his PhD and get’s through his infinite list of other priorities… or @asinghvi17 maybe we can add this to our list), this will become trivial:
Hi! If I understand that syntax correctly, then dtm[tracks.geometry[1],:] will select particular point from the dtm, am I right?
If so, this is not the scenario that we are dealing with. The intersection of dtm and tracks is empty and we need to somehow interpolate the altitude (and then slope) at points given by tracks.
It might be my misinterpretation of what dtm[tracks.geometry[1],:] does.
This will return the rows of dtm that intersect with first geometry of the tracks geotable. If the first geometry is a Segment or a Rope then it will return all pixels that intersect with the Rope.
It seems like you have an easier problem then. You can look into Interpolate and InterpolateNeighbors and ask to interpolate along the tracks. Perhaps you have a MWE to share? We can look into it.
That’s certainly one option I hadn’t considered. How is that different from the kriging method? As far as I know, kriging can be used for interpolation as well. I’ll post the MWE later. The MWE consist only loading 2 geojson files using GeoIO.jl, one for the DEM and one for the tram track.
The Interpolate and InterpolateNeighbors transforms accept any geostatistical model, including Kriging, IDW, LWR, Polynomial. They are higher-level constructs to perform interpolation over general domains, not just grids.
So using InterpolateNeighbors should do the work? That sounds way easier than I thought. In that case I’ll post update on that once I try this. Thank you
I am currently in a state where I have only selected points within a certain distance from the tram track (150 m) from the digital elevation map, as shown in the first picture (the red color is the tram track).
To my limited knowledge so far, the empirical variogram should have a plateau. Thus, a linear trend seems to be evident. After removing the linear trend using
Set the maxlag to a smaller value in order to focus on small lag distances. This is the most important part of the variogram that you need to capture well.
Try to detrend with a different degree, for example Detrend(2) to see if the variogram estimates change.
After that you can estimate the cross-validation error of this Kriging model against simpler models like NN or IDW. This is explained in the docs here:
The figure shows 2 empirical variograms: with and without the detrend function called.
After reducing the distance to 30 meters and performing detrend(2), the `linear character’ seems to persist. I’ll definitely check the provided paper.
Thanks for the paper, I tried to use fit function from GeoStatsFunctions to fit the empirical variogram to PowerVariogram using GeoStatsFunctions.fit(PowerVariogram, g) but it seems that fitting a power variogram is not yet supported.
That would be great. We only have fitting methods implemented for stationary models so far. It should be easy to add a new method for PowerVariogram. Please let us know if you need help.
However, when I use only plot(γ) the variogram is shown only to 3 m, therefore the plot does not really match with plot(g), where g is the EmpiricalVariogram with h going up to 30 m.
You can do plot(gamma, maxlag=30m) to set the default maximum lag in the plot. By default we plot theoretical models up to 3 times the range. In the case of the PowerVariogram we plot up to 3m by default given that it doesn’t have a range.
