diff --git a/src/data/profile_features.py b/src/data/profile_features.py
index 5ed2018..5825b43 100644
--- a/src/data/profile_features.py
+++ b/src/data/profile_features.py
@@ -19,10 +19,12 @@ def shapes_from_shp(shp_file):
     """
     shapes = []
     ids = []
-    for feat in fiona.open(shp_file,'r'):
+    properties = []
+    for feat in fiona.open(shp_file, 'r'):
         shapes.append(shape(feat['geometry']))
         ids.append(feat['id'])
-    return shapes, ids
+        properties.append(feat['properties'])
+    return shapes, ids, properties
 
 
 def convert_coord_systems(g1, in_coord_system='EPSG:4326', out_coord_system='EPSG:28356'):
@@ -98,12 +100,12 @@ def get_slope(x, z, top_elevation, btm_elevation, method='end_points'):
         return -slope
 
 
-def distance_to_intersection(lat, lon, orientation, line_strings):
+def distance_to_intersection(lat, lon, landward_orientation, beach, line_strings, line_properties):
     """
     Returns the distance at whjch a line drawn from a lat/lon at an orientation intersects a line stinrg
     :param lat:
     :param lon:
-    :param orientation: Angle, clockwise positive from true north in degrees, of the tangent to the shoreline facing
+    :param landward_orientation: Angle, anticlockwise positive from east in degrees, towards the land
     towards the
     land.
     :param line_string:
@@ -113,15 +115,26 @@ def distance_to_intersection(lat, lon, orientation, line_strings):
     start_point = convert_coord_systems(start_point)
 
     distance = 1000  # m look up to 1000m for an intersection
-    new_point = Point(start_point.coords.xy[0] + distance * np.cos(np.deg2rad(orientation)),
-                      start_point.coords.xy[1] + distance * np.sin(np.deg2rad(orientation)))
-    profile_line = LineString([start_point, new_point])
-
-    # Check whether profile_line intersects with any lines in line_string
+    landward_point = Point(start_point.coords.xy[0] + distance * np.cos(np.deg2rad(landward_orientation)),
+                           start_point.coords.xy[1] + distance * np.sin(np.deg2rad(landward_orientation)))
+    landward_line = LineString([start_point, landward_point])
+    seaward_point = Point(start_point.coords.xy[0] - distance * np.cos(np.deg2rad(landward_orientation)),
+                          start_point.coords.xy[1] - distance * np.sin(np.deg2rad(landward_orientation)))
+    seaward_line = LineString([start_point, seaward_point])
+
+    # Look at relevant line_strings which have the same beach property in order to reduce computation time
+    line_strings = [s for s, p in zip(line_strings, line_properties) if p['beach'] == beach]
+
+    # Check whether profile_line intersects with any lines in line_string. If intersection point is landwards,
+    # consider this negative, otherwise seawards is positive.
     for line_string in line_strings:
-        intersection_points = profile_line.intersection(line_string)
-        if not intersection_points.is_empty:
-            return intersection_points.distance(start_point)
+        land_intersect_points = landward_line.intersection(line_string)
+        if not land_intersect_points.is_empty:
+            return -land_intersect_points.distance(start_point)
+
+        sea_intersect_points = seaward_line.intersection(line_string)
+        if not sea_intersect_points.is_empty:
+            return sea_intersect_points.distance(start_point)
 
     # If no intersections are found, return nothing.
     return None