Merge pull request #167 from IGNF/feat_color_offline

Feat color offline

Author: nlenglet-ign

.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:

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)

test/data/color/test_data_irc.tif

@@ -7,13 +7,19 @@
 import pytest
 
 from pdaltools import color
+from pdaltools.color import argument_parser
 
 cwd = os.getcwd()
 
 TEST_PATH = os.path.dirname(os.path.abspath(__file__))
 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"

test/data/color/test_data_rgb.tif

@@ -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)