Get point selection in line with what FUSION does

src/silvimetric/commands/shatter.py

@@ -42,7 +42,7 @@ def arrange(points: pd.DataFrame, leaf, attrs: list[str]):
     if points.size == 0:
         return None
 
-    points = points.loc[points.Y > leaf.bounds.miny]
+    points = points.loc[points.Y < leaf.bounds.maxy]
     points = points.loc[points.X < leaf.bounds.maxx, [*attrs, 'xi', 'yi']]
 
     points.loc[:, 'xi'] = da.floor(points.xi)

src/silvimetric/resources/data.py

@@ -61,14 +61,9 @@ def make_pipeline(self) -> pdal.Pipeline:
         class_zero = pdal.Filter.assign(value="Classification = 0")
         rn = pdal.Filter.assign(value="ReturnNumber = 1 WHERE ReturnNumber < 1")
         nor = pdal.Filter.assign(value="NumberOfReturns = 1 WHERE NumberOfReturns < 1")
-        ferry = pdal.Filter.ferry(dimensions="X=>xi, Y=>yi")
-
-        # assign_x = pdal.Filter.assign(
-        #     value=f"xi = (X - {self.root.minx}) / {resolution}")
-        # assign_y = pdal.Filter.assign(
-        #     value=f"yi = ({self.root.maxy} - Y) / {resolution}")
+        # ferry = pdal.Filter.ferry(dimensions="X=>xi, Y=>yi")
 
-        return reader | class_zero | rn | nor #| ferry | assign_x | assign_y #| smrf | hag
+        return reader | class_zero | rn | nor
 
     def get_pipeline(self) -> pdal.Pipeline:
         """Fetch the pipeline for the instance
@@ -122,7 +117,7 @@ def get_pipeline(self) -> pdal.Pipeline:
         # Add xi and yi – only need this for PDAL < 2.6
         ferry = pdal.Filter.ferry(dimensions="X=>xi, Y=>yi")
         assign_x = pdal.Filter.assign(value=f"xi = (X - {self.storageconfig.root.minx}) / {resolution}")
-        assign_y = pdal.Filter.assign(value=f"yi = ({self.storageconfig.root.maxy} - Y) / {resolution}")
+        assign_y = pdal.Filter.assign(value=f"yi = (Y - {self.storageconfig.root.miny}) / {resolution}")
         stages.append(ferry)
         stages.append(assign_x)
         stages.append(assign_y)

src/silvimetric/resources/extents.py

@@ -41,11 +41,11 @@ def __init__(self,
 
         self.x1 = math.floor((minx - self.root.minx) / resolution)
         """Minimum X index"""
-        self.y1 = math.floor((self.root.maxy - maxy) / resolution)
+        self.y1 = math.floor((miny - self.root.miny) / resolution)
         """Minimum Y index"""
         self.x2 = math.ceil((maxx - self.root.minx) / resolution)
         """Maximum X index"""
-        self.y2 = math.ceil((self.root.maxy - miny) / resolution)
+        self.y2 = math.ceil((maxy - self.root.miny) / resolution)
         """Maximum Y index"""
         self.domain: IndexDomainList = ((self.x1, self.x2), (self.y1, self.y2))
         """Minimum bounding rectangle of this Extents"""
@@ -111,8 +111,9 @@ def chunk(self, data: Data, res_threshold=100, pc_threshold=600000,
         # make bounds in scale with the desired resolution
         minx = bminx + (self.x1 * self.resolution)
         maxx = bminx + (self.x2 * self.resolution)
-        miny = bmaxy - (self.y2 * self.resolution)
-        maxy = bmaxy - (self.y1 * self.resolution)
+
+        miny = bminy + (self.y1 * self.resolution)
+        maxy = bminy + (self.y2 * self.resolution)
 
         chunk = Extents(Bounds(minx, miny, maxx, maxy), self.resolution, r)
 
@@ -151,7 +152,7 @@ def split(self):
         y_adjusted = math.floor((maxy - miny) / 2 / self.resolution)
 
         midx = minx + (x_adjusted * self.resolution)
-        midy = maxy - (y_adjusted * self.resolution)
+        midy = miny + (y_adjusted * self.resolution)
 
         return [
             Extents(Bounds(minx, miny, midx, midy), self.resolution, self.root), #lower left
@@ -237,11 +238,17 @@ def get_leaf_children(self, tile_size):
             ], dtype=np.float64)
         dx = (res * local_xs) + self.root.minx
 
+        # local_ys = np.array([
+        #         [min(y+ynum, self.y2), y]
+        #         for y in range(self.y1, self.y2, int(ynum))
+        #     ], dtype=np.float64)
+        # dy = (res * local_ys) + self.root.miny
+
         local_ys = np.array([
-                [min(y+ynum, self.y2), y]
+                [y, min(y+ynum, self.y2)]
                 for y in range(self.y1, self.y2, int(ynum))
             ], dtype=np.float64)
-        dy = self.root.maxy - (res * local_ys)
+        dy = (res * local_ys) + self.root.miny
 
         coords_list = np.array([[*x,*y] for x in dx for y in dy],dtype=np.float64)
         yield from [

tests/create_test_data.py

@@ -1,6 +1,6 @@
 import pdal
 import numpy as np
-from math import sqrt, ceil
+from math import sqrt, floor
 
 
 # Create a pointcloud dataset that follows expected paramters for data testing
@@ -28,10 +28,10 @@
 diff_maker = 0
 
 x_pos = np.arange(minx, maxx, interval, dtype=np.float32)
-y_pos = np.arange(maxy, miny, -1*interval, dtype=np.float32)
+y_pos = np.arange(miny, maxy, interval, dtype=np.float32)
 # positions = pos[np.where(pos % cell_size != 0)]
 
-alg = lambda y: ceil(y/cell_size) + diff_maker
+alg = lambda y: floor(y/cell_size) + diff_maker
 
 data = np.array([(x, y, alg(y), alg(y), alg(y), alg(y))
                   for x in x_pos for y in y_pos],

tests/test_cli.py

@@ -1,5 +1,4 @@
 import numpy as np
-import pytest
 
 from silvimetric.cli import cli
 from silvimetric.commands import shatter, info
@@ -13,7 +12,7 @@ def test_cli_init(self, tdb_filepath, runner, bounds):
                 catch_exceptions=False)
         assert res.exit_code == 0
 
-    def test_cli_shatter(self, runner, storage, maxy, date, tdb_filepath,
+    def test_cli_shatter(self, runner, storage, miny, date, tdb_filepath,
             copc_filepath):
 
         res = runner.invoke(cli.cli, args=["-d", tdb_filepath,
@@ -31,7 +30,7 @@ def test_cli_shatter(self, runner, storage, maxy, date, tdb_filepath,
             for xi in range(x):
                 for yi in range(y):
                     a[xi, yi]['Z'].size == 1
-                    assert bool(np.all(a[xi, yi]['Z'][0] == ((maxy/storage.config.resolution)-yi)))
+                    assert bool(np.all( a[xi, yi]['Z'][0] == ((miny/storage.config.resolution) + yi) ))
 
     def test_cli_scan(self, tdb_filepath, runner, copc_filepath, storage):
         res = runner.invoke(cli.cli, args=['-d', tdb_filepath, '--scheduler', 'single-threaded', 'scan', copc_filepath])

tests/test_extract.py

@@ -27,7 +27,7 @@ def tif_test(extract_config):
         assert raster.RasterYSize == ysize
 
         r = raster.ReadAsArray()
-        assert all([ r[y,x] == ((maxy/resolution)-(y+yimin))  for y in range(ysize) for x in range(xsize)])
+        assert all([ r[y,x] == ((miny/resolution)+(y+yimin))  for y in range(ysize) for x in range(xsize)])
 
 class Test_Extract(object):
 

tests/test_shatter.py

@@ -24,7 +24,7 @@ def write(x,y,val, s:Storage, attrs, dims, metrics):
 
 class Test_Shatter(object):
 
-    def test_shatter(self, shatter_config, storage: Storage, maxy):
+    def test_shatter(self, shatter_config, storage: Storage, miny):
         shatter(shatter_config)
         with storage.open('r') as a:
             y = a[:,0]['Z'].shape[0]
@@ -34,7 +34,7 @@ def test_shatter(self, shatter_config, storage: Storage, maxy):
             for xi in range(x):
                 for yi in range(y):
                     a[xi, yi]['Z'].size == 1
-                    assert bool(np.all(a[xi, yi]['Z'][0] == ((maxy/storage.config.resolution)-yi)))
+                    assert bool(np.all( a[xi, yi]['Z'][0] == ((miny/storage.config.resolution) + yi) ))
 
         # change attributes to make it a new run
         shatter_config.name = uuid.uuid4()
@@ -53,7 +53,7 @@ def test_shatter(self, shatter_config, storage: Storage, maxy):
                 for yi in range(y):
                     a[xi, yi]['Z'].size == 2
                     assert bool(np.all(a[xi, yi]['Z'][1] == a[xi,yi]['Z'][0]))
-                    assert bool(np.all(a[xi, yi]['Z'][1] == ((maxy/storage.config.resolution)-yi)))
+                    assert bool(np.all(a[xi, yi]['Z'][1] == ((miny/storage.config.resolution) + yi)))
 
         m = info(storage.config.tdb_dir)
         assert len(m['history']) == 2
@@ -114,7 +114,7 @@ def test_sub_bounds(self, sh_cfg, test_point_count, request):
 
     def test_partial_overlap(self, partial_shatter_config, test_point_count):
         pc = shatter(partial_shatter_config)
-        assert pc == 23100
+        assert pc == 23716
 
     @pytest.mark.skipif(
         os.environ.get('AWS_SECRET_ACCESS_KEY') is None or
@@ -125,7 +125,7 @@ def test_remote_creation(self, s3_shatter_config, s3_storage):
         # need processes scheduler to accurately test bug fix
         dask.config.set(scheduler="processes")
         resolution = s3_storage.config.resolution
-        maxy = s3_storage.config.root.maxy
+        miny = s3_storage.config.root.miny
         shatter(s3_shatter_config)
         with s3_storage.open('r') as a:
             y = a[:,0]['Z'].shape[0]
@@ -135,4 +135,4 @@ def test_remote_creation(self, s3_shatter_config, s3_storage):
             for xi in range(x):
                 for yi in range(y):
                     a[xi, yi]['Z'].size == 1
-                    assert bool(np.all(a[xi, yi]['Z'][0] == ((maxy/resolution)-yi)))
+                    assert bool(np.all( a[xi, yi]['Z'][0] == ((miny/resolution) + yi) ))