diff --git a/src/data/profile_features.py b/src/data/profile_features.py
index f81a21a..f3da959 100644
--- a/src/data/profile_features.py
+++ b/src/data/profile_features.py
@@ -7,6 +7,8 @@ from shapely.ops import transform
 import pyproj
 from functools import partial
 import numpy as np
+from scipy.interpolate import interp1d
+from scipy import stats
 
 def shapes_from_shp(shp_file):
     """
@@ -37,6 +39,62 @@ def convert_coord_systems(g1, in_coord_system='EPSG:4326', out_coord_system='EPS
     g2 = transform(project, g1)  # apply projection
     return g2
 
+
+def get_slope(x, z, top_elevation, btm_elevation, method='end_points'):
+    """
+    Returns a slope from a bed profile
+    :param x: List of x bed profile coordinates
+    :param z: List of z bed profile coordinates
+    :param top_elevation: Top elevation of where to take the slope
+    :param btm_elevation: Bottom elevation of where to take the slope
+    :param method: Method used to calculate slope (end_points or least_squares)
+    :return:
+    """
+
+    profile_line = LineString([[x,z] for x,z in zip(x,z)])
+
+    end_points = {
+        'top': {
+            'elevation': top_elevation,
+        },
+        'btm': {
+            'elevation': btm_elevation,
+        }}
+
+    # If there are multiple intersections, take most seaward of landward point?
+    multiple_points = 'landward'
+
+    for end_type in end_points.keys():
+        elevation = end_points[end_type]['elevation']
+        elevation_line = LineString([[min(x), elevation], [max(x), elevation]])
+        intersection_points = profile_line.intersection(elevation_line)
+
+        if intersection_points.is_empty:
+            return None
+
+        if multiple_points == 'landward':
+            intersection_point = list(intersection_points)[0]
+        elif multiple_points == 'seaward':
+            intersection_point = list(intersection_points)[-1]
+
+        end_points[end_type]['x'] = intersection_point.xy[0][0]
+        end_points[end_type]['z'] = intersection_point.xy[1][0]
+
+    if method == 'end_points':
+        x_top = end_points['top']['x']
+        x_btm = end_points['btm']['x']
+        z_top = end_points['top']['z']
+        z_btm = end_points['btm']['z']
+        return -(z_top - z_btm) / (x_top - x_btm)
+
+    elif method == 'least_squares':
+        profile_mask = [True if end_points['top']['x'] < pts < end_points['btm']['x'] else False for pts in x ]
+        profile_x = np.array(x)[profile_mask].tolist()
+        profile_z = np.array(z)[profile_mask].tolist()
+        slope, _, _, _, _ = stats.linregress(profile_x, profile_z)
+        return -slope
+
+
 def distance_to_intersection(lat,lon,orientation,line_strings):
     """
     Returns the distance at whjch a line drawn from a lat/lon at an orientation intersects a line stinrg