diff --git a/outputs_2017088_Survey2.py b/outputs_2017088_Survey2.py
index ef1eb2f..29f3dd4 100644
--- a/outputs_2017088_Survey2.py
+++ b/outputs_2017088_Survey2.py
@@ -96,21 +96,57 @@ def make_raster(las, output, lastools_loc, keep_only_ground=False, step=0.2):
         return None
 
 
-def extract_pts(las, in_csv, tmp_csv, lastools_loc, keep_only_ground=True):
-    # output is the full path and filename (inc extension) to put in
-    path_2_lascontrol=lastools_loc+"\\bin\\lascontrol"
-
-    if keep_only_ground==True:
-        command="%s -i %s -cp %s -keep_class 2 -parse sxyz -cp_out %s" % (path_2_lascontrol, las, in_csv, tmp_csv)
-    else:
-        command="%s -i %s -cp %s -parse sxyz -cp_out %s" % (path_2_lascontrol, las, in_csv, tmp_csv)
-
-    returncode,output = check_output(command, False)
-    if returncode!= 0:
-        print("Error. lascontrol failed on %s" % las.split('//')[-1].split('.')[0])
-        print(output)
-    else:
-        return None
+def extract_pts(las_in, cp_in, survey_date, keep_only_ground=True):
+    """Extract elevations from a las surface based on x and y coordinates.
+
+    Requires lastools in system path.
+
+    Args:
+        las_in:            input point cloud (las)
+        cp_in:             point coordinates with columns: id, x, y, z (csv)
+        survey_date:       survey date string, e.g. '19700101'
+        keep_only_ground:  only keep points classified as 'ground' (boolean)
+
+    Returns:
+        Dataframe containing input coordinates with extracted elevations
+    """
+
+    cmd = ['lascontrol', '-i', las_in, '-cp', cp_in, '-parse', 'sxyz']
+
+    if keep_only_ground == True:
+        cmd += ['-keep_class', '2']
+
+    # Call lastools
+    process = subprocess.Popen(
+        cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
+    stdout, stderr = process.communicate()
+    errcode = process.returncode
+
+    # Handle errors, if detected
+    if errcode != 0:
+        print("Error. lascontrol failed on {}".format(
+            os.path.basename(las_in)))
+        print(stderr.decode())
+
+    # Load result into pandas dataframe
+    df = pd.read_csv(io.BytesIO(stdout))
+
+    # Tidy up dataframe
+    df = df.drop(columns=['diff'])
+    df['lidar_z'] = pd.to_numeric(df['lidar_z'], errors='coerce')
+    df['Beach'] = beach
+    df = df[[
+        'Beach', 'ProfileNum', 'Easting', 'Northing', 'Chainage', 'lidar_z'
+    ]]
+
+    # Rename columns
+    new_names = {
+        'ProfileNum': 'Profile',
+        'lidar_z': 'Elevation_{}'.format(survey_date),
+    }
+    df = df.rename(columns=new_names)
+
+    return df