diff --git a/opendm/dem/ground_rectification/io/las_io.py b/opendm/dem/ground_rectification/io/las_io.py
index ca832e9c..918486d1 100755
--- a/opendm/dem/ground_rectification/io/las_io.py
+++ b/opendm/dem/ground_rectification/io/las_io.py
@@ -12,72 +12,40 @@ def read_cloud(point_cloud_path):
 
     # Open point cloud and read its properties using pdal
     pipeline = pdal.Pipeline('[{"type":"readers.las","filename":"%s"}]' % point_cloud_path)
-    cnt = pipeline.execute()
+    pipeline.execute()
 
-    log.ODM_INFO("pdal arrays: %s" % pipeline.arrays)
+    metadata = pipeline.metadata
+    arrays = pipeline.arrays
 
-    dimensions = pipeline.schema['schema']['dimensions']
-    #log.ODM_INFO("Type: %s" % type(pipeline.schema))
-    log.ODM_INFO("Dimensions: %s" % dimensions)
-
-    # The x column index is the index of the object with the name 'X'
-    x_index = next((index for (index, d) in enumerate(dimensions) if d['name'] == 'X'), None)
-    y_index = next((index for (index, d) in enumerate(dimensions) if d['name'] == 'Y'), None)
-    z_index = next((index for (index, d) in enumerate(dimensions) if d['name'] == 'Z'), None)
-    classification_index = next((index for (index, d) in enumerate(dimensions) if d['name'] == 'Classification'), None)
-    red_index = next((index for (index, d) in enumerate(dimensions) if d['name'] == 'Red'), None)
-    green_index = next((index for (index, d) in enumerate(dimensions) if d['name'] == 'Green'), None)
-    blue_index = next((index for (index, d) in enumerate(dimensions) if d['name'] == 'Blue'), None)
-
-    # Log indices
-    log.ODM_INFO("x_index: %s" % x_index)
-    log.ODM_INFO("y_index: %s" % y_index)
-    log.ODM_INFO("z_index: %s" % z_index)
-    log.ODM_INFO("classification_index: %s" % classification_index)
-    log.ODM_INFO("red_index: %s" % red_index)
-    log.ODM_INFO("green_index: %s" % green_index)
-    log.ODM_INFO("blue_index: %s" % blue_index)
-
-    pts = pipeline.arrays[0]
-    log.ODM_INFO("pts: %s" % pts)
-
-    x = (pt[x_index] for pt in pts)
-    y = (pt[y_index] for pt in pts)
-    z = (pt[z_index] for pt in pts)
-    classification = (pt[classification_index] for pt in pts)
-    red = (pt[red_index] for pt in pts)
-    green = (pt[green_index] for pt in pts)
-    blue = (pt[blue_index] for pt in pts)
+    # Extract point coordinates, classification, and RGB values
+    x = arrays[0]["X"]
+    y = arrays[0]["Y"]
+    z = arrays[0]["Z"]
+    classification = arrays[0]["Classification"].astype(np.uint8)
+    red = arrays[0]["Red"]
+    green = arrays[0]["Green"]
+    blue = arrays[0]["Blue"]
 
+    # Create PointCloud object
     cloud = PointCloud.with_dimensions(x, y, z, classification, red, green, blue)
 
     # Return the result
-    return pipeline.metadata, cloud
+    return metadata, cloud
 
-def write_cloud(header, point_cloud, output_point_cloud_path, write_extra_dimensions=False):
-    # Open output file
-    output_las_file = laspy.LasData(header)
+def write_cloud(metadata, point_cloud, output_point_cloud_path, write_extra_dimensions=False):
 
-    if write_extra_dimensions:
-        extra_dims = [laspy.ExtraBytesParams(name=name, type=dimension.get_las_type(), description="Dimension added by Ground Extend") for name, dimension in point_cloud.extra_dimensions_metadata.items()]
-        output_las_file.add_extra_dims(extra_dims)
-        # Assign dimension values
-        for dimension_name, values in point_cloud.extra_dimensions.items():
-            setattr(output_las_file, dimension_name, values)
+    # Create PDAL pipeline to write point cloud
+    pipeline = pdal.Pipeline('[{"type": "writers.las","filename": "%s","compression": "laszip","extra_dims": %s}]' %
+                             (output_point_cloud_path, str(write_extra_dimensions).lower()))
 
     # Adapt points to scale and offset
     [x, y] = np.hsplit(point_cloud.xy, 2)
-    output_las_file.x = x.ravel()
-    output_las_file.y = y.ravel()
-    output_las_file.z = point_cloud.z
+    z = point_cloud.z
 
     # Set color
     [red, green, blue] = np.hsplit(point_cloud.rgb, 3)
-    output_las_file.red = red.ravel()
-    output_las_file.green = green.ravel()
-    output_las_file.blue = blue.ravel()
-
     # Set classification
-    output_las_file.classification = point_cloud.classification.astype(np.uint8)
+    classification = point_cloud.classification.astype(np.uint8)
 
-    output_las_file.write(output_point_cloud_path)
\ No newline at end of file
+    # Write point cloud with PDAL
+    pipeline.execute(np.column_stack((x, y, z, red, green, blue, classification)))