Merge pull request #170 from eEcoLiDAR/development

Development

CHANGELOG.md

@@ -6,6 +6,17 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## Unreleased
 
+## 0.4.2 - 2020-09-18
+## Added:
+- faster implementation of spatial selection
+
+## Changed:
+- all elements in WKT and shapefiles containing multiple polygons are considered (not only the first element)
+- the validity conditions for multipolygons have been relaxed: valid adjacent polygons are now accepted
+
+## Fixed:
+- bug in copy_point_cloud with masks and zero-filled feature-arrays
+
 ## 0.4.1 - 2020-05-20
 ## Added:
 - select_equal filter accepts list of values to compare to the points' attributes

CITATION.cff

@@ -52,7 +52,7 @@ authors:
     family-names: Koma
     given-names: Zsófia
 cff-version: "1.0.3"
-date-released: 2020-05-20
+date-released: 2020-09-18
 doi: "10.5281/zenodo.1219422"
 keywords:
   - "airborne laser scanning"
@@ -62,5 +62,5 @@ keywords:
 license: "Apache-2.0"
 message: "If you use this software, please cite it using these metadata."
 title: "Laserchicken: toolkit for ALS point clouds"
-version: "0.4.1"
+version: "0.4.2"
 ...

laserchicken/_version.txt

@@ -1 +1 @@
-0.4.1
+0.4.2

laserchicken/filter.py

@@ -9,6 +9,7 @@
 from shapely.errors import WKTReadingError
 from shapely.wkt import loads
 from shapely.geometry import box
+from shapely.vectorized import contains
 import numpy as np
 
 from laserchicken.keys import point
@@ -115,9 +116,9 @@ def select_polygon(point_cloud, polygon_string, read_from_file=False, return_mas
     else:
         polygon = _load_polygon(polygon_string)
 
-    if isinstance(polygon, shapely.geometry.polygon.Polygon) and polygon.is_valid:
+    if isinstance(polygon, shapely.geometry.polygon.Polygon):
         points_in = _contains(point_cloud, polygon)
-    elif isinstance(polygon,shapely.geometry.multipolygon.MultiPolygon) and polygon.is_valid:
+    elif isinstance(polygon,shapely.geometry.multipolygon.MultiPolygon):
         points_in = []
         count=1
         for poly in polygon:
@@ -145,17 +146,17 @@ def _read_wkt_file(path):
     with open(path) as f:
         content = f.readlines()
 
-    content = [x.strip() for x in content]
-    return _load_polygon(content[0])
+    content = [_load_polygon(x.strip()) for x in content]
+    geom = shapely.geometry.MultiPolygon(content) if len(content) > 1 else content[0]
+    return geom
 
 
 def _read_shp_file(path):
     shape = shapefile.Reader(path)
-    # first feature of the shapefile
-    feature = shape.shapeRecords()[0]
-    first = feature.shape.__geo_interface__
-    # or shp_geom = shape(first) with PyShp)
-    shp_geom = shapely.geometry.shape(first)
+    features = shape.shapeRecords()
+    shp_geoms = [shapely.geometry.shape(feature.shape.__geo_interface__)
+                 for feature in features]
+    shp_geom = shapely.geometry.MultiPolygon(shp_geoms) if len(shp_geoms) > 1 else shp_geoms[0]
     return shp_geom
 
 
@@ -192,6 +193,9 @@ def _contains(pc, polygon):
     y = pc[point]['y']['data']
     points_in = []
 
+    if not polygon.is_valid:
+        raise ValueError('Invalid polygon in input')
+
     mbr = polygon.envelope
     point_box = box(np.min(x), np.min(y), np.max(x), np.max(y))
 
@@ -204,10 +208,8 @@ def _contains(pc, polygon):
         tree = kd_tree.get_kdtree_for_pc(pc)
         indices = np.sort(tree.query_ball_point(x=p, r=rad))
 
-        point_id = 0
-        for i in indices:            
-            if polygon.contains(Point(x[i], y[i])):
-                points_in.append(i)
-                point_id += 1
+        if len(indices) > 0:
+            mask = contains(polygon, x[indices], y[indices])
+            points_in.extend(indices[mask])
 
     return points_in

laserchicken/test_filter.py

@@ -220,13 +220,6 @@ def test_points_in_polygon_wkt_MultiLine():
         with pytest.raises(ValueError):
             select_polygon(pc_in, "MULTILINESTRING((3 4,10 50,20 25),(-5 -8,-10 -8,-15 -4))")
 
-    @staticmethod
-    def test_points_in_polygon_wkt_MultiPolygon():
-        pc_in = load("testdata/AHN2.las")
-        with pytest.raises(ValueError):
-            select_polygon(pc_in,
-                           "MULTIPOLYGON(((1 1,5 1,5 5,1 5,1 1),(2 2, 3 2, 3 3, 2 3,2 2)),((3 3,6 2,6 4,3 3)))")
-
     @staticmethod
     def test_points_in_polygon_wkt_Collection():
         pc_in = load("testdata/AHN2.las")
@@ -330,12 +323,6 @@ def test_points_in_polygon_wkt_MultiLine():
         with pytest.raises(ValueError):
             select_polygon(pc_in, "testdata/ahn2_geometries_wkt/multiline.wkt", read_from_file=True)
 
-    @staticmethod
-    def test_points_in_polygon_wkt_MultiPolygon():
-        pc_in = load("testdata/AHN2.las")
-        with pytest.raises(ValueError):
-            select_polygon(pc_in, "testdata/ahn2_geometries_wkt/multipolygon.wkt", read_from_file=True)
-
     @staticmethod
     def test_points_in_polygon_wkt_Collection():
         pc_in = load("testdata/AHN2.las")

laserchicken/utils.py

@@ -122,7 +122,8 @@ def copy_point_cloud(source_point_cloud, array_mask=None):
             new_value = copy_point_cloud(value, array_mask)
         elif isinstance(value, np.ndarray):
             if array_mask is not None:
-                new_value = value[array_mask] if any(value) else np.copy(value)
+                new_value = (value[array_mask]
+                             if value.size > 0 else np.copy(value))
             else:
                 new_value = np.copy(value)
         else: