From bb9e7b89c155d17f0423d118a6f081a493788cae Mon Sep 17 00:00:00 2001 From: Andrew Kirmse Date: Mon, 18 Mar 2024 19:51:15 -0700 Subject: [PATCH] Generated tiles always have elevation in meters. --- scripts/run_lidar_prominence.py | 42 +++++++++++++++++++++------------ 1 file changed, 27 insertions(+), 15 deletions(-) diff --git a/scripts/run_lidar_prominence.py b/scripts/run_lidar_prominence.py index 34541f9..6d73225 100644 --- a/scripts/run_lidar_prominence.py +++ b/scripts/run_lidar_prominence.py @@ -1,6 +1,6 @@ # Processes a set of Lidar tiles: # - Converts them to lat-long projection -# - Resamples to tiles that are of fixed size and resolution +# - Resamples to tiles that are of fixed size and resolution, with elevations in meters # # Writes a shapefiles "boundary.shp" to the tile dir, showing the nodata # boundary of the input. @@ -78,6 +78,18 @@ def write_multipolygon_to_file(multipolygon, filename): feature = None ds = None +def polygon_for_tile(x, y, size): + """Return a square with (x,y) as its lower left corner, size degrees in both directions""" + box = ogr.Geometry(ogr.wkbLinearRing) + box.AddPoint(x, y) + box.AddPoint(x, y + size) + box.AddPoint(x + size, y + size) + box.AddPoint(x + size, y) + box.AddPoint(x, y) + polybox = ogr.Geometry(ogr.wkbPolygon) + polybox.AddGeometry(box) + return polybox + def create_vrts(tile_dir, input_files): """ Creates virtual rasters (VRTs) for the input files on disk. @@ -163,13 +175,16 @@ def create_vrts(tile_dir, input_files): @handle_ctrl_c def process_tile(args): - (x, y, vrt_filename, output_filename) = args + """scale is multiplied by each input value""" + (x, y, vrt_filename, output_filename, scale) = args print(f"Processing {x:.2f}, {y:.2f}") gdal.UseExceptions() + translate_options = gdal.TranslateOptions( format = "EHdr", width = TILE_SIZE_SAMPLES, height = TILE_SIZE_SAMPLES, projWin = [x, y + TILE_SIZE_DEGREES, x + TILE_SIZE_DEGREES, y], + scaleParams = [[-30000, 30000, -30000*scale, 30000*scale]], callback=gdal.TermProgress_nocb) gdal.Translate(output_filename, vrt_filename, options = translate_options) @@ -191,7 +206,7 @@ def main(): parser.add_argument('--threads', default=1, type=int, help="Number of threads to use in computing prominence") parser.add_argument('--min_prominence', default=100, type=float, - help="Filter to this minimum prominence (units are --input_units") + help="Filter to this minimum prominence in meters") parser.add_argument('input_files', type=str, nargs='+', help='Input Lidar tiles, or GDAL VRT of tiles') @@ -233,17 +248,16 @@ def main(): x = xmin while x <= xmax - epsilon: # Skip this tile if it doesn't overlap any data in the source - box = ogr.Geometry(ogr.wkbLinearRing) - box.AddPoint(x, y) - box.AddPoint(x, y + TILE_SIZE_DEGREES) - box.AddPoint(x + TILE_SIZE_DEGREES, y + TILE_SIZE_DEGREES) - box.AddPoint(x + TILE_SIZE_DEGREES, y) - box.AddPoint(x, y) - polybox = ogr.Geometry(ogr.wkbPolygon) - polybox.AddGeometry(box) + polybox = polygon_for_tile(x, y, TILE_SIZE_DEGREES) if polybox.Intersects(bounds): output_filename = os.path.join(args.tile_dir, filename_for_coordinates(x, y)) - process_args.append((x, y, warped_vrt_filename, output_filename)) + # If input is in feet, will need to scale tile to meters + if args.input_units == "feet": + scale = 0.3048 + else: + scale = 1 + + process_args.append((x, y, warped_vrt_filename, output_filename, scale)) x += TILE_SIZE_DEGREES y += TILE_SIZE_DEGREES @@ -272,9 +286,7 @@ def main(): merge_inputs = os.path.join(args.output_dir, "*_*.dvt") merge_output = os.path.join(args.output_dir, "results") units_multiplier = 1 - if args.input_units == "feet" and args.output_units == "meters": - units_multiplier = 1.0 / 3.28084 - if args.input_units == "meters" and args.output_units == "feet": + if args.output_units == "feet": units_multiplier = 3.28084 merge_command = f"{merge_binary} --v=1 -f" + \ f" -t {args.threads} -m {args.min_prominence} -s {units_multiplier}" + \