class GeometryTypeError(TypeError):
    """Raised when geometry type does not fit.""" class MapcheteIOError(IOError):
    """Raised when mapchete cannot read a file.""" MAPCHETE_IO_RETRY_SETTINGS = {
    "tries": int(os.environ.get("MAPCHETE_IO_RETRY_TRIES", 3)),
    "delay": float(os.environ.get("MAPCHETE_IO_RETRY_DELAY", 1)),
    "backoff": float(os.environ.get("MAPCHETE_IO_RETRY_BACKOFF", 1)),
} def fs_from_path(path, timeout=5, session=None, username=None, password=None, **kwargs):
    """Guess fsspec FileSystem from path and initialize using the desired options."""
    if path.startswith("s3://"):
        return fsspec.filesystem(
            "s3",
            requester_pays=os.environ.get("AWS_REQUEST_PAYER") == "requester",
            config_kwargs=dict(connect_timeout=timeout, read_timeout=timeout),
            session=session,
            client_kwargs=kwargs,
        )
    elif path.startswith(("http://", "https://")):
        if username:  # pragma: no cover
            from aiohttp import BasicAuth

            auth = BasicAuth(username, password)
        else:
            auth = None
        return fsspec.filesystem("https", auth=auth, asynchronous=False, **kwargs)
    else:
        return fsspec.filesystem("file", asynchronous=False, **kwargs) def path_exists(path, fs=None, **kwargs):
    """
    Check if file exists either remote or local.

    Parameters
    ----------
    path : path to file

    Returns
    -------
    exists : bool
    """
    fs = fs or fs_from_path(path, **kwargs)
    fs.invalidate_cache(path=path)
    return fs.exists(path) def makedirs(path):
    """
    Silently create all subdirectories of path if path is local.

    Parameters
    ----------
    path : path
    """
    fs = fs_from_path(path)
    # create parent directories on local filesystems
    if fs.protocol == "file":
        fs.makedirs(path, exist_ok=True) def copy(src_path, dst_path, src_fs=None, dst_fs=None, overwrite=False):
    """Copy path from one place to the other."""
    src_fs = src_fs or fs_from_path(src_path)
    dst_fs = dst_fs or fs_from_path(dst_path)

    if not overwrite and dst_fs.exists(dst_path):
        raise IOError(f"{dst_path} already exists")

    # create parent directories on local filesystems
    makedirs(os.path.dirname(dst_path))

    # copy either within a filesystem or between filesystems
    if src_fs == dst_fs:
        src_fs.copy(src_path, dst_path)
    else:
        with src_fs.open(src_path, "rb") as src:
            with dst_fs.open(dst_path, "wb") as dst:
                dst.write(src.read()) def reproject_geometry(
    geometry,
    src_crs=None,
    dst_crs=None,
    clip_to_crs_bounds=True,
    error_on_clip=False,
    segmentize_on_clip=False,
    segmentize=False,
    segmentize_fraction=100,
    validity_check=True,
    antimeridian_cutting=False,
    retry_with_clip=True,
):
    """
    Reproject a geometry to target CRS.

    Also, clips geometry if it lies outside the destination CRS boundary.
    Supported destination CRSes for clipping: 4326 (WGS84), 3857 (Spherical
    Mercator) and 3035 (ETRS89 / ETRS-LAEA).

    Parameters
    ----------
    geometry : ``shapely.geometry``
    src_crs : ``rasterio.crs.CRS`` or EPSG code
        CRS of source data
    dst_crs : ``rasterio.crs.CRS`` or EPSG code
        target CRS
    error_on_clip : bool
        raises a ``RuntimeError`` if a geometry is outside of CRS bounds
        (default: False)
    validity_check : bool
        checks if reprojected geometry is valid and throws ``TopologicalError``
        if invalid (default: True)
    clip_to_crs_bounds : bool
        Always clip geometries to CRS bounds. (default: True)
    antimeridian_cutting : bool
        cut geometry at Antimeridian; can result in a multipart output geometry

    Returns
    -------
    geometry : ``shapely.geometry``
    """
    logger.debug("reproject geometry from %s to %s", src_crs, dst_crs)
    src_crs = validate_crs(src_crs)
    dst_crs = validate_crs(dst_crs)

    def _reproject_geom(geometry, src_crs, dst_crs):
        if geometry.is_empty:
            return geometry
        else:
            out_geom = to_shape(
                transform_geom(
                    src_crs.to_dict(),
                    dst_crs.to_dict(),
                    mapping(geometry),
                    antimeridian_cutting=antimeridian_cutting,
                )
            )
            return _repair(out_geom) if validity_check else out_geom

    def _segmentize_value(geometry, segmentize_fraction):
        height = geometry.bounds[3] - geometry.bounds[1]
        width = geometry.bounds[2] - geometry.bounds[0]
        return min([height, width]) / segmentize_fraction

    # return repaired geometry if no reprojection needed
    if src_crs == dst_crs or geometry.is_empty:
        return _repair(geometry)
    # geometry needs to be clipped to its CRS bounds
    elif (
        clip_to_crs_bounds
        and dst_crs.is_epsg_code
        and dst_crs.get("init")
        != "epsg:4326"  # and is not WGS84 (does not need clipping)
        and (
            dst_crs.get("init") in CRS_BOUNDS
            or pyproj.CRS(dst_crs.to_epsg()).area_of_use.bounds
        )
    ):
        wgs84_crs = CRS().from_epsg(4326)
        # get dst_crs boundaries
        crs_bbox = box(
            *CRS_BOUNDS.get(
                dst_crs.get("init"), pyproj.CRS(dst_crs.to_epsg()).area_of_use.bounds
            )
        )
        # reproject geometry to WGS84
        geometry_4326 = _reproject_geom(geometry, src_crs, wgs84_crs)
        # raise error if geometry has to be clipped
        if error_on_clip and not geometry_4326.within(crs_bbox):
            raise RuntimeError("geometry outside target CRS bounds")

        clipped = crs_bbox.intersection(geometry_4326)

        # segmentize clipped geometry using one 100th of with or height depending on
        # which is shorter
        if segmentize_on_clip or segmentize:
            clipped = segmentize_geometry(
                clipped, _segmentize_value(clipped, segmentize_fraction)
            )

        # clip geometry dst_crs boundaries and return
        return _reproject_geom(clipped, wgs84_crs, dst_crs)

    # return without clipping if destination CRS does not have defined bounds
    else:
        try:
            if segmentize:
                return _reproject_geom(
                    segmentize_geometry(
                        geometry, _segmentize_value(geometry, segmentize_fraction)
                    ),
                    src_crs,
                    dst_crs,
                )
            else:
                return _reproject_geom(geometry, src_crs, dst_crs)
        except TopologicalError:  # pragma: no cover
            raise
        except ValueError as exc:  # pragma: no cover
            if retry_with_clip:
                logger.error(
                    "error when transforming %s from %s to %s: %s, trying to use CRS bounds clip",
                    geometry,
                    src_crs,
                    dst_crs,
                    exc,
                )
                return reproject_geometry(
                    geometry,
                    src_crs=src_crs,
                    dst_crs=dst_crs,
                    clip_to_crs_bounds=True,
                    error_on_clip=error_on_clip,
                    segmentize_on_clip=segmentize_on_clip,
                    segmentize=segmentize,
                    segmentize_fraction=segmentize_fraction,
                    validity_check=validity_check,
                    antimeridian_cutting=antimeridian_cutting,
                    retry_with_clip=False,
                )
            else:
                raise def segmentize_geometry(geometry, segmentize_value):
    """
    Segmentize Polygon outer ring by segmentize value.

    Just Polygon geometry type supported.

    Parameters
    ----------
    geometry : ``shapely.geometry``
    segmentize_value: float

    Returns
    -------
    geometry : ``shapely.geometry``
    """
    if geometry.geom_type != "Polygon":
        raise TypeError("segmentize geometry type must be Polygon")

    return Polygon(
        LinearRing(
            [
                p
                # pick polygon linestrings
                for l in map(
                    lambda x: LineString([x[0], x[1]]),
                    zip(geometry.exterior.coords[:-1], geometry.exterior.coords[1:]),
                )
                # interpolate additional points in between and don't forget end point
                for p in [
                    l.interpolate(segmentize_value * i).coords[0]
                    for i in range(int(l.length / segmentize_value))
                ]
                + [l.coords[1]]
            ]
        )
    ) def to_shape(geom):
    """
    Convert geometry to shapely geometry if necessary.

    Parameters
    ----------
    geom : shapely geometry or GeoJSON mapping

    Returns
    -------
    shapely geometry
    """
    return shape(geom) if isinstance(geom, dict) else geom def multipart_to_singleparts(geom):
    """
    Yield single part geometries if geom is multipart, otherwise yield geom.

    Parameters
    ----------
    geom : shapely geometry

    Returns
    -------
    shapely single part geometries
    """
    if isinstance(geom, base.BaseGeometry):
        if hasattr(geom, "geoms"):
            for subgeom in geom.geoms:
                yield subgeom
        else:
            yield geom def clean_geometry_type(
    geometry, target_type, allow_multipart=True, raise_exception=True
):
    """
    Return geometry of a specific type if possible.

    Filters and splits up GeometryCollection into target types. This is
    necessary when after clipping and/or reprojecting the geometry types from
    source geometries change (i.e. a Polygon becomes a LineString or a
    LineString becomes Point) in some edge cases.

    Parameters
    ----------
    geometry : ``shapely.geometry``
    target_type : string
        target geometry type
    allow_multipart : bool
        allow multipart geometries (default: True)

    Returns
    -------
    cleaned geometry : ``shapely.geometry``
        returns None if input geometry type differs from target type

    Raises
    ------
    GeometryTypeError : if geometry type does not match target_type
    """
    if target_type == "Unknown":  # pragma: no cover
        return geometry

    multipart_geoms = {
        "Point": MultiPoint,
        "LineString": MultiLineString,
        "Polygon": MultiPolygon,
        "MultiPoint": MultiPoint,
        "MultiLineString": MultiLineString,
        "MultiPolygon": MultiPolygon,
    }
    if target_type not in multipart_geoms.keys():
        raise TypeError("target type is not supported: %s" % target_type)

    if geometry.geom_type == target_type:
        return geometry

    elif allow_multipart:
        target_multipart_type = multipart_geoms[target_type]
        if geometry.geom_type == "GeometryCollection":
            return target_multipart_type(
                [
                    clean_geometry_type(
                        g, target_type, allow_multipart, raise_exception=raise_exception
                    )
                    for g in geometry.geoms
                ]
            )
        elif (
            isinstance(geometry, target_multipart_type)
            or multipart_geoms[geometry.geom_type] == target_multipart_type
        ):
            return geometry

    if raise_exception:
        raise GeometryTypeError(
            "geometry type does not match: %s, %s" % (geometry.geom_type, target_type)
        )

    else:
        return GeometryCollection() def _repair(geom):
    repaired = geom.buffer(0) if geom.geom_type in ["Polygon", "MultiPolygon"] else geom
    if repaired.is_valid:
        return repaired
    else:
        raise TopologicalError(
            "geometry is invalid (%s) and cannot be repaired"
            % explain_validity(repaired)
        ) def validate_bounds(bounds):
    """
    Return validated bounds.

    Bounds must be a list or tuple with exactly four elements.

    Parameters
    ----------
    bounds : list or tuple

    Returns
    -------
    Bounds

    Raises
    ------
    TypeError if type is invalid.
    """
    if not isinstance(bounds, (tuple, list)):
        raise TypeError("bounds must be either a tuple or a list: %s" % str(bounds))
    if len(bounds) != 4:
        raise ValueError("bounds has to have exactly four values: %s" % str(bounds))
    return Bounds(*bounds)
from contextlib import ExitStack
from fiona.errors import DriverError, FionaError, FionaValueError
from fiona.io import MemoryFile
from retry import retry
from rasterio.crs import CRS
from shapely.geometry import box, mapping, shape
from shapely.errors import TopologicalError
from tilematrix import clip_geometry_to_srs_bounds
from itertools import chain
from mapchete.errors import GeometryTypeError, MapcheteIOError
from mapchete.io._misc import MAPCHETE_IO_RETRY_SETTINGS
from mapchete.io._path import fs_from_path, path_exists, makedirs, copy
from mapchete.io._geometry_operations import (
    reproject_geometry,
    segmentize_geometry,
    to_shape,
    multipart_to_singleparts,
    clean_geometry_type,
    _repair,
)
from mapchete.validate import validate_bounds
                from rtree import index
import os
import logging
import fiona
import json
import warnings
            import geobuf
            )
            if allow_multipart_geometries:
                cleaned_output_fetures = [clipped]
            else:
                cleaned_output_fetures = multipart_to_singleparts(clipped)
            for out_geom in cleaned_output_fetures:
                if out_geom.is_empty:  # pragma: no cover
                    continue

                out_features.append(
                    {"geometry": mapping(out_geom), "properties": feature["properties"]}
                )
        except Exception as e:
            logger.warning("failed to prepare geometry for writing: %s", e)
            continue

    # write if there are output features
    if out_features:
        try:
            if out_driver.lower() in ["geobuf"]:
                # write data to remote file
                with VectorWindowMemoryFile(
                    tile=out_tile,
                    features=out_features,
                    schema=out_schema,
                    driver=out_driver,
                ) as memfile:
                    logger.debug((out_tile.id, "write tile", out_path))
                    with fs_from_path(out_path).open(out_path, "wb") as dst:
                        dst.write(memfile)
            else:  # pragma: no cover
                # write data to local file
                # this part is not covered by tests as we now try to let fiona directly
                # write to S3
                with fiona.open(
                    out_path,
                    "w",
                    schema=out_schema,
                    driver=out_driver,
                    crs=out_tile.crs.to_dict(),
                ) as dst:
                    logger.debug((out_tile.id, "write tile", out_path))
                    dst.writerecords(out_features)
        except Exception as e:
            logger.error("error while writing file %s: %s", out_path, e)
            raise

    else:
        logger.debug((out_tile.id, "nothing to write", out_path))


class VectorWindowMemoryFile:
    """Context manager around fiona.io.MemoryFile."""

    def __init__(self, tile=None, features=None, schema=None, driver=None):
        """Prepare data & profile."""
        self.tile = tile
        self.schema = schema
        self.driver = driver
        self.features = features

    def __enter__(self):
        """Open MemoryFile, write data and return."""
        if self.driver.lower() == "geobuf":

            return geobuf.encode(
                dict(
                    type="FeatureCollection",
                    features=[dict(f, type="Feature") for f in self.features],
                )
            )
        else:  # pragma: no cover
            # this part is excluded now for tests as we try to let fiona write directly
            # to S3
            self.fio_memfile = MemoryFile()
            with self.fio_memfile.open(
                schema=self.schema, driver=self.driver, crs=self.tile.crs
            ) as dst:
                dst.writerecords(self.features)
            return self.fio_memfile.getbuffer()

    def __exit__(self, *args):
        """Make sure MemoryFile is closed."""
        try:
            self.fio_memfile.close()
        except AttributeError:
            pass


@retry(
    logger=logger,
    exceptions=(DriverError, FionaError, FionaValueError),
    **MAPCHETE_IO_RETRY_SETTINGS,
)
def _get_reprojected_features(
    inp=None,
    dst_bounds=None,
    dst_crs=None,
    validity_check=False,
    clip_to_crs_bounds=False,
):
    logger.debug("reading %s", inp)
    with ExitStack() as exit_stack:
        if isinstance(inp, str):
            try:
                src = exit_stack.enter_context(fiona.open(inp, "r"))
                src_crs = CRS(src.crs)
            except Exception as e:
                if path_exists(inp):
                    logger.error("error while reading file %s: %s", inp, e)
                    raise e
                else:
                    raise FileNotFoundError("%s not found" % inp)
        else:
            src = inp
            src_crs = inp.crs
        # reproject tile bounding box to source file CRS for filter
        if src_crs == dst_crs:
            dst_bbox = box(*dst_bounds)
        else:

            dst_bbox = reproject_geometry(