Add EPT reader support to read_lidar function

pyforestscan/calculate.py

@@ -148,16 +148,20 @@ def calculate_pai(pad, min_height=1, max_height=None):
 
 def calculate_fhd(voxel_returns):
     """
-    Calculate the Frequency Histogram Diversity (FHD) for a given set of voxel returns.
+    Calculate the Foliage Height Diversity (FHD) for a given set of voxel returns.
 
-    This function computes the Frequency Histogram Diversity by calculating the entropy
-    of the voxel return proportions along the z-axis.
+    This function computes Foliage Height Diversity by calculating the entropy
+    of the voxel return proportions along the z-axis, which represents vertical structure
+    in the canopy.
 
     :param voxel_returns:
-        A numpy array of shape (x, y, z) representing voxel returns.
+        A numpy array of shape (x, y, z) representing voxel returns, where x and y are spatial
+        dimensions, and z represents height bins (or layers along the vertical axis).
 
     :return:
-        A numpy array of shape (x, y"""
+        A numpy array of shape (x, y) representing the FHD values for each (x, y) location.
+        Areas with no voxel returns will have NaN values.
+    """
     sum_counts = np.sum(voxel_returns, axis=2, keepdims=True)
 
     with np.errstate(divide='ignore', invalid='ignore'):

pyforestscan/handlers.py

@@ -158,18 +158,19 @@ def validate_crs(crs_list):
     return True
 
 
-def read_lidar(input_file, srs, thin_radius=None, hag=False, hag_dtm=False, dtm=None, crop_poly=False, poly=None):
+def read_lidar(input_file, srs, bounds=None, thin_radius=None, hag=False, hag_dtm=False, dtm=None, crop_poly=False, poly=None):
     """
     Reads and processes a LiDAR point cloud file using PDAL based on specified options.
 
-    :param srs:
     :param input_file: str, The path to the input LiDAR file. Supported formats are .las, .laz, .copc, and .copc.laz.
+    :param srs: str, The Spatial Reference System (SRS) of the point cloud.
+    :param bounds:  Bounds within which to crop the data. Must be of the form: ([xmin, xmax], [ymin, ymax], [zmin, zmax])
     :param thin_radius: float, optional, The radius for thinning the point cloud. Must be a positive number.
     :param hag: bool, optional, If True, calculate Height Above Ground (HAG) using Delaunay triangulation.
     :param hag_dtm: bool, optional, If True, calculate Height Above Ground (HAG) using a DTM file.
     :param dtm: str, optional, The path to the DTM file used when hag_dtm is True. Must be a .tif file.
     :param crop_poly: bool, optional, If True, crop the point cloud using the polygon defined in the poly file.
-    :param poly: str, optional, The path to the polygon file used for cropping.
+    :param poly: str, optional, The path to the polygon file used for cropping OR the WKT of the Polygon geometry.
 
     :return: numpy.ndarray, The processed point cloud data or None if no data is retrieved.
 
@@ -183,14 +184,19 @@ def read_lidar(input_file, srs, thin_radius=None, hag=False, hag_dtm=False, dtm=
 
     las_extensions = ('.las', '.laz')
     copc_extensions = ('.copc', '.copc.laz')
+    ept_file = ('ept.json')
 
     file_lower = input_file.lower()
     if file_lower.endswith(las_extensions):
         reader = 'readers.las'
     elif file_lower.endswith(copc_extensions):
         reader = 'readers.copc'
+    elif file_lower.endswith(ept_file):
+        reader = 'readers.ept'
     else:
-        raise ValueError("The input file must be a .las, .laz, .copc, or .copc.laz file.")
+        raise ValueError(
+            "The input file must be a .las, .laz, .copc, .copc.laz file, or an ept.json file."
+        )
 
     if hag and hag_dtm:
         raise ValueError("Cannot use both 'hag' and 'hag_dtm' options at the same time.")
@@ -199,10 +205,14 @@ def read_lidar(input_file, srs, thin_radius=None, hag=False, hag_dtm=False, dtm=
     crs_list = []
 
     if crop_poly:
-        if not poly or not os.path.isfile(poly):
-            raise FileNotFoundError(f"No such polygon file: '{poly}'")
-        polygon_wkt, crs_vector = load_polygon_from_file(poly)
-        crs_list.append(crs_vector)
+        if not poly:
+            raise ValueError(f"Must provide a polygon or polygon wkt if cropping to a polygon.")
+        if poly.strip().startswith(('POLYGON', 'MULTIPOLYGON')):
+            polygon_wkt = poly
+        else:
+            if not poly or not os.path.isfile(poly):
+                raise FileNotFoundError(f"No such polygon file: '{poly}'")
+            polygon_wkt, crs_vector = load_polygon_from_file(poly)
         pipeline_stages.append(_crop_polygon(polygon_wkt))
 
     if thin_radius is not None:
@@ -224,39 +234,16 @@ def read_lidar(input_file, srs, thin_radius=None, hag=False, hag_dtm=False, dtm=
         crs_list.append(crs_raster)
         pipeline_stages.append(_hag_raster(dtm))
 
-    metadata_pipeline = pdal.Pipeline(
-        json.dumps({
-            "pipeline": [
-                {
-                    "type": reader,
-                    "spatialreference": srs, #todo: make srs optional
-                    "filename": input_file
-                },
-                {
-                    "type": "filters.info"
-                }
-            ]
-        })
-    )
-    metadata_pipeline.execute()
-
-    try:
-        crs_pointcloud = metadata_pipeline.metadata['metadata']['readers.las']['spatialreference']
-    except KeyError:
-        try:
-            crs_pointcloud = metadata_pipeline.metadata['metadata']['readers.copc']['spatialreference']
-        except KeyError:
-            raise ValueError("Unable to retrieve spatial reference from the point cloud metadata.")
-    crs_list.append(crs_pointcloud)
-
-    validate_crs(crs_list)
-
+    base_pipeline = {
+        "type": reader,
+        "spatialreference": srs,
+        "filename": input_file
+    }
+    if bounds:
+        base_pipeline["bounds"] = f"{bounds}"
     main_pipeline_json = {
         "pipeline": [
-            {
-                "type": reader,
-                "filename": input_file
-            }
+            base_pipeline
         ] + pipeline_stages
     }
 
@@ -296,12 +283,14 @@ def write_las(arrays, output_file, srs=None, compress=True):
     if srs:
         pipeline_steps.append({
             "type": "filters.reprojection",
+            "in_srs": srs,
             "out_srs": srs
         })
 
     pipeline_steps.append({
         "type": output_format,
         "filename": output_file,
+        "minor_version": "4",
         "extra_dims": "all"
     })
 

setup.py

@@ -5,7 +5,7 @@
 
 setuptools.setup(
     name="pyforestscan",
-    version="0.1.6",
+    version="0.1.7",
     author="Joseph Emile Honour Percival",
     author_email="ipercival@gmail.com",
     description="Analyzing forest structure using aerial LiDAR data",

tests/test_handlers.py

@@ -221,7 +221,7 @@ def test_read_lidar_missing_dtm_file():
 def test_read_lidar_unsupported_input_extension():
     input_file = get_test_data_path("input.txt")
     srs = "EPSG:4326"
-    with pytest.raises(ValueError, match="The input file must be a .las, .laz, .copc, or .copc.laz file."):
+    with pytest.raises(ValueError, match="The input file must be a .las, .laz, .copc, .copc.laz file, or an ept.json file."):
         read_lidar(
             input_file=input_file,
             srs=srs