[GIS] Slope at a point on a Line (LineString) using Shapely and python

pythonshapely

I've got a method for making cross-section points along a shapely LineString object that I'm not totally happy with.

My strategy is:

  1. Create regular points every 50cm along a LineString
  2. Create an offset point 1cm away from the 50.1cm point
  3. Use these two points to find a slope
# Get 50cm spaced points                    
points = [linestr.interpolate(dist) for dist in np.arange(0, linestr.length, 0.5)]

# Get points at 50.1cm
_offsetpts = [linestr.interpolate(currDist+0.001) for currDist in np.arange(0, linestr.length, 0.5)]

# Use these two points to do rise/run and get a slope
slopes = [((points[idx].coords[0][1] - _offsetpts[idx].coords[0][1]) /
           (points[idx].coords[0][0] - _offsetpts[idx].coords[0][0])) for idx, pt in enumerate(_offsetpts)]

It works but it's not very elegant and I worry it will break for edge cases (like when 1 extra cm runs over the edge of the line)

Does anyone know a better way to get the slope of a line a known distance along a LineString object using python and shapely?

Best Answer

You don't need a second list of points to find the slopes

import itertools
# iterate by pairs of points
slopes = [(second.y-first.y)/(second.x-first.x)for first, second in itertools.izip(points, points[1:])]

If you want the mean slope of the slopes

import numpy as np
def mean_slope(xs,ys):
    m = (((np.mean(xs)*np.mean(ys)) - np.mean(xs*ys)) /((np.mean(xs)**2) - np.mean(xs*xs)))
    return m
print mean_slope(slopes)
-2.2934712193638029

You can also use NumPy or SciPy

x, y = zip(*[(pt.x, pt.y) for pt in points])
slope, intercept = np.polyfit(x,y,1)
print slope
-2.2934712193638029
from scipy import stats
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
print slope
-2.2934712193638029

If the LineString is straight all the values of slopes = slope

If you want the slope between two point

first = points[5]
second = points[6]
print (second.y-first.y)/(second.x-first.x) 
-2.2934712193638029 #straight line in my case
print np.polyfit([second.x,first.x],[second.y,first.y],1)[0]
-2.2934712193638029 
print stats.linregress([second.x,first.x],[second.y,first.y])[0]
-2.2934712193638029
Related Question