working on #364

cate/util/im/cmap_lc.py

@@ -0,0 +1,88 @@
+# The MIT License (MIT)
+# Copyright (c) 2016, 2017 by the ESA CCI Toolbox development team and contributors
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy of
+# this software and associated documentation files (the "Software"), to deal in
+# the Software without restriction, including without limitation the rights to
+# use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
+# of the Software, and to permit persons to whom the Software is furnished to do
+# so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in all
+# copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+
+import matplotlib
+import matplotlib.cm as cm
+
+__author__ = "Norman Fomferra (Brockmann Consult GmbH)"
+
+LAND_COVER_CCI_CMAP = 'land_cover_cci'
+
+
+def register_lc_color_map():
+    lc_color_mappings = [
+        (0, dict(r=0, g=0, b=0)),
+        (10, dict(r=255, g=255, b=100)),
+        (11, dict(r=255, g=255, b=100)),
+        (12, dict(r=255, g=255, b=0)),
+        (20, dict(r=170, g=240, b=240)),
+        (30, dict(r=220, g=240, b=100)),
+        (40, dict(r=200, g=200, b=100)),
+        (50, dict(r=0, g=100, b=0)),
+        (60, dict(r=0, g=160, b=0)),
+        (61, dict(r=0, g=160, b=0)),
+        (62, dict(r=170, g=200, b=0)),
+        (70, dict(r=0, g=60, b=0)),
+        (71, dict(r=0, g=60, b=0)),
+        (72, dict(r=0, g=80, b=0)),
+        (80, dict(r=40, g=80, b=0)),
+        (81, dict(r=40, g=80, b=0)),
+        (82, dict(r=40, g=100, b=0)),
+        (90, dict(r=120, g=130, b=0)),
+        (100, dict(r=140, g=160, b=0)),
+        (110, dict(r=190, g=150, b=0)),
+        (120, dict(r=150, g=100, b=0)),
+        (121, dict(r=120, g=75, b=0)),
+        (122, dict(r=150, g=100, b=0)),
+        (130, dict(r=255, g=180, b=50)),
+        (140, dict(r=255, g=220, b=210)),
+        (150, dict(r=255, g=235, b=175)),
+        (151, dict(r=255, g=205, b=120)),
+        (152, dict(r=255, g=210, b=120)),
+        (153, dict(r=255, g=235, b=175)),
+        (160, dict(r=0, g=120, b=190)),
+        (170, dict(r=0, g=150, b=120)),
+        (180, dict(r=0, g=220, b=130)),
+        (190, dict(r=195, g=20, b=0)),
+        (200, dict(r=255, g=245, b=215)),
+        (201, dict(r=220, g=220, b=220)),
+        (202, dict(r=255, g=245, b=215)),
+        (210, dict(r=0, g=70, b=200)),
+        (220, dict(r=255, g=255, b=255)),
+    ]
+
+    classes = {lc: color for lc, color in lc_color_mappings}
+
+    invalid_rgba = (0, 0, 0, 0.5)
+    class_0_rgba = (0, 0, 0, 0)
+
+    rgba_list = []
+    num_entries = 256
+    last_rgba = invalid_rgba
+    for i in range(num_entries):
+        color = classes.get(i)
+        if color:
+            last_rgba = (color['r'] / 255, color['g'] / 255, color['b'] / 255, 1.0)
+        rgba_list.append(last_rgba)
+    rgba_list[0] = class_0_rgba
+
+    cmap = matplotlib.colors.ListedColormap(rgba_list, name=LAND_COVER_CCI_CMAP, N=num_entries)
+    cm.register_cmap(cmap=cmap)

cate/util/im/cmaps.py

@@ -28,6 +28,8 @@
 import numpy as np
 from PIL import Image
 
+from cate.util.im.cmap_lc import register_lc_color_map, LAND_COVER_CCI_CMAP
+
 __author__ = "Norman Fomferra (Brockmann Consult GmbH)"
 
 # Have colormaps separated into categories:
@@ -68,7 +70,7 @@
             'Pastel2', 'Set1', 'Set2', 'Set3')),
           ('Miscellaneous',
            'Colormaps that don\'t fit into the categories above.',
-           ('gist_earth', 'terrain', 'ocean', 'gist_stern',
+           (LAND_COVER_CCI_CMAP, 'gist_earth', 'terrain', 'ocean', 'gist_stern',
             'brg', 'CMRmap', 'cubehelix',
             'gnuplot', 'gnuplot2', 'gist_ncar',
             'nipy_spectral', 'jet', 'rainbow',
@@ -85,10 +87,14 @@ def get_cmaps():
     <cbar-png-bytes> are encoded PNG images of size 256 x 2 pixels,
     :return: all known matplotlib color maps
     """
+
     global _CBARS_LOADED, _CMAPS
     if not _CBARS_LOADED:
         _LOCK.acquire()
         _CBARS_LOADED = True
+
+        register_lc_color_map()
+
         new_cmaps = []
         for cmap_category, cmap_description, cmap_names in _CMAPS:
             cbar_list = []
@@ -129,8 +135,6 @@ def get_cmaps():
                         if a > 1.0:
                             a = 1.0
                     new_cmap = matplotlib.colors.ListedColormap(new_colors, name=new_name)
-                    print("WARNING: could not create colormap '{}' because '{}' is has type ListedColormap"
-                          .format(new_name, cmap.name))
                     cm.register_cmap(cmap=new_cmap)
                 else:
                     new_name = cmap.name + '_alpha' if hasattr(cmap, 'name') else 'unknown'
@@ -156,6 +160,7 @@ def get_cmaps():
 
                 cbar_list.append((cmap_name, cbar_png_bytes))
             new_cmaps.append((cmap_category, cmap_description, tuple(cbar_list)))
+
         _CMAPS = tuple(new_cmaps)
         _LOCK.release()
         # import pprint

test/util/im/test_utils.py

@@ -83,6 +83,94 @@ def test_num_levels_min(self):
             utils.cardinal_log2(n * ts)
             # print(s, ts, n, n * ts - s, l2)
 
+    def test_f_and_g(self):
+        import numba as nb
+        from numpy.testing import assert_almost_equal
+
+        @nb.jit(nopython=True)
+        def pyramid_subdivision_count(s_max: int, ts: int, ntl0_max: int = 1):
+            """
+            Compute number of times *w* can be divided by 2 without remainder and while the result is still
+            integer-dividable by *ts*.
+            """
+            count = 0
+            s = s_max
+            while s % 2 == 0 and s % ts == 0 and (s // ts) % 2 == 0 and (s // ts) > ntl0_max:
+                s //= 2
+                count += 1
+            return count
+
+        @nb.jit(nopython=True)
+        def pyramid_subdivision(w_max: int, h_max: int,
+                                ts_min: int, ts_max: int,
+                                tw_out, th_out,
+                                ntl0x_max: int = 1,
+                                ntl0y_max: int = 1):
+            size = ts_max - ts_min + 1
+
+            cx = np.empty(ts_max - ts_min + 1, dtype=np.int32)
+            cy = np.empty(ts_max - ts_min + 1, dtype=np.int32)
+            for i in range(size):
+                ts = ts_min + i
+                cx[i] = pyramid_subdivision_count(w_max, ts, ntl0_max=ntl0x_max)
+                cy[i] = pyramid_subdivision_count(h_max, ts, ntl0_max=ntl0y_max)
+
+            cx_max = -1
+            cy_max = -1
+            for i in range(size):
+                cx_max = max(cx[i], cx_max)
+                cy_max = max(cy[i], cy_max)
+
+            c = min(cx_max, cy_max)
+
+            for ix in range(tw_out.size):
+                tw_out[ix] = 0
+            for iy in range(th_out.size):
+                th_out[iy] = 0
+
+            if c <= 0:
+                return 0
+
+            ix = 0
+            iy = 0
+            for i in range(size):
+                if cx[i] >= c and ix < tw_out.size:
+                    tw_out[ix] = ts_min + i
+                    ix += 1
+                if cy[i] >= c and iy < th_out.size:
+                    th_out[iy] = ts_min + i
+                    iy += 1
+
+            return c
+
+        w = 129600
+        h = w // 2
+
+        self.assertEqual(pyramid_subdivision_count(w, 45), 6)
+        self.assertEqual(pyramid_subdivision_count(w, 90), 5)
+        self.assertEqual(pyramid_subdivision_count(w, 100), 4)
+        self.assertEqual(pyramid_subdivision_count(w, 180), 4)
+        self.assertEqual(pyramid_subdivision_count(w, 225), 6)
+        self.assertEqual(pyramid_subdivision_count(w, 256), 0)
+        self.assertEqual(pyramid_subdivision_count(w, 405), 6)
+        self.assertEqual(pyramid_subdivision_count(w, 675), 6)
+        self.assertEqual(pyramid_subdivision_count(w, 1024), 0)
+        self.assertEqual(pyramid_subdivision_count(w, 1800), 3)
+        self.assertEqual(pyramid_subdivision_count(w, 2048), 0)
+
+        tw_out = np.zeros(10, dtype=np.int32)
+        th_out = np.zeros(10, dtype=np.int32)
+
+        c = pyramid_subdivision(w, h, 180, 2048, tw_out, th_out, ntl0x_max=1, ntl0y_max=1)
+        self.assertEqual(c, 5)
+        assert_almost_equal(tw_out, np.array([225, 270, 405, 450, 675, 810, 1350, 2025, 0, 0]))
+        assert_almost_equal(th_out, np.array([225, 405, 675, 2025, 0, 0, 0, 0, 0, 0]))
+
+        c = pyramid_subdivision(w, h, 180, 2048, tw_out, th_out, ntl0x_max=2, ntl0y_max=1)
+        self.assertEqual(c, 5)
+        assert_almost_equal(tw_out, np.array([225, 270, 405, 450, 675, 810, 1350, 2025, 0, 0]))
+        assert_almost_equal(th_out, np.array([225, 405, 675, 2025, 0, 0, 0, 0, 0, 0]))
+
     def test_chunk_size(self):
         # print('----------- test_chunk_size:')
         for s in range(270, 64 * 270, 270):