🐛 fix grid creation bit precision (close #210); 🐛 fix rio cogeo lib (c…

…lose #209) and migrate to https bdc-catalog library (close #211)

cube_builder/controller.py

@@ -6,27 +6,25 @@
 # under the terms of the MIT License; see LICENSE file for more details.
 #
 
+
 """Define Cube Builder business interface."""
 from copy import deepcopy
 from datetime import datetime
 from typing import Tuple, Union
 
 # 3rdparty
-import rasterio
 import sqlalchemy
 from bdc_catalog.models import (Band, BandSRC, Collection, CompositeFunction, GridRefSys, Item, MimeType, Quicklook,
                                 ResolutionUnit, SpatialRefSys, Tile, db)
 from geoalchemy2 import func
-from geoalchemy2.shape import from_shape
 from rasterio.crs import CRS
-from rasterio.warp import transform
-from shapely.geometry import Polygon
 from werkzeug.exceptions import NotFound, abort
 
 from .constants import (CLEAR_OBSERVATION_ATTRIBUTES, CLEAR_OBSERVATION_NAME, COG_MIME_TYPE, DATASOURCE_ATTRIBUTES,
                         PROVENANCE_ATTRIBUTES, PROVENANCE_NAME, SRID_ALBERS_EQUAL_AREA, TOTAL_OBSERVATION_ATTRIBUTES,
                         TOTAL_OBSERVATION_NAME)
 from .forms import CollectionForm
+from .grids import create_grids
 from .models import Activity, CubeParameters
 from .utils.image import validate_merges
 from .utils.processing import get_cube_parts, get_or_create_model
@@ -513,114 +511,30 @@ def get_grs_schema(cls, grs_id):
     @classmethod
     def create_grs_schema(cls, name, description, projection, meridian, degreesx, degreesy, bbox, srid=100001):
         """Create a Brazil Data Cube Grid Schema."""
-        bbox = bbox.split(',')
-        bbox_obj = {
-            "w": float(bbox[0]),
-            "n": float(bbox[1]),
-            "e": float(bbox[2]),
-            "s": float(bbox[3])
-        }
-        tile_srs_p4 = "+proj=longlat +ellps=GRS80 +no_defs"
-        if projection == 'aea':
-            tile_srs_p4 = "+proj=aea +lat_0=-12 +lon_0={} +lat_1=-2 +lat_2=-22 +x_0=5000000 +y_0=10000000 +ellps=GRS80 +units=m +no_defs".format(meridian)
-        elif projection == 'sinu':
-            tile_srs_p4 = "+proj=sinu +lon_0={} +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs".format(meridian)
-
-        # Number of tiles and base tile
-        num_tiles_x = int(360. / degreesx)
-        num_tiles_y = int(180. / degreesy)
-        h_base = num_tiles_x / 2
-        v_base = num_tiles_y / 2
-
-        # Tile size in meters (dx,dy) at center of system (argsmeridian,0.)
-        src_crs = '+proj=longlat +ellps=GRS80 +no_defs'
-        dst_crs = tile_srs_p4
-        xs = [(meridian - degreesx / 2), (meridian + degreesx / 2), meridian, meridian, 0.]
-        ys = [0., 0., -degreesy / 2, degreesy / 2, 0.]
-        out = transform(CRS.from_proj4(src_crs), CRS.from_proj4(dst_crs), xs, ys, zs=None)
-        x1 = out[0][0]
-        x2 = out[0][1]
-        y1 = out[1][2]
-        y2 = out[1][3]
-        dx = x2 - x1
-        dy = y2 - y1
-
-        # Coordinates of WRS center (0.,0.) - top left coordinate of (h_base,v_base)
-        x_center = out[0][4]
-        y_center = out[1][4]
-        # Border coordinates of WRS grid
-        x_min = x_center - dx * h_base
-        y_max = y_center + dy * v_base
-
-        # Upper Left is (xl,yu) Bottom Right is (xr,yb)
-        xs = [bbox_obj['w'], bbox_obj['e'], meridian, meridian]
-        ys = [0., 0., bbox_obj['n'], bbox_obj['s']]
-        out = transform(src_crs, dst_crs, xs, ys, zs=None)
-        xl = out[0][0]
-        xr = out[0][1]
-        yu = out[1][2]
-        yb = out[1][3]
-        h_min = int((xl - x_min) / dx)
-        h_max = int((xr - x_min) / dx)
-        v_min = int((y_max - yu) / dy)
-        v_max = int((y_max - yb) / dy)
-
-        tiles = []
-        features = []
-        dst_crs = '+proj=longlat +ellps=GRS80 +no_defs'
-        src_crs = tile_srs_p4
-
-        for ix in range(h_min, h_max+1):
-            x1 = x_min + ix*dx
-            x2 = x1 + dx
-            for iy in range(v_min, v_max+1):
-                y1 = y_max - iy*dy
-                y2 = y1 - dy
-                # Evaluate the bounding box of tile in longlat
-                xs = [x1, x2, x2, x1]
-                ys = [y1, y1, y2, y2]
-                out = rasterio.warp.transform(src_crs, dst_crs, xs, ys, zs=None)
-
-                polygon = from_shape(
-                    Polygon(
-                        [
-                            (x1, y2),
-                            (x2, y2),
-                            (x2, y1),
-                            (x1, y1),
-                            (x1, y2)
-                        ]
-                    ),
-                    srid=SRID_ALBERS_EQUAL_AREA
-                )
-
-                # Insert tile
-                tile_name = '{0:03d}{1:03d}'.format(ix, iy)
-                tiles.append(dict(
-                    name=tile_name
-                ))
-                features.append(dict(
-                    tile=tile_name,
-                    geom=polygon
-                ))
+        grid_mapping, proj4 = create_grids(names=[name], projection=projection, meridian=meridian,
+                                           degreesx=[degreesx], degreesy=[degreesy],
+                                           bbox=bbox, srid=srid)
+        grid = grid_mapping[name]
 
         with db.session.begin_nested():
-            crs = CRS.from_proj4(tile_srs_p4)
+            crs = CRS.from_proj4(proj4)
             data = dict(
                 auth_name='Albers Equal Area',
                 auth_srid=srid,
                 srid=srid,
                 srtext=crs.to_wkt(),
-                proj4text=tile_srs_p4
+                proj4text=proj4
             )
 
             spatial_index, _ = get_or_create_model(SpatialRefSys, defaults=data, srid=srid)
 
-            grs = GridRefSys.create_geometry_table(table_name=name, features=features, srid=SRID_ALBERS_EQUAL_AREA)
+            grs = GridRefSys.create_geometry_table(table_name=name,
+                                                   features=grid['features'],
+                                                   srid=srid)
             grs.description = description
             db.session.add(grs)
 
-            [db.session.add(Tile(**tile, grs=grs)) for tile in tiles]
+            [db.session.add(Tile(**tile, grs=grs)) for tile in grid['tiles']]
         db.session.commit()
 
         return 'Grid {} created with successfully'.format(name), 201

cube_builder/forms.py

@@ -36,7 +36,8 @@ class GridRefSysForm(SQLAlchemyAutoSchema):
     meridian = fields.Integer(required=True, load_only=True)
     degreesx = fields.Float(required=True, load_only=True)
     degreesy = fields.Float(required=True, load_only=True)
-    bbox = fields.String(required=True, load_only=True)
+    bbox = fields.List(fields.Float, required=True, load_only=True)
+    srid = fields.Integer(required=True, load_only=True)
 
     class Meta:
         """Internal meta information of form interface."""

cube_builder/grids.py

@@ -0,0 +1,159 @@
+#
+# This file is part of Python Module for Cube Builder.
+# Copyright (C) 2019-2021 INPE.
+#
+# Cube Builder is free software; you can redistribute it and/or modify it
+# under the terms of the MIT License; see LICENSE file for more details.
+#
+
+"""Module for Grid creation."""
+
+import math
+from typing import List, Tuple
+
+from bdc_catalog.utils import geom_to_wkb
+from rasterio.crs import CRS
+from rasterio.warp import transform
+from shapely.geometry import Polygon
+
+
+def _create_tiles(tile_size: Tuple[float, float],
+                  grid_x_min: float, grid_y_max: float,
+                  bbox: Tuple[float, float, float, float],
+                  srid: float):
+    tile_size_x, tile_size_y = tile_size
+
+    tiles = []
+    features = []
+    xmin, xmax, ymin, ymax = bbox
+
+    h_min = int((xmin - grid_x_min) / tile_size_x)
+    h_max = int((xmax - grid_x_min) / tile_size_x)
+    v_min = int((grid_y_max - ymax) / tile_size_y)
+    v_max = int((grid_y_max - ymin) / tile_size_y)
+
+    for ix in range(h_min, h_max + 1):
+        x1 = grid_x_min + ix * tile_size_x
+        x2 = x1 + tile_size_x
+        for iy in range(v_min, v_max + 1):
+            y1 = grid_y_max - iy * tile_size_y
+            y2 = y1 - tile_size_y
+
+            polygon = geom_to_wkb(
+                Polygon(
+                    [
+                        (x1, y2),
+                        (x2, y2),
+                        (x2, y1),
+                        (x1, y1),
+                        (x1, y2)
+                    ]
+                ),
+                srid=srid
+            )
+
+            # Insert tile
+            tile_name = '{0:03d}{1:03d}'.format(ix, iy)
+            tiles.append(dict(
+                name=tile_name
+            ))
+            features.append(dict(
+                tile=tile_name,
+                geom=polygon
+            ))
+
+    return tiles, features
+
+
+def create_grids(names: List[str], projection, meridian,
+                 degreesx: List[float], degreesy: List[float], bbox, srid=100001):
+    """Create a list of hierarchically grids for Brazil Data Cube.
+
+    With a meridian, bounding box and list of degrees, this function creates a list of grids
+    based in the lowest degree. After that, the other grids are generated inherit the lowest
+    to be hierarchical.
+
+    For example, the grids BDC_SM_V2, BDC_MD_V2 and BDC_LG_V2 are hierarchically identical.
+        - BDC_SM_V2 consists in tiles of 1.5 by 1.0 degrees.
+        - BDC_MD_V2 represents 4 tiles of BDC_SM_V2 (2x2).
+        - BDC_LG_V2 represents 16 tiles of BDC_SM_V2 (4x4) or 4 tiles of BDC_MD_V2 (2x2)
+
+    Note:
+        Keep the order of names according degree's
+
+    Args:
+        names (List[str]): List of grid names
+        projection (str): The Grid projection kind - "sinu" for sinusoidal grids and "aea" for Albers Equal Area.
+        meridian (float): The center pixel for grid as reference.
+        degreesx (List[float]): List the degree in X
+        degreesy (List[float]): List the degree in Y
+        bbox (Tuple[float, float, float, float]): The bounding box limits (minx, miny, maxx, maxy).
+        srid (int): The projection SRID to create.
+    """
+    bbox_obj = {
+        "w": float(bbox[0]),
+        "n": float(bbox[1]),
+        "e": float(bbox[2]),
+        "s": float(bbox[3])
+    }
+    tile_srs_p4 = "+proj=longlat +ellps=GRS80 +no_defs"
+    if projection == 'aea':
+        tile_srs_p4 = "+proj=aea +lat_0=-12 +lon_0={} +lat_1=-2 +lat_2=-22 +x_0=5000000 +y_0=10000000 +ellps=GRS80 +units=m +no_defs".format(
+            meridian)
+    elif projection == 'sinu':
+        tile_srs_p4 = "+proj=sinu +lon_0={} +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs".format(
+            meridian)
+
+    ref_degree_x = degreesx[0]
+    ref_degree_y = degreesy[0]
+
+    # Tile size in meters (dx,dy) at center of system (argsmeridian,0.)
+    src_crs = '+proj=longlat +ellps=GRS80 +no_defs'
+    dst_crs = tile_srs_p4
+    xs = [(meridian - ref_degree_x / 2.), (meridian + ref_degree_x / 2.), meridian, meridian, 0.]
+    ys = [0., 0., -ref_degree_y / 2., ref_degree_y / 2., 0.]
+    out = transform(CRS.from_proj4(src_crs), CRS.from_proj4(dst_crs), xs, ys, zs=None)
+    xmin_center_tile = out[0][0]
+    xmax_center_tile = out[0][1]
+    ymin_center_tile = out[1][2]
+    ymax_center_tile = out[1][3]
+    tile_size_x = xmax_center_tile - xmin_center_tile
+    tile_size_y = ymax_center_tile - ymin_center_tile
+
+    bbox_alber_envelope_xs = [bbox_obj['w'], bbox_obj['e'], meridian, meridian]
+    bbox_alber_envelope_ys = [0., 0., bbox_obj['n'], bbox_obj['s']]
+    bbox_alber_envelope = transform(src_crs, dst_crs, bbox_alber_envelope_xs, bbox_alber_envelope_ys, zs=None)
+
+    total_tiles_left = math.ceil(abs((xmin_center_tile - bbox_alber_envelope[0][0])) / tile_size_x)
+    total_tiles_upper = math.ceil(abs((ymax_center_tile - bbox_alber_envelope[1][2])) / tile_size_y)
+
+    # Border coordinates of WRS grid
+    x_min = xmin_center_tile - (tile_size_x * total_tiles_left)
+    y_max = ymax_center_tile + (tile_size_y * total_tiles_upper)
+
+    # Upper Left is (xl,yu) Bottom Right is (xr,yb)
+    xs = [bbox_obj['w'], bbox_obj['e'], meridian, meridian]
+    ys = [0., 0., bbox_obj['n'], bbox_obj['s']]
+    out = transform(src_crs, dst_crs, xs, ys, zs=None)
+    xl = out[0][0]
+    xr = out[0][1]
+    yu = out[1][2]
+    yb = out[1][3]
+
+    grids = {}
+
+    for grid_name, res_x, res_y in zip(names, degreesx, degreesy):
+        factor_x = res_x / ref_degree_x
+        factor_y = res_y / ref_degree_y
+        grid_tile_size_x = tile_size_x * factor_x
+        grid_tile_size_y = tile_size_y * factor_y
+        tiles, features = _create_tiles(tile_size=(grid_tile_size_x, grid_tile_size_y),
+                                        grid_x_min=x_min, grid_y_max=y_max,
+                                        bbox=(xl, xr, yb, yu),
+                                        srid=srid)
+        grids[grid_name] = dict(
+            tiles=tiles,
+            features=features
+        )
+
+    return grids, tile_srs_p4

cube_builder/maestro.py

@@ -198,7 +198,10 @@ def orchestrate(self):
         dstart = self.params['start_date']
         dend = self.params['end_date']
 
-        timeline = Timeline(**temporal_schema, start_date=dstart, end_date=dend).mount()
+        if self.datacube.composite_function.alias == 'IDT':
+            timeline = [[dstart, dend]]
+        else:
+            timeline = Timeline(**temporal_schema, start_date=dstart, end_date=dend).mount()
 
         where = [Tile.grid_ref_sys_id == self.datacube.grid_ref_sys_id]
 
@@ -316,23 +319,6 @@ def datacube_bands(self) -> List[Band]:
             return list(filter(lambda band: band.name in self.params['bands'], self.bands))
         return [b for b in self.bands if b.name != DATASOURCE_NAME]
 
-    @staticmethod
-    def get_bbox(tile_id: str, grs: GridRefSys) -> str:
-        """Retrieve the bounding box representation as string."""
-        geom_table = grs.geom_table
-
-        bbox_result = db.session.query(
-            geom_table.c.tile,
-            func.ST_AsText(func.ST_BoundingDiagonal(func.ST_Transform(geom_table.c.geom, 4326)))
-        ).filter(
-            geom_table.c.tile == tile_id
-        ).first()
-
-        bbox = bbox_result[1][bbox_result[1].find('(') + 1:bbox_result[0].find(')')]
-        bbox = bbox.replace(' ', ',')
-
-        return bbox
-
     def dispatch_celery(self):
         """Dispatch datacube generation on celery workers.
 

cube_builder/utils/image.py

@@ -98,6 +98,7 @@ def validate_merges(images: ResultProxy, num_threads: int = Config.MAX_THREADS_I
         futures = executor.map(validate, images)
 
         output = dict()
+        error_map = dict()
 
         for row, errors in futures:
             if row is None:
@@ -112,9 +113,10 @@ def validate_merges(images: ResultProxy, num_threads: int = Config.MAX_THREADS_I
 
             output[row.date]['file'] = row.file
             output[row.date]['errors'].extend(errors)
-            if row.traceback:
+            if row.traceback and row.date not in error_map:
                 output[row.date]['errors'].append(dict(message=row.traceback, band=row.band,
                                                        filename=_file_name(row.link)))
+                error_map[row.date] = True
 
             output[row.date]['bands'].setdefault(row.band, list())
             output[row.date]['bands'][row.band].append(row.link)

cube_builder/utils/processing.py

@@ -374,12 +374,6 @@ def merge(merge_file: str, mask: dict, assets: List[dict], band: str, band_map:
     merge_file = Path(merge_file)
     merge_file.parent.mkdir(parents=True, exist_ok=True)
 
-    template.update({
-        'tiled': True,
-        "interleave": "pixel",
-        'driver': 'GTiff'
-    })
-
     options = dict(
         file=str(merge_file),
         efficacy=efficacy,