diff --git a/src/data/profile_features.py b/src/data/profile_features.py
index 3aa0c66..742ee2b 100644
--- a/src/data/profile_features.py
+++ b/src/data/profile_features.py
@@ -165,13 +165,47 @@ def get_sites_dune_crest_toe():
     :return:
     """
 
-    # Get site information
+    # Get site information. Base our profile features on each site
     data_folder = './data/interim'
     df_sites = pd.read_csv(os.path.join(data_folder, 'sites.csv'), index_col=[0])
+    df_profiles = pd.read_csv(os.path.join(data_folder, 'profiles.csv'), index_col=[0, 1, 2])
+    df_profile_features = df_sites
+
+    features = {
+        'dune_crest':
+            {
+                'file': './data/raw/profile_features/dune_crests.shp'
+            },
+        'dune_toe':
+            {
+                'file': './data/raw/profile_features/dune_toes.shp'
+            },
+    }
 
     # Import our dune crest and toes
-    for f in ['./data/raw/profile_features/dune_crests.shp']:
-        shapes, _ = shapes_from_shp(f)
+    for feat in features.keys():
+        shapes, _, properties = shapes_from_shp(features[feat]['file'])
         shapes = [convert_coord_systems(x) for x in shapes]
 
-    # Iterate through each site
+        # Figure out the x coordinates of our crest and toes, by looking at where our beach sections intersect our
+        # shape files.
+        col_name = '{}_x'.format(feat)
+        df_profile_features[col_name] = df_profile_features['profile_x_lat_lon'] + \
+                                        df_profile_features.apply(lambda row:
+                                                                  distance_to_intersection(
+                                                                      row['lat'], row['lon'], row['orientation'],
+                                                                      row['beach'], shapes, properties),
+                                                                  axis=1)
+        # Get the elevations of the crest and toe
+        col_name = '{}_z'.format(feat)
+        df_profile_features[col_name] = df_profile_features.apply(lambda row:
+                                                                  beach_profile_elevation(
+                                                                      row['{}_x'.format(feat)],
+                                                                      df_profiles,
+                                                                      'prestorm',
+                                                                      row.name),
+                                                                  axis=1)
+
+    df_profile_features = df_profile_features.drop(columns=['beach', 'lat','lon','orientation'])
+    df_profile_features.to_csv('./data/interim/profile_features.csv')
+