Support for orthophotos with more than 3 bands, uint16 cogeo workarounds

app/api/formulas.py

@@ -94,6 +94,13 @@ camera_filters = [
     'NRG',
     'NRB',
 
+    'RGBN',
+
+    'BGRNRe',
+    'BGRReN',
+    'RGBNRe',
+    'RGBReN',
+
     # more?
     # TODO: certain cameras have only two bands? eg. MAPIR NDVI BLUE+NIR
 ]
@@ -123,22 +130,25 @@ def lookup_formula(algo, band_order = 'RGB'):
     return expr, hrange
 
 @lru_cache(maxsize=2)
-def get_algorithm_list():
-    return [{'id': k, 'filters': get_camera_filters_for(algos[k]), **algos[k]} for k in algos if not k.startswith("_")]
+def get_algorithm_list(max_bands=3):
+    return [{'id': k, 'filters': get_camera_filters_for(algos[k], max_bands), **algos[k]} for k in algos if not k.startswith("_")]
 
-def get_camera_filters_for(algo):
+def get_camera_filters_for(algo, max_bands=3):
     result = []
     expr = algo['expr']
-    bands = list(set(re.findall("([A-Z]+?[a-z]*)", expr)))
+    pattern = re.compile("([A-Z]+?[a-z]*)")
+    bands = list(set(re.findall(pattern, expr)))
     for f in camera_filters:
         # Count bands that show up in the filter
         count = 0
-        for b in f:
+        fbands = list(set(re.findall(pattern, f)))
+
+        for b in fbands:
             if b in bands:
                 count += 1
 
         # If all bands are accounted for, this is a valid filter for this algo
-        if count >= len(bands):
+        if count >= len(bands) and len(fbands) <= max_bands:
             result.append(f)
 
     return result

app/api/tiler.py

@@ -1,11 +1,15 @@
+import sys
+sys.path.insert(0, "/webodm/tmp/rio-tiler")
+
 import rasterio
+from rasterio.enums import ColorInterp
 import urllib
 import os
 from django.http import HttpResponse
 from rio_tiler.errors import TileOutsideBounds
 from rio_tiler.mercator import get_zooms
 from rio_tiler import main
-from rio_tiler.utils import array_to_image, get_colormap, expression, linear_rescale, _chunks, _apply_discrete_colormap
+from rio_tiler.utils import array_to_image, get_colormap, expression, linear_rescale, _chunks, _apply_discrete_colormap, has_alpha_band
 from rio_tiler.profiles import img_profiles
 
 import numpy as np
@@ -158,7 +162,13 @@ class Metadata(TaskNestedView):
             raise exceptions.NotFound()
 
         try:
-            info = main.metadata(raster_path, pmin=pmin, pmax=pmax, histogram_bins=255, histogram_range=hrange, expr=expr)
+            with rasterio.open(raster_path, "r") as src:
+                band_count = src.meta['count']
+                if has_alpha_band(src):
+                    band_count -= 1
+                print(band_count)
+
+                info = main.metadata(src, pmin=pmin, pmax=pmax, histogram_bins=255, histogram_range=hrange, expr=expr)
         except IndexError as e:
             # Caught when trying to get an invalid raster metadata
             raise exceptions.ValidationError("Cannot retrieve raster metadata: %s" % str(e))
@@ -186,7 +196,7 @@ class Metadata(TaskNestedView):
             colormaps = ['jet', 'terrain', 'gist_earth', 'pastel1']
         elif formula and bands:
             colormaps = ['rdylgn', 'spectral', 'rdylgn_r', 'spectral_r']
-            algorithms = *get_algorithm_list(),
+            algorithms = *get_algorithm_list(band_count),
 
         info['color_maps'] = []
         info['algorithms'] = algorithms
@@ -308,9 +318,29 @@ class Tiles(TaskNestedView):
 
         with rasterio.open(url) as src:
             minzoom, maxzoom = get_zoom_safe(src)
+            has_alpha = has_alpha_band(src)
             if z < minzoom - ZOOM_EXTRA_LEVELS or z > maxzoom + ZOOM_EXTRA_LEVELS:
                 raise exceptions.NotFound()
 
+            # Handle N-bands datasets
+            if tile_type == 'orthophoto':
+                ci = src.colorinterp
+
+                # More than 4 bands?
+                if len(ci) > 4:
+                    # Try to find RGBA band order
+                    if ColorInterp.red in ci and \
+                        ColorInterp.green in ci and \
+                        ColorInterp.blue in ci: # and ColorInterp.alpha in ci:
+                        indexes = (ci.index(ColorInterp.red) + 1,
+                                   ci.index(ColorInterp.green) + 1,
+                                   ci.index(ColorInterp.blue) + 1,)
+                        # TODO: adding alpha band should fix black backgrounds
+                        # but the tiles disappear. Probable bug in rasterio/GDAL
+                    else:
+                        # Fallback to first four
+                        indexes = (1, 2, 3, ) # , 4, )
+
         resampling="nearest"
         padding=0
         if tile_type in ["dsm", "dtm"]:
@@ -330,7 +360,7 @@ class Tiles(TaskNestedView):
             raise exceptions.NotFound("Outside of bounds")
 
         # Use alpha channel for transparency, don't use the mask if one is provided (redundant)
-        if tile.shape[0] == 4:
+        if has_alpha and expr is None:
             mask = None
 
         if color_map:

app/cogeo.py

@@ -4,6 +4,7 @@ import tempfile
 import shutil
 import rasterio
 from rio_cogeo.cogeo import cog_validate, cog_translate
+from rio_tiler.utils import has_alpha_band
 from webodm import settings
 
 logger = logging.getLogger('app.logger')
@@ -55,7 +56,11 @@ def assure_cogeo(src_path):
             GDAL_TIFF_OVR_BLOCKSIZE="128",
         )
 
-        cog_translate(dst, tmpfile, output_profile,
+        nodata = None
+        if has_alpha_band(dst) and dst.meta['dtype'] == 'uint16':
+            nodata = 0.0 # Hack to workaround https://github.com/cogeotiff/rio-cogeo/issues/112
+
+        cog_translate(dst, tmpfile, output_profile, nodata=nodata,
                       config=config, in_memory=False,
                       quiet=True, web_optimized=True)
 

app/tests/test_formulas.py

@@ -1,3 +1,4 @@
+import re
 from django.test import TestCase
 from app.api.formulas import lookup_formula, get_algorithm_list, get_camera_filters_for, algos
 
@@ -25,10 +26,26 @@ class TestFormulas(TestCase):
         self.assertTrue(lookup_formula("_TESTFUNC", "RGB")[1] == None)
 
     def test_algo_list(self):
-        list = get_algorithm_list()
+        al = get_algorithm_list()
 
-        # Do not show test algos
-        for i in list:
+        pattern = re.compile("([A-Z]+?[a-z]*)")
+        for i in al:
+            # Do not show test algos
             self.assertFalse(i['id'].startswith("_"))
 
-        self.assertTrue(get_camera_filters_for(algos['VARI']) == ['RGB'])
+            # Filters are less than 3 bands
+            for f in i['filters']:
+                bands = list(set(re.findall(pattern, f)))
+                self.assertTrue(len(bands) <= 3)
+
+        self.assertTrue(get_camera_filters_for(algos['VARI']) == ['RGB'])
+
+        # Request algorithms with more band filters
+        al = get_algorithm_list(max_bands=5)
+
+        pattern = re.compile("([A-Z]+?[a-z]*)")
+        for i in al:
+            # Filters are less than 5 bands
+            for f in i['filters']:
+                bands = list(set(re.findall(pattern, f)))
+                self.assertTrue(len(bands) <= 5)