Merge pull request #45 from IGNF/fix_color_reunion

Fix color reunion

CHANGELOG.md

@@ -1,3 +1,7 @@
+# 1.6.0
+- color: choose streams for RGB colorization, and IRC colorization (doc https://geoservices.ign.fr/services-web-experts-ortho)
+- color: detect white images.
+
 # 1.5.2
 - refactor tool to propagate header infos from one pipeline to another to use it by itself
 

pdaltools/_version.py

@@ -1,4 +1,4 @@
-__version__ = "1.5.2"
+__version__ = "1.6.0"
 
 
 if __name__ == "__main__":

pdaltools/color.py

@@ -3,8 +3,10 @@
 import time
 from math import ceil
 
+import numpy as np
 import pdal
 import requests
+from osgeo import gdal_array
 
 import pdaltools.las_info as las_info
 from pdaltools.unlock_file import copy_and_hack_decorator
@@ -58,7 +60,15 @@ def newfn(*args, **kwargs):
     return decorator
 
 
-def download_image_from_geoplateforme(proj, layer, minx, miny, maxx, maxy, pixel_per_meter, outfile, timeout):
+def is_image_white(filename: str):
+    raster_array = gdal_array.LoadFile(filename)
+    band_is_white = [np.all(band == 255) for band in raster_array]
+    return np.all(band_is_white)
+
+
+def download_image_from_geoplateforme(
+    proj, layer, minx, miny, maxx, maxy, pixel_per_meter, outfile, timeout, check_images
+):
     # Give single-point clouds a width/height of at least one pixel to have valid BBOX and SIZE
     if minx == maxx:
         maxx = minx + 1 / pixel_per_meter
@@ -88,6 +98,9 @@ def download_image_from_geoplateforme(proj, layer, minx, miny, maxx, maxy, pixel
     print(f"Ecriture du fichier: {outfile}")
     open(outfile, "wb").write(req.content)
 
+    if check_images and is_image_white(outfile):
+        raise ValueError(f"Downloaded image is white, with stream: {layer}")
+
 
 @copy_and_hack_decorator
 def color(
@@ -99,6 +112,9 @@ def color(
     color_rvb_enabled=True,
     color_ir_enabled=True,
     veget_index_file="",
+    check_images=False,
+    stream_RGB="ORTHOIMAGERY.ORTHOPHOTOS",
+    stream_IRC="ORTHOIMAGERY.ORTHOPHOTOS.IRC",
 ):
     metadata = las_info.las_info_metadata(input_file)
     minx, maxx, miny, maxy = las_info.get_bounds_from_header_info(metadata)
@@ -122,8 +138,9 @@ def color(
     if color_rvb_enabled:
         tmp_ortho = tempfile.NamedTemporaryFile()
         download_image_from_geoplateforme_retrying(
-            proj, "ORTHOIMAGERY.ORTHOPHOTOS", minx, miny, maxx, maxy, pixel_per_meter, tmp_ortho.name, timeout_second
+            proj, stream_RGB, minx, miny, maxx, maxy, pixel_per_meter, tmp_ortho.name, timeout_second, check_images
         )
+
         pipeline |= pdal.Filter.colorization(
             raster=tmp_ortho.name, dimensions="Red:1:256.0, Green:2:256.0, Blue:3:256.0"
         )
@@ -132,16 +149,9 @@ def color(
     if color_ir_enabled:
         tmp_ortho_irc = tempfile.NamedTemporaryFile()
         download_image_from_geoplateforme_retrying(
-            proj,
-            "ORTHOIMAGERY.ORTHOPHOTOS.IRC",
-            minx,
-            miny,
-            maxx,
-            maxy,
-            pixel_per_meter,
-            tmp_ortho_irc.name,
-            timeout_second,
+            proj, stream_IRC, minx, miny, maxx, maxy, pixel_per_meter, tmp_ortho_irc.name, timeout_second, check_images
         )
+
         pipeline |= pdal.Filter.colorization(raster=tmp_ortho_irc.name, dimensions="Infrared:1:256.0")
 
     pipeline |= pdal.Writer.las(
@@ -158,7 +168,7 @@ def color(
 
 
 def parse_args():
-    parser = argparse.ArgumentParser("Colorize tool")
+    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(
@@ -171,9 +181,39 @@ def parse_args():
     parser.add_argument(
         "--vegetation", type=str, default="", help="Vegetation file, value will be stored in Deviation field"
     )
+    parser.add_argument("--check-images", "-c", action="store_true", help="Check that downloaded image is not white")
+    parser.add_argument(
+        "--stream-RGB",
+        type=str,
+        default="ORTHOIMAGERY.ORTHOPHOTOS",
+        help="""WMS raster stream for RGB colorization:
+default stream (ORTHOIMAGERY.ORTHOPHOTOS) let the server choose the resolution
+for 20cm resolution rasters, use HR.ORTHOIMAGERY.ORTHOPHOTOS
+for 50 cm resolution rasters, use ORTHOIMAGERY.ORTHOPHOTOS.BDORTHO""",
+    )
+    parser.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""",
+    )
+
     return parser.parse_args()
 
 
 if __name__ == "__main__":
     args = parse_args()
-    color(args.input, args.output, args.proj, args.resolution, args.timeout, args.rvb, args.ir, args.vegetation)
+    color(
+        input_file=args.input,
+        output_file=args.output,
+        proj=args.proj,
+        pixel_per_meter=args.resolution,
+        timeout_second=args.timeout,
+        color_rvb_enabled=args.rvb,
+        color_ir_enabled=args.ir,
+        veget_index_file=args.vegetation,
+        check_images=args.check_images,
+        stream_RGB=args.stream_RGB,
+        stream_IRC=args.stream_IRC,
+    )

script/test/test_color_fail_white_image.sh

@@ -0,0 +1,19 @@
+mkdir tmp
+
+echo
+echo First test should faild with white image error.
+echo
+
+python -m pdaltools.color \
+-i ./test/data/sample_lareunion_epsg2975.laz \
+-o ./tmp/output.tif \
+--rvb -c
+
+echo 
+echo Second test should succeed, because we use RGB stream of 20 cm resolution.
+echo
+
+python -m pdaltools.color \
+-i ./test/data/sample_lareunion_epsg2975.laz \
+-o ./tmp/output.tif \
+--rvb -c --stream-RGB HR.ORTHOIMAGERY.ORTHOPHOTOS

script/test/test_download_reunion.sh

@@ -0,0 +1,19 @@
+mkdir tmp
+
+# test flux sur la reunion, flux 20cm:  HR.ORTHOIMAGERY.ORTHOPHOTOS => ok
+curl -o tmp/reunion1_20cm_ok.tif "https://data.geopf.fr/wms-r/wms?LAYERS=HR.ORTHOIMAGERY.ORTHOPHOTOS&EXCEPTIONS=text/xml&FORMAT=image/geotiff&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&STYLES=&CRS=EPSG:2975&BBOX=377000,7654000,377100,7654100&WIDTH=500&HEIGHT=500"
+
+# test flux sur la reunion, flux 50cm: ORTHOIMAGERY.ORTHOPHOTOS.BDORTHO => ok
+curl -o tmp/reunion1_50cm_ok.tif "https://data.geopf.fr/wms-r/wms?LAYERS=ORTHOIMAGERY.ORTHOPHOTOS.BDORTHO&EXCEPTIONS=text/xml&FORMAT=image/geotiff&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&STYLES=&CRS=EPSG:2975&BBOX=377000,7654000,377100,7654100&WIDTH=500&HEIGHT=500"
+
+# test flux sur la reunion, flux qui choisit: ORTHOIMAGERY.ORTHOPHOTOS => timeout
+curl -o tmp/reunion1_choix_timeout.txt "https://data.geopf.fr/wms-r/wms?LAYERS=ORTHOIMAGERY.ORTHOPHOTOS&EXCEPTIONS=text/xml&FORMAT=image/geotiff&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&STYLES=&CRS=EPSG:2975&BBOX=377000,7654000,377100,7654100&WIDTH=500&HEIGHT=500"
+
+# test flux sur la reunion, flux qui choisit: ORTHOIMAGERY.ORTHOPHOTOS => blanc
+curl -o tmp/reunion2_choix_blanc.tif "https://data.geopf.fr/wms-r/wms?LAYERS=ORTHOIMAGERY.ORTHOPHOTOS&EXCEPTIONS=text/xml&FORMAT=image/geotiff&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&STYLES=&CRS=EPSG:2975&BBOX=377000,7655950,377050,7655999.99&WIDTH=250&HEIGHT=250"
+
+# test flux sur la reunion, flux 20cm: HR.ORTHOIMAGERY.ORTHOPHOTOS => ok
+curl -o tmp/reunion2_20cm_ok.tif "https://data.geopf.fr/wms-r/wms?LAYERS=HR.ORTHOIMAGERY.ORTHOPHOTOS&EXCEPTIONS=text/xml&FORMAT=image/geotiff&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&STYLES=&CRS=EPSG:2975&BBOX=377000,7655950,377050,7655999.99&WIDTH=250&HEIGHT=250"
+
+# test flux sur la reunion, flux 50cm: ORTHOIMAGERY.ORTHOPHOTOS.BDORTHO => ok
+curl -o tmp/reunion2_50cm_ok.tif "https://data.geopf.fr/wms-r/wms?LAYERS=ORTHOIMAGERY.ORTHOPHOTOS.BDORTHO&EXCEPTIONS=text/xml&FORMAT=image/geotiff&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&STYLES=&CRS=EPSG:2975&BBOX=377000,7655950,377050,7655999.99&WIDTH=250&HEIGHT=250"

test/test_color.py

@@ -47,7 +47,7 @@ def test_epsg_fail():
 
 @pytest.mark.geopf
 def test_color_and_keeping_orthoimages():
-    tmp_ortho, tmp_ortho_irc = color.color(INPUT_PATH, OUTPUT_FILE, epsg)
+    tmp_ortho, tmp_ortho_irc = color.color(INPUT_PATH, OUTPUT_FILE, epsg, check_images=True)
     assert Path(tmp_ortho.name).exists()
     assert Path(tmp_ortho_irc.name).exists()
 
@@ -63,7 +63,7 @@ def test_color_narrow_cloud():
 @pytest.mark.geopf
 def test_download_image_ok():
     tif_output = os.path.join(TMPDIR, "download_image.tif")
-    color.download_image_from_geoplateforme(epsg, layer, minx, miny, maxx, maxy, pixel_per_meter, tif_output, 15)
+    color.download_image_from_geoplateforme(epsg, layer, minx, miny, maxx, maxy, pixel_per_meter, tif_output, 15, True)
 
 
 @pytest.mark.geopf
@@ -74,6 +74,27 @@ def test_color_epsg_2975_forced():
     color.color(input_path, output_path, 2975)
 
 
+def test_is_image_white_true():
+    input_path = os.path.join(TEST_PATH, "data/image/white.tif")
+    assert color.is_image_white(input_path), "This image should be detected as white"
+
+
+def test_is_image_white_false():
+    input_path = os.path.join(TEST_PATH, "data/image/colored.tif")
+    assert not color.is_image_white(input_path), "This image should NOT be detected as white"
+
+
+@pytest.mark.geopf
+def test_color_raise_for_white_image():
+    input_path = os.path.join(TEST_PATH, "data/sample_lareunion_epsg2975.laz")
+    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)
+
+    assert "Downloaded image is white" in str(excinfo.value)
+
+
 @pytest.mark.geopf
 def test_color_epsg_2975_detected():
     input_path = os.path.join(TEST_PATH, "data/sample_lareunion_epsg2975.laz")
@@ -86,22 +107,22 @@ def test_color_epsg_2975_detected():
 def test_download_image_raise1():
     retry_download = color.retry(2, 5)(color.download_image_from_geoplateforme)
     with pytest.raises(requests.exceptions.HTTPError):
-        retry_download(epsg, "MAUVAISE_COUCHE", minx, miny, maxx, maxy, pixel_per_meter, OUTPUT_FILE, 15)
+        retry_download(epsg, "MAUVAISE_COUCHE", minx, miny, maxx, maxy, pixel_per_meter, OUTPUT_FILE, 15, True)
 
 
 @pytest.mark.geopf
 def test_download_image_raise2():
     retry_download = color.retry(2, 5)(color.download_image_from_geoplateforme)
     with pytest.raises(requests.exceptions.HTTPError):
-        retry_download("9001", layer, minx, miny, maxx, maxy, pixel_per_meter, OUTPUT_FILE, 15)
+        retry_download("9001", layer, minx, miny, maxx, maxy, pixel_per_meter, OUTPUT_FILE, 15, True)
 
 
 def test_retry_on_server_error():
     with requests_mock.Mocker() as mock:
         mock.get(requests_mock.ANY, status_code=502, reason="Bad Gateway")
         with pytest.raises(requests.exceptions.HTTPError):
             retry_download = color.retry(2, 1, 2)(color.download_image_from_geoplateforme)
-            retry_download(epsg, layer, minx, miny, maxx, maxy, pixel_per_meter, OUTPUT_FILE, 15)
+            retry_download(epsg, layer, minx, miny, maxx, maxy, pixel_per_meter, OUTPUT_FILE, 15, True)
         history = mock.request_history
         assert len(history) == 3
 
@@ -111,7 +132,7 @@ def test_retry_on_connection_error():
         mock.get(requests_mock.ANY, exc=requests.exceptions.ConnectionError)
         with pytest.raises(requests.exceptions.ConnectionError):
             retry_download = color.retry(2, 1)(color.download_image_from_geoplateforme)
-            retry_download(epsg, layer, minx, miny, maxx, maxy, pixel_per_meter, OUTPUT_FILE, 15)
+            retry_download(epsg, layer, minx, miny, maxx, maxy, pixel_per_meter, OUTPUT_FILE, 15, True)
 
         history = mock.request_history
         assert len(history) == 3