Update bounds.py align_to_cell_lines() to match FUSION alignment (#112)

* Update bounds.py align_to_cell_lines() to match FUSION alignment

Modified code to align cell centers with multiple of resolution. It might be better to add an option to the database for pixel-is-point or pixel-is-cell and rename the new version of align_to_cell_lines() to align_to_cell_centers(). There are just a few places where the function is called that would need to add code to use the db option.

* Update bounds.py to include commented original code

* adding option for cell alignment to DB

* making adjust_alignment self-contained

* Updating tests and CLI to use alignment

* Setting default alignment to pixel_is_cell to pass tests

* added alignment fixture

* fixing duplicate parameter

* Fixed default alignment in CLI

* added alignment to s3_storage_config fixture

* Changing pixel_is_point and pixel_is_cell to pixelispoint and pixelisarea. Adding shatter test for pixelispoint case

* Correcting alignment in tests

* testing alignment

* Changing test_command_point() to modify storage.config with pixelispoint alignment

* Removing test for pixelispoint...more complicated than I thought!

* Adding CRS to bounds for pipeline reader

* Adding CRS to bounds when pipeline is created

* backed out CRS addition to bounds

* Updated to use AlignToCorner and AlignToCenter value for pixel alignment

* Parametrized Tests for alignment addition (#115)

* tests working after addition of parametrized fixture for alignment

* remove debugging line

* Fixing GDAL data type mapping and adding check for what to make nanvalue if dtype isn't float

* Grid metrics (#121)

* adding grid_metrics and corresponding filters

* adding lfs, grid_metrics, and a test for fusion comparisons

* adding fusion test case file

* adding fusion comparison rasters to test data

* wip, debugging differences in intensity tif output

* Add linting to keep consistency (#120)

* linting the project according to ruff, added config to pyproject.toml

* wip, getting tests outside of fusion working again

* fixing tests, skipping fusion test

* fixing linting errors

* remove debug noise, fix filter param name in scan

* fixed lmoments and cv metrics, now matching FUSION

* changing mode to mimic how FUSION is doing it

* added rest of fusion metrics

* added methods for counts and covers so elevation key can be applied correctly

* removing alignment from tests it's unnecessary in, fixing tests that relied on checking data

* adding docstrings, updating cli docs

* attempt to add fusion test back with lfs in gh action workflow

---------

Co-authored-by: Kyle Mann <kyle@hobu.co>

Author: bmcgaughey1

.gitattributes

@@ -0,0 +1 @@
+tests/data/NoCAL_PlumasNF_B2_2018_TestingData_FUSIONNormalized.copc.laz filter=lfs diff=lfs merge=lfs -text

.github/workflows/main.yml

@@ -30,6 +30,8 @@ jobs:
     steps:
       - name: Checkout
         uses: actions/checkout@v4
+        with:
+          lfs: 'true'
       - uses: mamba-org/setup-micromamba@v2
         with:
           init-shell: bash

docs/create.sh

@@ -1,4 +1,5 @@
 #!/bin/bash
 
 # This assumes you've already done pip install -r requirements.txt in docs dir
+conda run -n silvimetric jupyter-book config sphinx source
 python -m sphinx -T -b html -d _build/doctrees -D language=en source html

docs/source/_toc.yml

@@ -7,7 +7,6 @@ parts:
 - caption: Getting Started
   chapters:
   - file: about
-  - file: quickstart
 - caption: CLI
   chapters:
   - file: cli/extract

docs/source/tutorial.md

@@ -126,13 +126,17 @@ Usage:
 ```console
 Usage: silvimetric initialize [OPTIONS]
 
+  Initialize silvimetrics DATABASE
+
 Options:
---bounds BOUNDS         Root bounds that encapsulates all data  [required]
---crs CRS               Coordinate system of data  [required]
--a, --attributes ATTRS  List of attributes to include in Database
--m, --metrics METRICS   List of metrics to include in Database
---resolution FLOAT      Summary pixel resolution
---help                  Show this message and exit.
+  --bounds BOUNDS         Root bounds that encapsulates all data  [required]
+  --crs CRS               Coordinate system of data  [required]
+  -a, --attributes ATTRS  List of attributes to include in Database
+  -m, --metrics METRICS   List of metrics to include in output, eg. '-m
+                          stats,percentiles'
+  --resolution FLOAT      Summary pixel resolution
+  --alignment TEXT        Pixel alignment: 'AlignToCenter' or 'AlignToCorner'
+  --help                  Show this message and exit.
 ```
 
 (udm)=
@@ -212,16 +216,17 @@ Usage:
 ```console
 Usage: silvimetric scan [OPTIONS] POINTCLOUD
 
-Scan point cloud and determine the optimal tile size.
+  Scan point cloud, output information on it, and determine the optimal tile
+  size.
 
 Options:
---resolution FLOAT     Summary pixel resolution
---filter               Remove empty space in computation. Will take extra
-                        time.
---point_count INTEGER  Point count threshold.
---depth INTEGER        Quadtree depth threshold.
---bounds BOUNDS        Bounds to scan.
---help                 Show this message and exit.
+  --resolution FLOAT     Summary pixel resolution
+  --filter_empty         Remove empty space in computation. Will take extra
+                         time.
+  --point_count INTEGER  Point count threshold.
+  --depth INTEGER        Quadtree depth threshold.
+  --bounds BOUNDS        Bounds to scan.
+  --help                 Show this message and exit.
 ```
 
 Example:
@@ -246,14 +251,14 @@ Usage:
 Usage: silvimetric shatter [OPTIONS] POINTCLOUD
 
 Options:
---bounds BOUNDS                 Bounds for data to include in processing
---tilesize INTEGER              Number of cells to include per tile
---report                        Whether or not to write a report of the
-                                process, useful for debugging
---date [%Y-%m-%d|%Y-%m-%dT%H:%M:%SZ]
-                                Date the data was produced.
---dates <DATETIME DATETIME>...  Date range the data was produced during
---help                          Show this message and exit.
+  --bounds BOUNDS                 Bounds for data to include in processing
+  --tilesize INTEGER              Number of cells to include per tile
+  --report                        Whether or not to write a report of the
+                                  process, useful for debugging
+  --date [%Y-%m-%d|%Y-%m-%dT%H:%M:%SZ]
+                                  Date the data was produced.
+  --dates <DATETIME DATETIME>...  Date range the data was produced during
+  --help                          Show this message and exit.
 ```
 
 Example:
@@ -279,12 +284,12 @@ Usage:
 Usage: silvimetric info [OPTIONS]
 
 Options:
---bounds BOUNDS                 Bounds to filter by
---date [%Y-%m-%d|%Y-%m-%dT%H:%M:%SZ]
-                                Select processes with this date
---dates <DATETIME DATETIME>...  Select processes within this date range
---name TEXT                     Select processes with this name
---help                          Show this message and exit.
+  --bounds BOUNDS                 Bounds to filter by
+  --date [%Y-%m-%d|%Y-%m-%dT%H:%M:%SZ]
+                                  Select processes with this date
+  --dates <DATETIME DATETIME>...  Select processes within this date range
+  --name TEXT                     Select processes with this name
+  --help                          Show this message and exit.
 ```
 
 Example:
@@ -359,11 +364,11 @@ Usage: silvimetric extract [OPTIONS]
 Extract silvimetric metrics from DATABASE
 
 Options:
--a, --attributes ATTRS  List of attributes to include output
--m, --metrics METRICS   List of metrics to include in output
---bounds BOUNDS         Bounds for data to include in output
--o, --outdir PATH       Output directory.  [required]
---help                  Show this message and exit.
+  -a, --attributes ATTRS  List of attributes to include output
+  -m, --metrics METRICS   List of metrics to include in output
+  --bounds BOUNDS         Bounds for data to include in output
+  -o, --outdir PATH       Output directory.  [required]
+  --help                  Show this message and exit.
 ```
 
 Example:

environment.yml

@@ -17,5 +17,4 @@ dependencies:
   - python-json-logger
   - dill
   - pandas
-  - lmoments3
-  - typing_extensions
+  - typing_extensions

pyproject.toml

@@ -47,6 +47,8 @@ docstring-code-format = true
 select=[
   "A",
   "B",
+  "D101",
+  "D201",
   "E",
   "F",
   "ISC",

src/silvimetric/cli/cli.py

@@ -255,6 +255,12 @@ def scan_cmd(
 @click.option(
     '--resolution', type=float, default=30.0, help='Summary pixel resolution'
 )
+@click.option(
+    '--alignment',
+    type=str,
+    default='AlignToCenter',
+    help="Pixel alignment: 'AlignToCenter' or 'AlignToCorner'",
+)
 @click.pass_obj
 def initialize_cmd(
     app: ApplicationConfig,
@@ -263,6 +269,7 @@ def initialize_cmd(
     attributes: list[Attribute],
     resolution: float,
     metrics: list[Metric],
+    alignment: str
 ):
     """Initialize silvimetrics DATABASE"""
 
@@ -274,6 +281,7 @@ def initialize_cmd(
         attrs=attributes,
         metrics=metrics,
         resolution=resolution,
+        alignment=alignment
     )
     return initialize.initialize(storageconfig)
 

src/silvimetric/cli/common.py

@@ -1,12 +1,10 @@
 import click
 import pyproj
 import webbrowser
-from contextlib import nullcontext
 
 import dask
 from dask.diagnostics import ProgressBar
 from dask.distributed import Client, LocalCluster
-from distributed.client import _get_global_client as get_client
 
 from ..resources.metrics import (
     l_moments,
@@ -20,6 +18,14 @@
 
 
 class BoundsParamType(click.ParamType):
+    """Click parameter type for the Bounds class.
+
+    Accepts bounds as a stringified JSON object or bbox array. Examples:
+    "([1,101],[2,102],[3,103])"
+    "{\"minx\": 1,\"miny\": 2,\"maxx\": 101,\"maxy\": 102}"
+    "[1,2,101,102]"
+    "[1,2,3,101,102,103]"
+    """
     name = 'Bounds'
 
     def convert(self, value, param, ctx):
@@ -31,6 +37,10 @@ def convert(self, value, param, ctx):
 
 
 class CRSParamType(click.ParamType):
+    """Click parameter type for the Coordinate Reference System of a project.
+
+    Accepts a string and returns an instance of the pyproj.CRS class.
+    """
     name = 'CRS'
 
     def convert(self, value, param, ctx) -> pyproj.CRS:
@@ -42,6 +52,10 @@ def convert(self, value, param, ctx) -> pyproj.CRS:
 
 
 class AttrParamType(click.ParamType):
+    """Click parameter for SilviMetric Attributes.
+
+    Returns list of PDAL dimensions that match the strings input.
+    """
     name = 'Attrs'
 
     # TODO add import similar to metrics
@@ -56,10 +70,15 @@ def convert(self, value, param, ctx) -> list[Attribute]:
         elif isinstance(value, str):
             return Attributes[value]
         else:
-            self.fail(f'{value!r} is of an invalid type, {e}', param, ctx)
+            self.fail(f'{value!r} is of an invalid type.', param, ctx)
 
 
 class MetricParamType(click.ParamType):
+    """Custom Click parameter type.
+
+    This param accepts names of metric groups or a path to a file containing
+    custom metrics.
+    """
     name = 'metrics'
 
     def convert(self, value, param, ctx) -> list[Metric]:
@@ -80,7 +99,7 @@ def convert(self, value, param, ctx) -> list[Metric]:
                     if not p.exists():
                         self.fail(
                             'Failed to find import file for metrics at'
-                            f' {str(p)}',
+                            f' {p}',
                             param,
                             ctx,
                         )
@@ -93,7 +112,7 @@ def convert(self, value, param, ctx) -> list[Metric]:
                     ms = user_metrics.metrics()
                 except Exception as e:
                     self.fail(
-                        f'Failed to import metrics from {str(p)} with error {e}',
+                        f'Failed to import metrics from {p} with error {e}',
                         param,
                         ctx,
                     )
@@ -126,7 +145,7 @@ def convert(self, value, param, ctx) -> list[Metric]:
                             metrics.add(m)
                         else:
                             metrics.udpate(list(m))
-                except Exception as e:
+                except Exception:
                     self.fail(
                         f'{val!r} is not available in Metrics', param, ctx
                     )

src/silvimetric/commands/extract.py

@@ -13,7 +13,8 @@
 
 np_to_gdal_types = {
     np.dtype(np.byte).str: gdal.GDT_Byte,
-    np.dtype(np.int8).str: gdal.GDT_Int16,
+    np.dtype(np.uint8).str: gdal.GDT_Byte,
+    np.dtype(np.int8).str: gdal.GDT_Int8,
     np.dtype(np.uint16).str: gdal.GDT_UInt16,
     np.dtype(np.int16).str: gdal.GDT_Int16,
     np.dtype(np.uint32).str: gdal.GDT_UInt32,
@@ -26,7 +27,12 @@
 
 
 def write_tif(
-    xsize: int, ysize: int, data: np.ndarray, name: str, config: ExtractConfig
+    xsize: int,
+    ysize: int,
+    data: np.ndarray,
+    nan_val: float | int,
+    name: str,
+    config: ExtractConfig,
 ) -> None:
     """
     Write out a raster with GDAL
@@ -55,11 +61,17 @@ def write_tif(
 
     driver = gdal.GetDriverByName('GTiff')
     gdal_type = np_to_gdal_types[np.dtype(data.dtype).str]
-    tif = driver.Create(str(path), int(xsize), int(ysize), 1, gdal_type)
+    tif = driver.Create(
+        str(path),
+        int(xsize),
+        int(ysize),
+        1,
+        gdal_type,
+    )
     tif.SetGeoTransform(transform)
     tif.SetProjection(srs.ExportToWkt())
     tif.GetRasterBand(1).WriteArray(data)
-    tif.GetRasterBand(1).SetNoDataValue(np.nan)
+    tif.GetRasterBand(1).SetNoDataValue(nan_val)
     tif.FlushCache()
     tif = None
 
@@ -84,8 +96,8 @@ def expl(x):
 
     attrs = [a.name for a in storage.config.attrs if a.name not in ['X', 'Y']]
 
-    # set index so we can apply to the whole dataset without needing to skip X and Y
-    # then reset in the index because that's what metric.do expects
+    # set index so we can apply to the whole dataset without needing to skip X
+    # and Y then reset in the index because that's what metric.do expects
     exploded = data_in.set_index(['X', 'Y']).apply(expl)[attrs].reset_index()
     metric_data = dask.persist(
         *[m.do(exploded) for m in storage.config.metrics]
@@ -99,9 +111,9 @@ def handle_overlaps(
     config: ExtractConfig, storage: Storage, indices: np.ndarray
 ) -> pd.DataFrame:
     """
-    Handle cells that have overlapping data. We have to re-run metrics over these
-    cells as there's no other accurate way to determined metric values. If there
-    are no overlaps, this will do nothing.
+    Handle cells that have overlapping data. We have to re-run metrics over
+    these cells as there's no other accurate way to determined metric values.
+    If there are no overlaps, this will do nothing.
 
     :param config: ExtractConfig.
     :param storage: Database storage object.
@@ -143,7 +155,7 @@ def handle_overlaps(
             else:
                 data = pd.concat([data, d])
 
-        # should find values that are not unique, meaning they have multiple entries
+        # find values that are not unique, means they have multiple entries
         data = data.set_index(['X', 'Y'])
         redo_indices = data.index[data.index.duplicated(keep='first')]
         if redo_indices.empty:
@@ -179,7 +191,12 @@ def extract(config: ExtractConfig) -> None:
     config.log.debug(f'Extracting metrics {[m for m in ma_list]}')
     root_bounds = storage.config.root
 
-    e = Extents(config.bounds, config.resolution, root=root_bounds)
+    e = Extents(
+        config.bounds,
+        config.resolution,
+        storage.config.alignment,
+        root=root_bounds,
+    )
     i = e.get_indices()
     xsize = e.x2
     ysize = e.y2
@@ -191,11 +208,17 @@ def extract(config: ExtractConfig) -> None:
     config.log.info(f'Writing rasters to {config.out_dir}')
     for ma in ma_list:
         # TODO should output in sections so we don't run into memory problems
-        m_data = np.full(
-            shape=(ysize, xsize), fill_value=np.nan, dtype=final[ma].dtype
-        )
+        dtype = final[ma].dtype
+        if dtype.kind in ['u', 'i']:
+            nan_val = np.iinfo(dtype).max
+        elif dtype.kind == 'f':
+            nan_val = np.nan
+        else:
+            raise ValueError('Invalid Raster data type {dtype}.')
+
+        m_data = np.full(shape=(ysize, xsize), fill_value=nan_val, dtype=dtype)
         a = final[['X', 'Y', ma]].to_numpy()
         for x, y, md in a[:]:
             m_data[int(y)][int(x)] = md
 
-        write_tif(xsize, ysize, m_data, ma, config)
+        write_tif(xsize, ysize, m_data, nan_val, ma, config)