diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 2fe0fe5..9aa3596 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -3,7 +3,7 @@ repos: rev: 23.12.0 hooks: - id: black - language_version: python3.11 + language_version: python3.12 - repo: https://github.com/pycqa/flake8 rev: 6.1.0 hooks: diff --git a/pdaltools/color.py b/pdaltools/color.py index 7f1750e..4a771f3 100644 --- a/pdaltools/color.py +++ b/pdaltools/color.py @@ -31,7 +31,7 @@ def match_min_max_with_pixel_size(min_d: float, max_d: float, pixel_per_meter: f return min_d, max_d -def color( +def color_from_stream( input_file: str, output_file: str, proj="", @@ -156,28 +156,58 @@ def color( return tmp_ortho, tmp_ortho_irc -def parse_args(): - parser = argparse.ArgumentParser("Colorize tool", formatter_class=argparse.RawTextHelpFormatter) - parser.add_argument("--input", "-i", type=str, required=True, help="Input file") - parser.add_argument("--output", "-o", type=str, default="", help="Output file") - parser.add_argument( - "--proj", "-p", type=str, default="", help="Projection, default will use projection from metadata input" +def color_from_files( + input_file: str, + output_file: str, + rgb_image: str, + irc_image: str, + color_rvb_enabled=True, + color_ir_enabled=True, + veget_index_file="", + vegetation_dim="Deviation", +): + pipeline = pdal.Reader.las(filename=input_file) + + writer_extra_dims = "all" + + if veget_index_file and veget_index_file != "": + print(f"Remplissage du champ {vegetation_dim} à partir du fichier {veget_index_file}") + pipeline |= pdal.Filter.colorization(raster=veget_index_file, dimensions=f"{vegetation_dim}:1:256.0") + writer_extra_dims = [f"{vegetation_dim}=ushort"] + + # Warning: the initial color is multiplied by 256 despite its initial 8-bits encoding + # which turns it to a 0 to 255*256 range. + # It is kept this way because of other dependencies that have been tuned to fit this range + if color_rvb_enabled: + pipeline |= pdal.Filter.colorization(raster=rgb_image, dimensions="Red:1:256.0, Green:2:256.0, Blue:3:256.0") + if color_ir_enabled: + pipeline |= pdal.Filter.colorization(raster=irc_image, dimensions="Infrared:1:256.0") + + pipeline |= pdal.Writer.las( + filename=output_file, extra_dims=writer_extra_dims, minor_version="4", dataformat_id="8", forward="all" ) - parser.add_argument("--resolution", "-r", type=float, default=5, help="Resolution, in pixel per meter") - parser.add_argument("--timeout", "-t", type=int, default=300, help="Timeout, in seconds") - parser.add_argument("--rvb", action="store_true", help="Colorize RVB") - parser.add_argument("--ir", action="store_true", help="Colorize IR") - parser.add_argument( - "--vegetation", - type=str, - default="", - help="Vegetation file (raster), value will be stored in 'vegetation_dim' field", + + print("Traitement du nuage de point") + pipeline.execute() + + +def argument_parser(): + parser = argparse.ArgumentParser("Colorize tool") + subparsers = parser.add_subparsers(required=True) + + # first command is 'from_stream' + from_stream = subparsers.add_parser("from_stream", help="Images are downloaded from streams") + from_stream.add_argument( + "--proj", "-p", type=str, default="", help="Projection, default will use projection from metadata input" ) - parser.add_argument( - "--vegetation_dim", type=str, default="Deviation", help="name of the extra_dim uses for the vegetation value" + from_stream.add_argument("--timeout", "-t", type=int, default=300, help="Timeout, in seconds") + from_stream.add_argument("--rvb", action="store_true", help="Colorize RVB") + from_stream.add_argument("--ir", action="store_true", help="Colorize IR") + from_stream.add_argument("--resolution", "-r", type=float, default=5, help="Resolution, in pixel per meter") + from_stream.add_argument( + "--check-images", "-c", action="store_true", help="Check that downloaded image is not white" ) - parser.add_argument("--check-images", "-c", action="store_true", help="Check that downloaded image is not white") - parser.add_argument( + from_stream.add_argument( "--stream-RGB", type=str, default="ORTHOIMAGERY.ORTHOPHOTOS", @@ -186,27 +216,49 @@ def parse_args(): for 20cm resolution rasters, use HR.ORTHOIMAGERY.ORTHOPHOTOS for 50 cm resolution rasters, use ORTHOIMAGERY.ORTHOPHOTOS.BDORTHO""", ) - parser.add_argument( + from_stream.add_argument( "--stream-IRC", type=str, default="ORTHOIMAGERY.ORTHOPHOTOS.IRC", help="""WMS raster stream for IRC colorization. Default to ORTHOIMAGERY.ORTHOPHOTOS.IRC Documentation about possible stream : https://geoservices.ign.fr/services-web-experts-ortho""", ) - parser.add_argument( + from_stream.add_argument( "--size-max-GPF", type=int, default=5000, help="Maximum edge size (in pixels) of downloaded images." " If input file needs more, several images are downloaded and merged.", ) + add_common_options(from_stream) + from_stream.set_defaults(func=from_stream_func) - return parser.parse_args() + # second command is 'from_files' + from_files = subparsers.add_parser("from_files", help="Images are in directories from RGB/IRC") + from_files.add_argument("--image_RGB", type=str, required=True, help="RGB image filepath") + from_files.add_argument("--image_IRC", type=str, required=True, help="IRC image filepath") + add_common_options(from_files) + from_files.set_defaults(func=from_files_func) + return parser -if __name__ == "__main__": - args = parse_args() - color( + +def add_common_options(parser): + parser.add_argument("--input", "-i", type=str, required=True, help="Input file") + parser.add_argument("--output", "-o", type=str, default="", help="Output file") + parser.add_argument( + "--vegetation", + type=str, + default="", + help="Vegetation file (raster), value will be stored in 'vegetation_dim' field", + ) + parser.add_argument( + "--vegetation_dim", type=str, default="Deviation", help="name of the extra_dim uses for the vegetation value" + ) + + +def from_stream_func(args): + color_from_stream( input_file=args.input, output_file=args.output, proj=args.proj, @@ -221,3 +273,33 @@ def parse_args(): stream_IRC=args.stream_IRC, size_max_gpf=args.size_max_GPF, ) + + +def from_files_func(args): + if args.image_RGB and args.image_RGB != "": + color_rvb_enabled = True + else: + color_rvb_enabled = False + if args.image_IRC and args.image_IRC != "": + color_irc_enabled = True + else: + color_irc_enabled = False + + if not color_rvb_enabled and not color_irc_enabled: + raise ValueError("At least one of --rvb or --ir must be provided") + + color_from_files( + input_file=args.input, + output_file=args.output, + rgb_image=args.image_RGB, + irc_image=args.image_IRC, + color_rvb_enabled=color_rvb_enabled, + color_ir_enabled=color_irc_enabled, + veget_index_file=args.vegetation, + vegetation_dim=args.vegetation_dim, + ) + + +if __name__ == "__main__": + args = argument_parser.parse_args() + args.func(args) diff --git a/test/data/color/test_data_irc.tif b/test/data/color/test_data_irc.tif new file mode 100644 index 0000000..20883e7 Binary files /dev/null and b/test/data/color/test_data_irc.tif differ diff --git a/test/data/color/test_data_rgb.tif b/test/data/color/test_data_rgb.tif new file mode 100644 index 0000000..2c43232 Binary files /dev/null and b/test/data/color/test_data_rgb.tif differ diff --git a/test/test_color.py b/test/test_color.py index 6c3b610..67d2161 100644 --- a/test/test_color.py +++ b/test/test_color.py @@ -7,6 +7,7 @@ import pytest from pdaltools import color +from pdaltools.color import argument_parser cwd = os.getcwd() @@ -14,6 +15,11 @@ TMPDIR = os.path.join(TEST_PATH, "tmp", "color") INPUT_PATH = os.path.join(TEST_PATH, "data/test_noepsg_043500_629205_IGN69.laz") +INPUT_PATH_TILE = os.path.join(TEST_PATH, "data/test_data_77055_627760_LA93_IGN69.laz") + +RGB_IMAGE = os.path.join(TEST_PATH, "data/color/test_data_rgb.tif") +IRC_IMAGE = os.path.join(TEST_PATH, "data/color/test_data_irc.tif") + OUTPUT_FILE = os.path.join(TMPDIR, "Semis_2021_0435_6292_LA93_IGN69.colorized.las") EPSG = "2154" @@ -33,12 +39,12 @@ def test_epsg_fail(): RuntimeError, match="EPSG could not be inferred from metadata: No 'srs' key in metadata.", ): - color.color(INPUT_PATH, OUTPUT_FILE, "", 0.1, 15) + color.color_from_stream(INPUT_PATH, OUTPUT_FILE, "", 0.1, 15) @pytest.mark.geopf def test_color_and_keeping_orthoimages(): - tmp_ortho, tmp_ortho_irc = color.color(INPUT_PATH, OUTPUT_FILE, EPSG, check_images=True) + tmp_ortho, tmp_ortho_irc = color.color_from_stream(INPUT_PATH, OUTPUT_FILE, EPSG, check_images=True) assert Path(tmp_ortho.name).exists() assert Path(tmp_ortho_irc.name).exists() @@ -63,7 +69,7 @@ def test_color_narrow_cloud(): input_path = os.path.join(TEST_PATH, "data/test_data_0436_6384_LA93_IGN69_single_point.laz") output_path = os.path.join(TMPDIR, "color_narrow_cloud_test_data_0436_6384_LA93_IGN69_single_point.colorized.laz") # Test that clouds that are smaller in width or height to 20cm are still colorized without an error. - color.color(input_path, output_path, EPSG) + color.color_from_stream(input_path, output_path, EPSG) with laspy.open(output_path, "r") as las: las_data = las.read() # Check all points are colored @@ -75,10 +81,9 @@ def test_color_narrow_cloud(): @pytest.mark.geopf def test_color_standard_cloud(): - input_path = os.path.join(TEST_PATH, "data/test_data_77055_627760_LA93_IGN69.laz") output_path = os.path.join(TMPDIR, "color_standard_cloud_test_data_77055_627760_LA93_IGN69.colorized.laz") # Test that clouds that are smaller in width or height to 20cm are still colorized without an error. - color.color(input_path, output_path, EPSG) + color.color_from_stream(INPUT_PATH_TILE, output_path, EPSG) with laspy.open(output_path, "r") as las: las_data = las.read() # Check all points are colored @@ -92,7 +97,7 @@ def test_color_epsg_2975_forced(): input_path = os.path.join(TEST_PATH, "data/sample_lareunion_epsg2975.laz") output_path = os.path.join(TMPDIR, "color_epsg_2975_forced_sample_lareunion_epsg2975.colorized.laz") - color.color(input_path, output_path, 2975) + color.color_from_stream(input_path, output_path, 2975) # the test is not working, the image is not detected as white @@ -104,7 +109,7 @@ def test_color_epsg_2975_forced(): # output_path = os.path.join(TMPDIR, "sample_lareunion_epsg2975.colorized.white.laz")# # with pytest.raises(ValueError) as excinfo: -# color.color(input_path, output_path, check_images=True) +# color.color_from_stream(input_path, output_path, check_images=True) # assert "Downloaded image is white" in str(excinfo.value) @@ -114,18 +119,17 @@ def test_color_epsg_2975_detected(): input_path = os.path.join(TEST_PATH, "data/sample_lareunion_epsg2975.laz") output_path = os.path.join(TMPDIR, "color_epsg_2975_detected_sample_lareunion_epsg2975.colorized.laz") # Test that clouds that are smaller in width or height to 20cm are still clorized without an error. - color.color(input_path, output_path) + color.color_from_stream(input_path, output_path) def test_color_vegetation_only(): - """Test the color() function with only vegetation""" - input_path = os.path.join(TEST_PATH, "data/test_data_77055_627760_LA93_IGN69.laz") + """Test the color_from_stream() function with only vegetation""" output_path = os.path.join(TMPDIR, "test_color_vegetation.colorized.las") vegetation_path = os.path.join(TEST_PATH, "data/mock_vegetation.tif") # Test with all parameters explicitly defined - color.color( - input_file=input_path, + color.color_from_stream( + input_file=INPUT_PATH_TILE, output_file=output_path, proj="2154", # EPSG:2154 (Lambert 93) color_rvb_enabled=False, # RGB enabled @@ -153,14 +157,13 @@ def test_color_vegetation_only(): @pytest.mark.geopf def test_color_with_all_parameters(): - """Test the color() function with all parameters specified""" - input_path = os.path.join(TEST_PATH, "data/test_data_77055_627760_LA93_IGN69.laz") + """Test the color_from_stream() function with all parameters specified""" output_path = os.path.join(TMPDIR, "test_color_all_params.colorized.las") vegetation_path = os.path.join(TEST_PATH, "data/mock_vegetation.tif") # Test with all parameters explicitly defined - tmp_ortho, tmp_ortho_irc = color.color( - input_file=input_path, + tmp_ortho, tmp_ortho_irc = color.color_from_stream( + input_file=INPUT_PATH_TILE, output_file=output_path, proj="2154", # EPSG:2154 (Lambert 93) pixel_per_meter=2.0, # custom resolution @@ -196,3 +199,56 @@ def test_color_with_all_parameters(): # Verify that the vegetation dimension is present assert "vegetation_dim" in las_data.point_format.dimension_names, "Vegetation dimension should be present" assert not np.all(las_data.vegetation_dim == 0), "Vegetation dimension should not be empty" + + +def test_color_from_files(): + output_path = os.path.join(TMPDIR, "color_standard_cloud_files_test_data_77055_627760_LA93_IGN69.colorized.laz") + + color.color_from_files(INPUT_PATH_TILE, output_path, RGB_IMAGE, IRC_IMAGE) + + assert os.path.exists(output_path) + + with laspy.open(output_path, "r") as las: + las_data = las.read() + + # Verify that all points have been colorized (no 0 values) + las_rgb_missing = (las_data.red == 0) & (las_data.green == 0) & (las_data.blue == 0) + assert not np.any(las_rgb_missing), f"No point should have missing RGB, found {np.count_nonzero(las_rgb_missing)}" + assert not np.any(las_data.nir == 0), "No point should have missing NIR" + + +@pytest.mark.geopf +def test_main_from_stream(): + output_file = os.path.join(TMPDIR, "main_from_stream", "output_main_from_stream.laz") + os.makedirs(os.path.dirname(output_file)) + cmd = f"from_stream -i {INPUT_PATH_TILE} -o {output_file} -p {EPSG} --rvb --ir".split() + args = argument_parser().parse_args(cmd) + args.func(args) + + assert os.path.exists(output_file) + + with laspy.open(output_file, "r") as las: + las_data = las.read() + + # Verify that all points have been colorized (no 0 values) + las_rgb_missing = (las_data.red == 0) & (las_data.green == 0) & (las_data.blue == 0) + assert not np.any(las_rgb_missing), f"No point should have missing RGB, found {np.count_nonzero(las_rgb_missing)}" + assert not np.any(las_data.nir == 0), "No point should have missing NIR" + + +def test_main_from_files(): + output_file = os.path.join(TMPDIR, "main_from_files", "output_main_from_files.laz") + os.makedirs(os.path.dirname(output_file)) + cmd = f"from_files -i {INPUT_PATH_TILE} -o {output_file} --image_RGB {RGB_IMAGE} --image_IRC {IRC_IMAGE}".split() + args = argument_parser().parse_args(cmd) + args.func(args) + + assert os.path.exists(output_file) + + with laspy.open(output_file, "r") as las: + las_data = las.read() + + # Verify that all points have been colorized (no 0 values) + las_rgb_missing = (las_data.red == 0) & (las_data.green == 0) & (las_data.blue == 0) + assert not np.any(las_rgb_missing), f"No point should have missing RGB, found {np.count_nonzero(las_rgb_missing)}" + assert not np.any(las_data.nir == 0), "No point should have missing NIR" diff --git a/test/test_unlock.py b/test/test_unlock.py index ecae35d..c41d4f8 100644 --- a/test/test_unlock.py +++ b/test/test_unlock.py @@ -4,7 +4,7 @@ import laspy import pytest -from pdaltools.color import color +from pdaltools.color import color_from_stream from pdaltools.las_info import las_info_metadata from pdaltools.unlock_file import copy_and_hack_decorator, unlock_file @@ -39,7 +39,7 @@ def test_copy_and_hack_decorator_color(): LAS_FILE = os.path.join(TMPDIR, "test_pdalfail_0643_6319_LA93_IGN69.las") # Color works only when an epsg is present in the header or as a parameter - color(LAZ_FILE, LAS_FILE, "2154", 1) + color_from_stream(LAZ_FILE, LAS_FILE, "2154", 1) las = laspy.read(LAS_FILE) print(las.header)