improve intensity normalization by stretching 3*stddev around the median

pyproject.toml

@@ -38,4 +38,5 @@ build-backend = "setuptools.build_meta"
 [project.scripts]
 trenchrun = "trenchrun.__main__:trenchrun"
 rls = "trenchrun.__main__:rls"
+trenchrun-builder = "trenchrun.__main__:pipline_builder"
 

src/trenchrun/__init__.py

@@ -1 +1 @@
-__version__='0.2.4'
+__version__='0.3.0'

src/trenchrun/__main__.py

@@ -6,4 +6,7 @@ def trenchrun():
     sys.exit(pl.trenchrun_main())
 
 def rls():
-    sys.exit(pl.rls_main())
+    sys.exit(pl.rls_main())
+
+def pipline_builder():
+    sys.exit(pl.pipeline_builder_main())

src/trenchrun/argparser.py

@@ -58,3 +58,24 @@ def get_rls_parser(args):
 
     args = parser.parse_args(args)
     return args
+
+def get_builder_parser(args):
+
+    import argparse
+
+    parser = argparse.ArgumentParser(description='Build a PDAL pipeline for a given EPT and GeoJSON geometry that is suitable for Trenchrun')
+    parser.add_argument('ept',
+                        help='EPT to use for fetching content', type=str)
+    parser.add_argument('geojson',
+                        help='GeoJSON geometry describing footprint', type=pathlib.Path)
+    parser.add_argument('--output',
+                        help='Output filename', type=pathlib.Path, default='output-pipeline.json')
+    parser.add_argument('--buffer',
+                        help='Buffer the geometry', type=float, default=float('nan') )
+    parser.add_argument('--debug',
+                        action='store_true',
+                        help='print debug messages to stderr')
+
+
+    args = parser.parse_args(args)
+    return args

src/trenchrun/blend.py

@@ -13,6 +13,7 @@
 
 from .data import Intensity, DSM, Daylight
 
+
 def readBand(filename,
              bandNumber,
              npDataType = np.float32,
@@ -23,7 +24,14 @@ def readBand(filename,
     projection = ds.GetProjection()
     array = band.ReadAsArray().astype(npDataType)
     if normalize:
-        array = (array - array.min()) / (array.max() - array.min())
+        norm = mpl.colors.Normalize()
+        med = np.median(array)
+        std = np.std(array)
+        factor = 3
+        min = med - factor*std
+        max = med + factor*std
+        norm = mpl.colors.Normalize(vmin=min, vmax=max, clip=True)
+        array = norm(array)
 
     # mask off any nodata values
     nodata = band.GetNoDataValue()
@@ -56,8 +64,7 @@ def write(self, args):
         dsm = readBand(str(self.dsm.path), 1)
 
         if intensity.shape != daylight.shape:
-            x, y = intensity.shape
-            daylight = daylight[0:x, 0:y]
+            raise Exception("Images are not the same size and shape!")
 
         logs.logger.info(f'Intensity shape {intensity.shape} ')
         logs.logger.info(f'Daylight shape {daylight.shape} ')
@@ -75,7 +82,10 @@ def write(self, args):
         cmap = mpl.cm.Greys_r
         daylight_RGB = cmap(daylight)
 
-        RGBA = intensity_RGBA * 0.5 + daylight_RGB * 0.5
+        try:
+            RGBA = intensity_RGBA * 0.5 + daylight_RGB * 0.5
+        except ValueError:
+            breakpoint()
 
         numBands = RGBA.shape[2]
 

src/trenchrun/builder.py

@@ -0,0 +1,51 @@
+
+from .logs import logger
+
+import pyproj
+from shapely.geometry import shape, Polygon, LineString
+from shapely.ops import transform
+import pdal
+
+import json
+import math
+
+def make_pipeline(args):
+    ept = pdal.Reader.ept(filename=args.ept, polygon = args.geometry.wkt, requests=16)
+    withheld = pdal.Filter.expression(expression="Withheld != 1")
+    no_noise = pdal.Filter.expression(expression = "!(Classification == 12 || Classification ==7)")
+    gdal = pdal.Writer.gdal(filename='dsm.tif', bounds = str(list(args.geometry.bounds)),
+                            data_type = "float32",
+                            dimension = "Z",
+                            output_type="idw",
+                            window_size=3,
+                            resolution = 0.5)
+    args.pipeline = ept | withheld | no_noise | gdal
+
+def reproject_geojson(args):
+    geojson_srs = pyproj.CRS("EPSG:4326")
+    ept_srs = pyproj.CRS("EPSG:3857")
+
+    geojson = args.geojson.open().read()
+    geojson = json.loads(geojson)
+    if 'features' in geojson:
+        geometry = geojson['features'][0]['geometry']
+    else:
+        geometry = geojson['geometry']
+    geojson = shape(geometry)
+
+    forward = pyproj.Transformer.from_crs(geojson_srs, ept_srs, always_xy=True).transform
+
+    geometry = transform(forward, geojson)
+    if (not math.isnan(args.buffer)):
+        geometry = geometry.buffer(args.buffer)
+    args.geometry = geometry
+
+def do(args):
+    logger.info(f"opening {args.ept} for processing")
+
+    reproject_geojson(args)
+    make_pipeline(args)
+    with args.output.open(mode='w') as f:
+        f.write(args.pipeline.pipeline)
+
+

src/trenchrun/data.py

@@ -7,6 +7,7 @@
 import pathlib
 
 import pdal
+from osgeo import gdal
 
 
 def run(command):
@@ -30,16 +31,20 @@ def __init__(self, reader):
         self.path = pathlib.Path(tempfile.NamedTemporaryFile(suffix='.tif', delete=False).name)
         self.validate()
 
+
     def __del__(self):
         self.path.unlink()
 
+    def getMetadata(self):
+        pass
+
     def getStage(self):
         pass
 
     def validate(self):
         pass
 
-    def process(self):
+    def process(self, *args, **kwargs):
         stage = self.getStage()
 
         pipeline = self.reader.get() | stage
@@ -49,7 +54,11 @@ def process(self):
             count = pipeline.execute_streaming(chunk_size=self.reader.args.chunk_size)
         else:
             count = pipeline.execute()
+
         logs.logger.info(f'Wrote {self.name} product for {count} points')
+        self.metadata = pipeline.metadata
+        self.getMetadata()
+
 
 class Intensity(Product):
     def __init__(self, reader):
@@ -63,6 +72,10 @@ def getStage(self):
             data_type='uint16_t',
             dimension='Intensity',
             output_type = 'idw',
+            origin_x = self.reader.args.origin_x,
+            origin_y = self.reader.args.origin_y,
+            width = self.reader.args.width,
+            height = self.reader.args.height,
             resolution=self.reader.args.resolution,
         )
 
@@ -90,6 +103,20 @@ def getStage(self):
                                  resolution=self.reader.args.resolution)
         return stage
 
+    def getMetadata(self):
+        ds = gdal.Open(self.path)
+        xmin, xpixel, _, ymax, _, ypixel = ds.GetGeoTransform()
+        width, height = ds.RasterXSize, ds.RasterYSize
+        # xmax = xmin + width * xpixel
+        ymin = ymax + height * ypixel
+
+        self.reader.args.origin_x = xmin
+        self.reader.args.origin_y = ymin
+        self.reader.args.width = width
+        self.reader.args.height = height
+
+
+
 class Daylight(object):
     def __init__(self, dsm: DSM):
         self.name = 'Daylight'

src/trenchrun/do.py

@@ -1,6 +1,6 @@
 from . import data
 from . import blend
-from .logs import logger 
+from .logs import logger
 
 
 def doIt(args):

src/trenchrun/main.py

@@ -1,8 +1,9 @@
 
-from . import logs 
+from . import logs
 from . import do
 from . import argparser
 from . import rls
+from . import builder
 import sys
 
 def rls_main():
@@ -21,6 +22,14 @@ def trenchrun_main():
 
     do.doIt(args)
 
+def pipeline_builder_main():
+    args = argparser.get_builder_parser(sys.argv[1:])
+    if args.debug:
+        logs.logger.setLevel(logs.logging.DEBUG)
+        logs.handler.setLevel(logs.logging.DEBUG)
+
+    builder.do(args)
+
 
 if __name__ == "__main__":
     sys.exit(trenchrun_main())

src/trenchrun/rls.py

@@ -15,8 +15,6 @@
 from shapely.ops import transform
 from shapely import wkt
 
-
-
 def get_polyline(args):
 
     line_dd = wkt.loads(args.line)