Merge 8c4dffc into b54e323

CHANGELOG.rst

@@ -2,6 +2,13 @@
 Changelog
 #########
 
+----
+0.27
+----
+* enable reading from output tile directories which have a different CRS
+* enable GeoPackage as single file input
+* fixed antimeridian shift check
+
 ----
 0.26
 ----

mapchete/__init__.py

@@ -14,4 +14,4 @@
 logging.getLogger("fiona").setLevel(logging.ERROR)
 
 
-__version__ = "0.26"
+__version__ = "0.27"

mapchete/formats/default/tile_directory.py

@@ -5,12 +5,12 @@
 import numpy as np
 import numpy.ma as ma
 import os
-from shapely.geometry import box
+from shapely.geometry import box, shape, mapping
 
 from mapchete.config import validate_values
 from mapchete.errors import MapcheteConfigError
 from mapchete.formats import base, load_output_writer
-from mapchete.io import path_exists, absolute_path, read_json
+from mapchete.io import path_exists, absolute_path, read_json, tile_to_zoom_level
 from mapchete.io.vector import reproject_geometry, read_vector_window
 from mapchete.io.raster import read_raster_window, create_mosaic, resample_from_array
 from mapchete.tile import BufferedTilePyramid
@@ -143,34 +143,93 @@ def __init__(self, input_params, **kwargs):
         else:
             self._profile = None
 
-    def open(self, tile, **kwargs):
+    def open(
+        self,
+        tile,
+        tile_directory_zoom=None,
+        matching_method="gdal",
+        matching_max_zoom=None,
+        matching_precision=8,
+        fallback_to_higher_zoom=False,
+        **kwargs
+    ):
         """
         Return InputTile object.
 
         Parameters
         ----------
         tile : ``Tile``
+        tile_directory_zoom : None
+            If set, data will be read from exactly this zoom level
+        matching_method : str ('gdal' or 'min') (default: 'gdal')
+            gdal: Uses GDAL's standard method. Here, the target resolution is calculated
+                by averaging the extent's pixel sizes over both x and y axes. This
+                approach returns a zoom level which may not have the best quality but will
+                speed up reading significantly.
+            min: Returns the zoom level which matches the minimum resolution of the
+                extents four corner pixels. This approach returns the zoom level with the
+                best possible quality but with low performance. If the tile extent is
+                outside of the destination pyramid, a TopologicalError will be raised.
+        matching_max_zoom : int (default: None)
+            If set, it will prevent reading from zoom levels above the maximum.
+        matching_precision : int
+            Round resolutions to n digits before comparing.
+        fallback_to_higher_zoom : bool (default: False)
+            In case no data is found at zoom level, try to read data from higher zoom
+            levels. Enabling this setting can lead to many IO requests in areas with no
+            data.
 
         Returns
         -------
         input tile : ``InputTile``
             tile view of input data
         """
+        # determine tile bounds in TileDirectory CRS
+        td_bounds = reproject_geometry(
+            tile.bbox,
+            src_crs=tile.tp.crs,
+            dst_crs=self.td_pyramid.crs
+        ).bounds
+
+        # find target zoom level
+        if tile_directory_zoom is not None:
+            zoom = tile_directory_zoom
+        else:
+            zoom = tile_to_zoom_level(
+                tile, dst_pyramid=self.td_pyramid, matching_method=matching_method,
+                precision=matching_precision
+            )
+            if matching_max_zoom is not None:
+                zoom = min([zoom, matching_max_zoom])
+
+        if fallback_to_higher_zoom:
+            tiles_paths = []
+            # check if tiles exist otherwise try higher zoom level
+            while len(tiles_paths) == 0 and zoom >= 0:
+                tiles_paths = _get_tiles_paths(
+                    basepath=self.path,
+                    ext=self._ext,
+                    pyramid=self.td_pyramid,
+                    bounds=td_bounds,
+                    zoom=zoom
+                )
+                logger.debug("%s existing tiles found at zoom %s", len(tiles_paths), zoom)
+                zoom -= 1
+        else:
+            tiles_paths = _get_tiles_paths(
+                basepath=self.path,
+                ext=self._ext,
+                pyramid=self.td_pyramid,
+                bounds=td_bounds,
+                zoom=zoom
+            )
+            logger.debug("%s existing tiles found at zoom %s", len(tiles_paths), zoom)
         return InputTile(
             tile,
-            tiles_paths=[
-                (_tile, _path)
-                for _tile, _path in [
-                    (t, os.path.join(*([
-                        self.path, str(t.zoom), str(t.row), str(t.col)
-                        ])) + "." + self._ext)
-                    for t in self.td_pyramid.tiles_from_bounds(
-                        tile.bounds, tile.zoom)
-                ]
-                if path_exists(_path)
-            ],
+            tiles_paths=tiles_paths,
             file_type=self._file_type,
             profile=self._profile,
+            td_crs=self.td_pyramid.crs,
             **kwargs
         )
 
@@ -195,6 +254,22 @@ def bbox(self, out_crs=None):
         )
 
 
+def _get_tiles_paths(basepath=None, ext=None, pyramid=None, bounds=None, zoom=None):
+    return [
+        (_tile, _path)
+        for _tile, _path in [
+            (
+                t,
+                "%s.%s" % (
+                    os.path.join(*([basepath, str(t.zoom), str(t.row), str(t.col)])), ext
+                )
+            )
+            for t in pyramid.tiles_from_bounds(bounds, zoom)
+        ]
+        if path_exists(_path)
+    ]
+
+
 class InputTile(base.InputTile):
     """
     Target Tile representation of input data.
@@ -216,10 +291,15 @@ def __init__(self, tile, **kwargs):
         self._tiles_paths = kwargs["tiles_paths"]
         self._file_type = kwargs["file_type"]
         self._profile = kwargs["profile"]
+        self._td_crs = kwargs["td_crs"]
 
     def read(
-        self, validity_check=False, indexes=None, resampling="nearest",
-        dst_nodata=None, gdal_opts=None
+        self,
+        validity_check=False,
+        indexes=None,
+        resampling="nearest",
+        dst_nodata=None,
+        gdal_opts=None
     ):
         """
         Read reprojected & resampled input data.
@@ -244,40 +324,37 @@ def read(
         -------
         data : list for vector files or numpy array for raster files
         """
+        logger.debug("reading data from CRS %s to CRS %s", self._td_crs, self.tile.tp.crs)
         if self._file_type == "vector":
             if self.is_empty():
                 return []
-            return list(chain.from_iterable([
-                read_vector_window(
-                    _path, self.tile, validity_check=validity_check)
-                for _, _path in self._tiles_paths
-            ]))
+            else:
+                return read_vector_window(
+                    [path for _, path in self._tiles_paths],
+                    self.tile,
+                    validity_check=validity_check
+                )
         else:
             if self.is_empty():
-                count = (len(indexes) if indexes else self._profile["count"], )
+                bands = len(indexes) if indexes else self._profile["count"]
                 return ma.masked_array(
                     data=np.full(
-                        count + self.tile.shape, self._profile["nodata"],
-                        dtype=self._profile["dtype"]),
+                        (bands, self.tile.height, self.tile.width),
+                        self._profile["nodata"],
+                        dtype=self._profile["dtype"]
+                    ),
                     mask=True
                 )
-            tiles = [
-                (
-                    _tile,
-                    read_raster_window(
-                        _path, _tile, indexes=indexes, resampling=resampling,
-                        src_nodata=self._profile["nodata"], dst_nodata=dst_nodata,
-                        gdal_opts=gdal_opts
-                    )
+            else:
+                return read_raster_window(
+                    [path for _, path in self._tiles_paths],
+                    self.tile,
+                    indexes=indexes,
+                    resampling=resampling,
+                    src_nodata=self._profile["nodata"],
+                    dst_nodata=dst_nodata,
+                    gdal_opts=gdal_opts
                 )
-                for _tile, _path in self._tiles_paths
-            ]
-            return resample_from_array(
-                in_raster=create_mosaic(tiles=tiles, nodata=self._profile["nodata"]),
-                out_tile=self.tile,
-                resampling=resampling,
-                nodataval=self._profile["nodata"]
-            )
 
     def is_empty(self):
         """

mapchete/io/__init__.py

@@ -4,6 +4,8 @@
 import logging
 import os
 import rasterio
+from rasterio.warp import calculate_default_transform
+from shapely.errors import TopologicalError
 from shapely.geometry import box
 from tilematrix import TilePyramid
 from urllib.request import urlopen
@@ -33,8 +35,7 @@ def get_best_zoom_level(input_file, tile_pyramid_type):
     Parameters
     ----------
     input_file : path to raster file
-    tile_pyramid_type : ``TilePyramid`` projection (``geodetic`` or
-        ``mercator``)
+    tile_pyramid_type : ``TilePyramid`` projection (``geodetic`` or``mercator``)
 
     Returns
     -------
@@ -89,6 +90,91 @@ def get_segmentize_value(input_file=None, tile_pyramid=None):
     return pixelsize * tile_pyramid.tile_size
 
 
+def tile_to_zoom_level(tile, dst_pyramid=None, matching_method="gdal", precision=8):
+    """
+    Determine the best zoom level in target TilePyramid from given Tile.
+
+
+    Parameters
+    ----------
+    tile : BufferedTile
+    dst_pyramid : BufferedTilePyramid
+    matching_method : str ('gdal' or 'min')
+        gdal: Uses GDAL's standard method. Here, the target resolution is calculated by
+            averaging the extent's pixel sizes over both x and y axes. This approach
+            returns a zoom level which may not have the best quality but will speed up
+            reading significantly.
+        min: Returns the zoom level which matches the minimum resolution of the extent's
+            four corner pixels. This approach returns the zoom level with the best
+            possible quality but with low performance. If the tile extent is outside of
+            the destination pyramid, a TopologicalError will be raised.
+    precision : int
+        Round resolutions to n digits before comparing.
+
+    Returns
+    -------
+    zoom : int
+    """
+    def width_height(bounds):
+        try:
+            l, b, r, t = reproject_geometry(
+                box(*bounds), src_crs=tile.crs, dst_crs=dst_pyramid.crs
+            ).bounds
+        except ValueError:
+            raise TopologicalError("bounds cannot be translated into target CRS")
+        return r - l, t - b
+
+    if tile.tp.crs == dst_pyramid.crs:
+        return tile.zoom
+    else:
+        if matching_method == "gdal":
+            # use rasterio/GDAL method to calculate default warp target properties
+            transform, width, height = calculate_default_transform(
+                tile.tp.crs,
+                dst_pyramid.crs,
+                tile.width,
+                tile.height,
+                *tile.bounds
+            )
+            # this is the resolution the tile would have in destination TilePyramid CRS
+            tile_resolution = round(transform[0], precision)
+        elif matching_method == "min":
+            # calculate the minimum pixel size from the four tile corner pixels
+            l, b, r, t = tile.bounds
+            x = tile.pixel_x_size
+            y = tile.pixel_y_size
+            res = []
+            for bounds in [
+                (l, t - y, l + x, t),  # left top
+                (l, b, l + x, b + y),  # left bottom
+                (r - x, b, r, b + y),  # right bottom
+                (r - x, t - y, r, t)   # right top
+            ]:
+                try:
+                    w, h = width_height(bounds)
+                    res.extend([w, h])
+                except TopologicalError:
+                    logger.debug("pixel outside of destination pyramid")
+            if res:
+                tile_resolution = round(min(res), precision)
+            else:
+                raise TopologicalError("tile outside of destination pyramid")
+        else:
+            raise ValueError("invalid method given: %s", matching_method)
+        logger.debug(
+            "we are looking for a zoom level interpolating to %s resolution",
+            tile_resolution
+        )
+        zoom = 0
+        while True:
+            td_resolution = round(dst_pyramid.pixel_x_size(zoom), precision)
+            if td_resolution <= tile_resolution:
+                break
+            zoom += 1
+        logger.debug("target zoom for %s: %s (%s)", tile_resolution, zoom, td_resolution)
+        return zoom
+
+
 def path_is_remote(path, s3=True):
     """
     Determine whether file path is remote or local.