diff --git a/pyproject.toml b/pyproject.toml index 8941e34..c8de4d7 100644 --- a/pyproject.toml +++ b/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" diff --git a/src/trenchrun/__init__.py b/src/trenchrun/__init__.py index 2de26d7..1ca87cf 100644 --- a/src/trenchrun/__init__.py +++ b/src/trenchrun/__init__.py @@ -1 +1 @@ -__version__='0.2.4' +__version__='0.3.0' diff --git a/src/trenchrun/__main__.py b/src/trenchrun/__main__.py index 340d168..b210546 100644 --- a/src/trenchrun/__main__.py +++ b/src/trenchrun/__main__.py @@ -6,4 +6,7 @@ def trenchrun(): sys.exit(pl.trenchrun_main()) def rls(): - sys.exit(pl.rls_main()) \ No newline at end of file + sys.exit(pl.rls_main()) + +def pipline_builder(): + sys.exit(pl.pipeline_builder_main()) \ No newline at end of file diff --git a/src/trenchrun/argparser.py b/src/trenchrun/argparser.py index 4ba5bfd..04e1a14 100644 --- a/src/trenchrun/argparser.py +++ b/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 \ No newline at end of file diff --git a/src/trenchrun/blend.py b/src/trenchrun/blend.py index a9ec96d..c6a5c8c 100644 --- a/src/trenchrun/blend.py +++ b/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] diff --git a/src/trenchrun/builder.py b/src/trenchrun/builder.py new file mode 100644 index 0000000..439c8e7 --- /dev/null +++ b/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) + + diff --git a/src/trenchrun/data.py b/src/trenchrun/data.py index 68184d6..298057d 100644 --- a/src/trenchrun/data.py +++ b/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' diff --git a/src/trenchrun/do.py b/src/trenchrun/do.py index 300610c..878b77e 100644 --- a/src/trenchrun/do.py +++ b/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): diff --git a/src/trenchrun/main.py b/src/trenchrun/main.py index a4fea7f..602d297 100644 --- a/src/trenchrun/main.py +++ b/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()) diff --git a/src/trenchrun/rls.py b/src/trenchrun/rls.py index 566142b..4ef939a 100644 --- a/src/trenchrun/rls.py +++ b/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)