rasterio/tests/test_warp.py

"""rasterio.warp module tests"""

import json
import logging
import os
import sys

from affine import Affine
import numpy as np
from numpy.testing import assert_almost_equal
import pytest

import rasterio
from rasterio.control import GroundControlPoint
from rasterio.crs import CRS
from rasterio.enums import Resampling
from rasterio.env import GDALVersion
from rasterio.errors import (
    CRSError,
    GDALVersionError,
    TransformError,
    WarpOperationError,
)
from rasterio.warp import (
    reproject,
    transform_geom,
    transform,
    transform_bounds,
    calculate_default_transform,
    aligned_target,
    SUPPORTED_RESAMPLING,
    GDAL2_RESAMPLING,
)
from rasterio import windows

from . import rangehttpserver
from .conftest import requires_gdal22, requires_gdal3, requires_gdal_lt_3

log = logging.getLogger(__name__)

gdal_version = GDALVersion.runtime()


DST_TRANSFORM = Affine(300.0, 0.0, -8789636.708, 0.0, -300.0, 2943560.235)


def flatten_coords(coordinates):
    """Yield a flat sequence of coordinates to help testing"""
    for elem in coordinates:
        if isinstance(elem, (float, int)):
            yield elem

        else:
            for x in flatten_coords(elem):
                yield x


reproj_expected = (
    ({"CHECK_WITH_INVERT_PROJ": False}, 6644), ({"CHECK_WITH_INVERT_PROJ": True}, 6644)
)


class ReprojectParams:
    """Class to assist testing reprojection by encapsulating parameters."""

    def __init__(self, left, bottom, right, top, width, height, src_crs, dst_crs):
        self.width = width
        self.height = height
        src_res = float(right - left) / float(width)
        self.src_transform = Affine(src_res, 0, left, 0, -src_res, top)
        self.src_crs = src_crs
        self.dst_crs = dst_crs

        dt, dw, dh = calculate_default_transform(
            src_crs, dst_crs, width, height, left, bottom, right, top
        )
        self.dst_transform = dt
        self.dst_width = dw
        self.dst_height = dh


def default_reproject_params():
    return ReprojectParams(
        left=-120,
        bottom=30,
        right=-80,
        top=70,
        width=80,
        height=80,
        src_crs=CRS.from_epsg(4326),
        dst_crs=CRS.from_epsg(2163),
    )


def uninvertable_reproject_params():
    return ReprojectParams(
        left=-120,
        bottom=30,
        right=-80,
        top=70,
        width=80,
        height=80,
        src_crs=CRS.from_epsg(4326),
        dst_crs=CRS.from_epsg(26836),
    )


WGS84_crs = CRS.from_epsg(4326)


class RangeRequestErrorHandler(rangehttpserver.RangeRequestHandler):
    """Return 500 for a range of bytes to simulate a malfunctioning server.

    The byte range is specific to rasterio's RGB.byte.tif file.

    """

    def send_head(self):
        if "Range" not in self.headers:
            self.range = None
            return super().send_head()
        try:
            self.range = rangehttpserver.parse_byte_range(self.headers["Range"])
        except ValueError as e:
            self.send_error(400, "Invalid byte range")
            return None
        first, last = self.range

        # Our "poison byte" is at position 1609000.
        if first <= 1609000 <= last:
            self.send_error(500, "Boom!")
            return None

        # Mirroring SimpleHTTPServer.py here
        path = self.translate_path(self.path)
        f = None
        ctype = self.guess_type(path)
        try:
            f = open(path, "rb")
        except IOError:
            self.send_error(404, "File not found")
            return None

        fs = os.fstat(f.fileno())
        file_len = fs[6]
        if first >= file_len:
            self.send_error(416, "Requested Range Not Satisfiable")
            return None

        self.send_response(206)
        self.send_header("Content-type", ctype)
        self.send_header("Accept-Ranges", "bytes")

        if last is None or last >= file_len:
            last = file_len - 1
        response_length = last - first + 1

        self.send_header(
            "Content-Range", "bytes %s-%s/%s" % (first, last, file_len)
        )
        self.send_header("Content-Length", str(response_length))
        self.send_header("Last-Modified", self.date_time_string(fs.st_mtime))
        self.end_headers()
        return f


def test_transform_src_crs_none():
    with pytest.raises(CRSError):
        transform(None, WGS84_crs, [1], [1])


def test_transform_dst_crs_none():
    with pytest.raises(CRSError):
        transform(WGS84_crs, None, [1], [1])


def test_transform_bounds_src_crs_none():
    with pytest.raises(CRSError):
        transform_bounds(None, WGS84_crs, 0, 0, 0, 0)


def test_transform_bounds_dst_crs_none():
    with pytest.raises(CRSError):
        transform_bounds(WGS84_crs, None, 0, 0, 0, 0)


def test_transform_geom_src_crs_none():
    with pytest.raises(CRSError):
        transform_geom(None, WGS84_crs, None)


def test_transform_geom_dst_crs_none():
    with pytest.raises(CRSError):
        transform_geom(WGS84_crs, None, None)


def test_reproject_src_crs_none():
    with pytest.raises(CRSError):
        reproject(
            np.ones((2, 2)),
            np.zeros((2, 2)),
            src_transform=Affine.identity(),
            dst_transform=Affine.identity(),
            dst_crs=WGS84_crs,
        )


def test_reproject_dst_crs_none():
    with pytest.raises(CRSError):
        reproject(
            np.ones((2, 2)),
            np.zeros((2, 2)),
            src_transform=Affine.identity(),
            dst_transform=Affine.identity(),
            src_crs=WGS84_crs,
        )


def test_transform():
    """2D and 3D."""
    WGS84_crs = CRS.from_epsg(4326)
    WGS84_points = ([12.492269], [41.890169], [48.])
    ECEF_crs = CRS.from_epsg(4978)
    ECEF_points = ([4642610.], [1028584.], [4236562.])
    ECEF_result = transform(WGS84_crs, ECEF_crs, *WGS84_points)
    assert np.allclose(np.array(ECEF_result), np.array(ECEF_points))

    UTM33_crs = CRS.from_epsg(32633)
    UTM33_points = ([291952], [4640623])
    UTM33_result = transform(WGS84_crs, UTM33_crs, *WGS84_points[:2])
    assert np.allclose(np.array(UTM33_result), np.array(UTM33_points))


def test_transform_bounds():
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        l, b, r, t = src.bounds
        assert np.allclose(
            transform_bounds(src.crs, CRS.from_epsg(4326), l, b, r, t),
            (
                -78.95864996545055,
                23.564991210854686,
                -76.57492370013823,
                25.550873767433984,
            ),
        )


def test_transform_bounds__esri_wkt():
    left, bottom, right, top = \
        (-78.95864996545055, 23.564991210854686,
         -76.57492370013823, 25.550873767433984)
    dst_projection_string = (
        'PROJCS["USA_Contiguous_Albers_Equal_Area_Conic_USGS_version",'
        'GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",'
        'SPHEROID["GRS_1980",6378137.0,298.257222101]],'
        'PRIMEM["Greenwich",0.0],'
        'UNIT["Degree",0.0174532925199433]],'
        'PROJECTION["Albers"],'
        'PARAMETER["false_easting",0.0],'
        'PARAMETER["false_northing",0.0],'
        'PARAMETER["central_meridian",-96.0],'
        'PARAMETER["standard_parallel_1",29.5],'
        'PARAMETER["standard_parallel_2",45.5],'
        'PARAMETER["latitude_of_origin",23.0],'
        'UNIT["Meter",1.0],'
        'VERTCS["NAVD_1988",'
        'VDATUM["North_American_Vertical_Datum_1988"],'
        'PARAMETER["Vertical_Shift",0.0],'
        'PARAMETER["Direction",1.0],UNIT["Centimeter",0.01]]]')
    assert np.allclose(
        transform_bounds(CRS.from_epsg(4326),
                         dst_projection_string,
                         left,
                         bottom,
                         right,
                         top),
        (
            1721263.7931814701,
            219684.49332178483,
            2002926.56696663,
            479360.16562217404),
    )


@pytest.mark.parametrize(
    "density,expected",
    [
        (0, (-1684649.41338, -350356.81377, 1684649.41338, 2234551.18559)),
        (100, (-1684649.41338, -555777.79210, 1684649.41338, 2234551.18559)),
    ],
)
def test_transform_bounds_densify(density, expected):
    # This transform is non-linear along the edges, so densification produces
    # a different result than otherwise
    src_crs = CRS.from_epsg(4326)
    dst_crs = CRS.from_epsg(2163)
    with rasterio.Env(OSR_USE_NON_DEPRECATED="NO"):
        assert np.allclose(
            expected,
            transform_bounds(src_crs, dst_crs, -120, 40, -80, 64, densify_pts=density),
        )


def test_transform_bounds_no_change():
    """Make sure that going from and to the same crs causes no change."""
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        l, b, r, t = src.bounds
        assert np.allclose(transform_bounds(src.crs, src.crs, l, b, r, t), src.bounds)


def test_transform_bounds_densify_out_of_bounds():
    with pytest.raises(ValueError):
        transform_bounds(
            CRS.from_epsg(4326),
            CRS.from_epsg(32610),
            -120,
            40,
            -80,
            64,
            densify_pts=-10,
        )


def test_calculate_default_transform():
    target_transform = Affine(
        0.0028535715391804096,
        0.0,
        -78.95864996545055,
        0.0,
        -0.0028535715391804096,
        25.550873767433984,
    )

    with rasterio.open("tests/data/RGB.byte.tif") as src:
        wgs84_crs = CRS.from_epsg(4326)
        dst_transform, width, height = calculate_default_transform(
            src.crs, wgs84_crs, src.width, src.height, *src.bounds
        )

        assert dst_transform.almost_equals(target_transform)
        assert width == 835
        assert height == 696


def test_calculate_default_transform_single_resolution():
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        target_resolution = 0.1
        target_transform = Affine(
            target_resolution,
            0.0,
            -78.95864996545055,
            0.0,
            -target_resolution,
            25.550873767433984,
        )
        dst_transform, width, height = calculate_default_transform(
            src.crs,
            CRS.from_epsg(4326),
            src.width,
            src.height,
            *src.bounds,
            resolution=target_resolution
        )

        assert dst_transform.almost_equals(target_transform)
        assert width == 24
        assert height == 20


def test_calculate_default_transform_multiple_resolutions():
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        target_resolution = (0.2, 0.1)
        target_transform = Affine(
            target_resolution[0],
            0.0,
            -78.95864996545055,
            0.0,
            -target_resolution[1],
            25.550873767433984,
        )

        dst_transform, width, height = calculate_default_transform(
            src.crs,
            CRS.from_epsg(4326),
            src.width,
            src.height,
            *src.bounds,
            resolution=target_resolution
        )

        assert dst_transform.almost_equals(target_transform)
        assert width == 12
        assert height == 20


def test_calculate_default_transform_dimensions():
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        dst_width, dst_height = (113, 103)
        target_transform = Affine(
            0.02108612597535966,
            0.0,
            -78.95864996545055,
            0.0,
            -0.0192823863230055,
            25.550873767433984,
        )

        dst_transform, width, height = calculate_default_transform(
            src.crs,
            CRS.from_epsg(4326),
            src.width,
            src.height,
            *src.bounds,
            dst_width=dst_width,
            dst_height=dst_height
        )

        assert dst_transform.almost_equals(target_transform)
        assert width == dst_width
        assert height == dst_height


def test_reproject_ndarray():
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        source = src.read(1)

    dst_crs = dict(
        proj="merc",
        a=6378137,
        b=6378137,
        lat_ts=0.0,
        lon_0=0.0,
        x_0=0.0,
        y_0=0,
        k=1.0,
        units="m",
        nadgrids="@null",
        wktext=True,
        no_defs=True,
    )
    out = np.empty(src.shape, dtype=np.uint8)
    reproject(
        source,
        out,
        src_transform=src.transform,
        src_crs=src.crs,
        dst_transform=DST_TRANSFORM,
        dst_crs=dst_crs,
        resampling=Resampling.nearest,
    )
    assert (out > 0).sum() == 438113


def test_reproject_view():
    """Source views are reprojected properly"""
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        source = src.read(1)

    window = windows.Window(100, 100, 500, 500)
    # window = windows.get_data_window(source)
    reduced_array = source[window.toslices()]
    reduced_transform = windows.transform(window, src.transform)

    # Assert that we're working with a view.
    assert reduced_array.base is source

    dst_crs = dict(
        proj="merc",
        a=6378137,
        b=6378137,
        lat_ts=0.0,
        lon_0=0.0,
        x_0=0.0,
        y_0=0,
        k=1.0,
        units="m",
        nadgrids="@null",
        wktext=True,
        no_defs=True,
    )

    out = np.empty(src.shape, dtype=np.uint8)

    reproject(
        reduced_array,
        out,
        src_transform=reduced_transform,
        src_crs=src.crs,
        dst_transform=DST_TRANSFORM,
        dst_crs=dst_crs,
        resampling=Resampling.nearest,
    )

    assert (out > 0).sum() == 299199


def test_reproject_epsg():
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        source = src.read(1)

    dst_crs = {"init": "epsg:3857"}
    out = np.empty(src.shape, dtype=np.uint8)
    reproject(
        source,
        out,
        src_transform=src.transform,
        src_crs=src.crs,
        dst_transform=DST_TRANSFORM,
        dst_crs=dst_crs,
        resampling=Resampling.nearest,
    )
    assert (out > 0).sum() == 438113


def test_reproject_epsg__simple_array():
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        source = src.read(1)

    dst_crs = {"init": "EPSG:3857"}
    out, dst_transform = reproject(
        source,
        src_transform=src.transform,
        src_crs=src.crs,
        dst_crs=dst_crs,
        resampling=Resampling.nearest,
    )
    assert (out > 0).sum() == 383077
    assert_almost_equal(tuple(dst_transform),
                        tuple(Affine(330.2992903555146, 0.0, -8789636.707871985,
                                     0.0, -330.2992903555146, 2943560.2346221623)),
                        decimal=5)


def test_reproject_epsg__simple_array_resolution():
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        source = src.read(1)

    dst_crs = {"init": "EPSG:3857"}
    out, dst_transform = reproject(
        source,
        src_transform=src.transform,
        src_crs=src.crs,
        dst_crs=dst_crs,
        dst_resolution=(300, 300),
        resampling=Resampling.nearest,
    )
    assert (out > 0).sum() == 464503
    assert_almost_equal(tuple(dst_transform),
                        tuple(Affine(300, 0.0, -8789636.707871985,
                                     0.0, -300, 2943560.2346221623)),
                        decimal=5)


def test_reproject_epsg__simple_array_dst():
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        source = src.read(1)

    dst_crs = {"init": "EPSG:3857"}
    dst_out = np.empty(src.shape, dtype=np.uint8)

    out, dst_transform = reproject(
        source,
        dst_out,
        src_transform=src.transform,
        src_crs=src.crs,
        dst_crs=dst_crs,
        resampling=Resampling.nearest,
    )
    assert (out > 0).sum() == 368123
    assert_almost_equal(tuple(dst_transform),
                        tuple(Affine(335.3101519032594, 0.0, -8789636.707871985,
                                     0.0, -338.579773957742, 2943560.2346221623)),
                        decimal=5)


def test_reproject_epsg__simple():
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        dst_crs = {"init": "EPSG:3857"}
        out, dst_transform = reproject(
            rasterio.band(src, 1),
            dst_crs=dst_crs,
            resampling=Resampling.nearest,
        )
    assert (out > 0).sum() == 383077
    assert_almost_equal(tuple(dst_transform),
                        tuple(Affine(330.2992903555146, 0.0, -8789636.707871985,
                                     0.0, -330.2992903555146, 2943560.2346221623)),
                        decimal=5)


def test_reproject_epsg__simple_resolution():
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        dst_crs = {"init": "EPSG:3857"}
        out, dst_transform = reproject(
            rasterio.band(src, 1),
            dst_crs=dst_crs,
            dst_resolution=(300, 300),
            resampling=Resampling.nearest,
        )
    assert (out > 0).sum() == 464503
    assert_almost_equal(tuple(dst_transform),
                        tuple(Affine(300.0, 0.0, -8789636.707871985,
                                     0.0, -300.0, 2943560.2346221623)),
                        decimal=5)


def test_reproject_no_destination_with_transform():
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        source = src.read(1)

    dst_crs = {"init": "EPSG:3857"}
    with pytest.raises(ValueError):
        reproject(
            source,
            src_transform=src.transform,
            src_crs=src.crs,
            dst_crs=dst_crs,
            dst_transform=DST_TRANSFORM,
            resampling=Resampling.nearest,
        )


def test_reproject_out_of_bounds():
    """Using EPSG code is not appropriate for the transform.

    Should return blank image.
    """
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        source = src.read(1)

    dst_crs = {"init": "epsg:32619"}
    out = np.zeros(src.shape, dtype=np.uint8)
    reproject(
        source,
        out,
        src_transform=src.transform,
        src_crs=src.crs,
        dst_transform=DST_TRANSFORM,
        dst_crs=dst_crs,
        resampling=Resampling.nearest,
    )
    assert not out.any()


@requires_gdal3
@pytest.mark.parametrize("options, expected", reproj_expected)
def test_reproject_nodata(options, expected):
    # Older combinations of GDAL and PROJ might have got this transformation wrong.
    # Results look better with GDAL 3.
    nodata = 215

    with rasterio.Env(**options):
        params = uninvertable_reproject_params()
        source = np.ones((params.width, params.height), dtype=np.uint8)
        out = np.zeros((params.dst_width, params.dst_height), dtype=source.dtype)
        out.fill(120)  # Fill with arbitrary value

        reproject(
            source,
            out,
            src_transform=params.src_transform,
            src_crs=params.src_crs,
            src_nodata=nodata,
            dst_transform=params.dst_transform,
            dst_crs=params.dst_crs,
            dst_nodata=nodata,
        )

        assert (out == 1).sum() == expected
        assert (out == nodata).sum() == (
            params.dst_width * params.dst_height - expected
        )


@requires_gdal3
@pytest.mark.parametrize("options, expected", reproj_expected)
def test_reproject_nodata_nan(options, expected):

    with rasterio.Env(**options):
        params = uninvertable_reproject_params()
        source = np.ones((params.width, params.height), dtype=np.float32)
        out = np.zeros((params.dst_width, params.dst_height), dtype=source.dtype)
        out.fill(120)  # Fill with arbitrary value

        reproject(
            source,
            out,
            src_transform=params.src_transform,
            src_crs=params.src_crs,
            src_nodata=np.nan,
            dst_transform=params.dst_transform,
            dst_crs=params.dst_crs,
            dst_nodata=np.nan,
        )

        assert (out == 1).sum() == expected
        assert np.isnan(out).sum() == (params.dst_width * params.dst_height - expected)


@requires_gdal3
@pytest.mark.parametrize("options, expected", reproj_expected)
def test_reproject_dst_nodata_default(options, expected):
    """If nodata is not provided, destination will be filled with 0."""

    with rasterio.Env(**options):
        params = uninvertable_reproject_params()
        source = np.ones((params.width, params.height), dtype=np.uint8)
        out = np.zeros((params.dst_width, params.dst_height), dtype=source.dtype)
        out.fill(120)  # Fill with arbitrary value

        reproject(
            source,
            out,
            src_transform=params.src_transform,
            src_crs=params.src_crs,
            dst_transform=params.dst_transform,
            dst_crs=params.dst_crs,
        )

        assert (out == 1).sum() == expected
        assert (out == 0).sum() == (params.dst_width * params.dst_height - expected)


def test_reproject_invalid_dst_nodata():
    """dst_nodata must be in value range of data type."""
    params = default_reproject_params()

    source = np.ones((params.width, params.height), dtype=np.uint8)
    out = source.copy()

    with pytest.raises(ValueError):
        reproject(
            source,
            out,
            src_transform=params.src_transform,
            src_crs=params.src_crs,
            src_nodata=0,
            dst_transform=params.dst_transform,
            dst_crs=params.dst_crs,
            dst_nodata=999999999,
        )


def test_reproject_invalid_src_nodata():
    """src_nodata must be in range for data type."""
    params = default_reproject_params()

    source = np.ones((params.width, params.height), dtype=np.uint8)
    out = source.copy()

    with pytest.raises(ValueError):
        reproject(
            source,
            out,
            src_transform=params.src_transform,
            src_crs=params.src_crs,
            src_nodata=999999999,
            dst_transform=params.dst_transform,
            dst_crs=params.dst_crs,
            dst_nodata=215,
        )


def test_reproject_init_nodata_tofile(tmpdir):
    """Test that nodata is being initialized."""
    params = default_reproject_params()

    tiffname = str(tmpdir.join("foo.tif"))

    source1 = np.zeros((params.width, params.height), dtype=np.uint8)
    source2 = source1.copy()

    # fill both sources w/ arbitrary values
    rows, cols = source1.shape
    source1[:rows // 2, :cols // 2] = 200
    source2[rows // 2:, cols // 2:] = 100

    kwargs = {
        "count": 1,
        "width": params.width,
        "height": params.height,
        "dtype": np.uint8,
        "driver": "GTiff",
        "crs": params.dst_crs,
        "transform": params.dst_transform,
    }

    with rasterio.open(tiffname, "w", **kwargs) as dst:
        reproject(
            source1,
            rasterio.band(dst, 1),
            src_transform=params.src_transform,
            src_crs=params.src_crs,
            src_nodata=0.0,
            dst_transform=params.dst_transform,
            dst_crs=params.dst_crs,
            dst_nodata=0.0,
        )

        # 200s should be overwritten by 100s
        reproject(
            source2,
            rasterio.band(dst, 1),
            src_transform=params.src_transform,
            src_crs=params.src_crs,
            src_nodata=0.0,
            dst_transform=params.dst_transform,
            dst_crs=params.dst_crs,
            dst_nodata=0.0,
        )

    with rasterio.open(tiffname) as src:
        assert src.read().max() == 100


def test_reproject_no_init_nodata_tofile(tmpdir):
    """Test that nodata is not being initialized."""
    params = default_reproject_params()

    tiffname = str(tmpdir.join("foo.tif"))

    source1 = np.zeros((params.width, params.height), dtype=np.uint8)
    source2 = source1.copy()

    # fill both sources w/ arbitrary values
    rows, cols = source1.shape
    source1[:rows // 2, :cols // 2] = 200
    source2[rows // 2:, cols // 2:] = 100

    kwargs = {
        "count": 1,
        "width": params.width,
        "height": params.height,
        "dtype": np.uint8,
        "driver": "GTiff",
        "crs": params.dst_crs,
        "transform": params.dst_transform,
    }

    with rasterio.open(tiffname, "w", **kwargs) as dst:
        reproject(
            source1,
            rasterio.band(dst, 1),
            src_transform=params.src_transform,
            src_crs=params.src_crs,
            src_nodata=0.0,
            dst_transform=params.dst_transform,
            dst_crs=params.dst_crs,
            dst_nodata=0.0,
        )

        reproject(
            source2,
            rasterio.band(dst, 1),
            src_transform=params.src_transform,
            src_crs=params.src_crs,
            src_nodata=0.0,
            dst_transform=params.dst_transform,
            dst_crs=params.dst_crs,
            dst_nodata=0.0,
            init_dest_nodata=False,
        )

    # 200s should remain along with 100s
    with rasterio.open(tiffname) as src:
        data = src.read()

    assert data.max() == 200


def test_reproject_no_init_nodata_toarray():
    """Test that nodata is being initialized."""
    params = default_reproject_params()

    source1 = np.zeros((params.width, params.height))
    source2 = source1.copy()
    out = source1.copy()

    # fill both sources w/ arbitrary values
    rows, cols = source1.shape
    source1[:rows // 2, :cols // 2] = 200
    source2[rows // 2:, cols // 2:] = 100

    reproject(
        source1,
        out,
        src_transform=params.src_transform,
        src_crs=params.src_crs,
        src_nodata=0.0,
        dst_transform=params.dst_transform,
        dst_crs=params.dst_crs,
        dst_nodata=0.0,
    )

    assert out.max() == 200
    assert out.min() == 0

    reproject(
        source2,
        out,
        src_transform=params.src_transform,
        src_crs=params.src_crs,
        src_nodata=0.0,
        dst_transform=params.dst_transform,
        dst_crs=params.dst_crs,
        dst_nodata=0.0,
        init_dest_nodata=False,
    )

    # 200s should NOT be overwritten by 100s
    assert out.max() == 200
    assert out.min() == 0


def test_reproject_multi():
    """Ndarry to ndarray."""
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        source = src.read()
    dst_crs = dict(
        proj="merc",
        a=6378137,
        b=6378137,
        lat_ts=0.0,
        lon_0=0.0,
        x_0=0.0,
        y_0=0,
        k=1.0,
        units="m",
        nadgrids="@null",
        wktext=True,
        no_defs=True,
    )
    destin = np.empty(source.shape, dtype=np.uint8)
    reproject(
        source,
        destin,
        src_transform=src.transform,
        src_crs=src.crs,
        dst_transform=DST_TRANSFORM,
        dst_crs=dst_crs,
        resampling=Resampling.nearest,
    )
    assert destin.any()


def test_warp_from_file():
    """File to ndarray."""
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        dst_crs = dict(
            proj="merc",
            a=6378137,
            b=6378137,
            lat_ts=0.0,
            lon_0=0.0,
            x_0=0.0,
            y_0=0,
            k=1.0,
            units="m",
            nadgrids="@null",
            wktext=True,
            no_defs=True,
        )
        destin = np.empty(src.shape, dtype=np.uint8)
        reproject(
            rasterio.band(src, 1), destin, dst_transform=DST_TRANSFORM, dst_crs=dst_crs
        )
    assert destin.any()


def test_warp_from_to_file(tmpdir):
    """File to file."""
    tiffname = str(tmpdir.join("foo.tif"))
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        dst_crs = dict(
            proj="merc",
            a=6378137,
            b=6378137,
            lat_ts=0.0,
            lon_0=0.0,
            x_0=0.0,
            y_0=0,
            k=1.0,
            units="m",
            nadgrids="@null",
            wktext=True,
            no_defs=True,
        )
        kwargs = src.meta.copy()
        kwargs.update(transform=DST_TRANSFORM, crs=dst_crs)
        with rasterio.open(tiffname, "w", **kwargs) as dst:
            for i in (1, 2, 3):
                reproject(rasterio.band(src, i), rasterio.band(dst, i))


def test_warp_from_to_file_multi(tmpdir):
    """File to file."""
    tiffname = str(tmpdir.join("foo.tif"))
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        dst_crs = dict(
            proj="merc",
            a=6378137,
            b=6378137,
            lat_ts=0.0,
            lon_0=0.0,
            x_0=0.0,
            y_0=0,
            k=1.0,
            units="m",
            nadgrids="@null",
            wktext=True,
            no_defs=True,
        )
        kwargs = src.meta.copy()
        kwargs.update(transform=DST_TRANSFORM, crs=dst_crs)
        with rasterio.open(tiffname, "w", **kwargs) as dst:
            for i in (1, 2, 3):
                reproject(rasterio.band(src, i), rasterio.band(dst, i), num_threads=2)


@pytest.fixture(scope="function")
def polygon_3373():
    """An EPSG:3373 polygon."""
    return {
        "type": "Polygon",
        "coordinates": (
            (
                (798842.3090855901, 6569056.500655151),
                (756688.2826828464, 6412397.888771972),
                (755571.0617232556, 6408461.009397383),
                (677605.2284582685, 6425600.39266733),
                (677605.2284582683, 6425600.392667332),
                (670873.3791649605, 6427248.603432341),
                (664882.1106069803, 6407585.48425362),
                (663675.8662823177, 6403676.990080649),
                (485120.71963574126, 6449787.167760638),
                (485065.55660851026, 6449802.826920689),
                (485957.03982722526, 6452708.625101285),
                (487541.24541826674, 6457883.292107048),
                (531008.5797472061, 6605816.560367976),
                (530943.7197027118, 6605834.9333479265),
                (531888.5010308184, 6608940.750411527),
                (533299.5981959199, 6613962.642851984),
                (533403.6388841148, 6613933.172096095),
                (576345.6064638699, 6761983.708069147),
                (577649.6721159086, 6766698.137844516),
                (578600.3589008929, 6770143.99782289),
                (578679.4732294685, 6770121.638265098),
                (655836.640492081, 6749376.357102599),
                (659913.0791150068, 6764770.1314677475),
                (661105.8478791204, 6769515.168134831),
                (661929.4670843681, 6772800.8565198565),
                (661929.4670843673, 6772800.856519875),
                (661975.1582566603, 6772983.354777632),
                (662054.7979028501, 6772962.86384242),
                (841909.6014891531, 6731793.200435557),
                (840726.455490463, 6727039.8672589315),
                (798842.3090855901, 6569056.500655151),
            ),
        ),
    }


def test_transform_geom_polygon_cutting(polygon_3373):
    geom = polygon_3373
    result = transform_geom("EPSG:3373", "EPSG:4326", geom, antimeridian_cutting=True)
    assert result["type"] == "MultiPolygon"
    assert len(result["coordinates"]) == 2


def test_transform_geom_polygon_offset(polygon_3373):
    geom = polygon_3373
    result = transform_geom(
        "EPSG:3373", "EPSG:4326", geom, antimeridian_cutting=True, antimeridian_offset=0
    )
    assert result["type"] == "MultiPolygon"
    assert len(result["coordinates"]) == 2


def test_transform_geom_polygon_precision(polygon_3373):
    geom = polygon_3373
    result = transform_geom(
        "EPSG:3373", "EPSG:4326", geom, precision=1, antimeridian_cutting=True
    )
    assert all(round(x, 1) == x for x in flatten_coords(result["coordinates"]))


def test_transform_geom_linestring_precision(polygon_3373):
    ring = polygon_3373["coordinates"][0]
    geom = {"type": "LineString", "coordinates": ring}
    result = transform_geom(
        "EPSG:3373", "EPSG:4326", geom, precision=1, antimeridian_cutting=True
    )
    assert all(round(x, 1) == x for x in flatten_coords(result["coordinates"]))


def test_transform_geom_linestring_precision_iso(polygon_3373):
    ring = polygon_3373["coordinates"][0]
    geom = {"type": "LineString", "coordinates": ring}
    result = transform_geom("EPSG:3373", "EPSG:3373", geom, precision=1)
    assert int(result["coordinates"][0][0] * 10) == 7988423


def test_transform_geom_linearring_precision(polygon_3373):
    ring = polygon_3373["coordinates"][0]
    geom = {"type": "LinearRing", "coordinates": ring}
    result = transform_geom(
        "EPSG:3373", "EPSG:4326", geom, precision=1, antimeridian_cutting=True
    )
    assert all(round(x, 1) == x for x in flatten_coords(result["coordinates"]))


def test_transform_geom_linestring_precision_z(polygon_3373):
    ring = polygon_3373["coordinates"][0]
    x, y = zip(*ring)
    ring = list(zip(x, y, [0.0 for i in range(len(x))]))
    geom = {"type": "LineString", "coordinates": ring}
    result = transform_geom("EPSG:3373", "EPSG:3373", geom, precision=1)
    assert int(result["coordinates"][0][0] * 10) == 7988423
    assert int(result["coordinates"][0][2] * 10) == 0


def test_transform_geom_multipolygon(polygon_3373):
    geom = {"type": "MultiPolygon", "coordinates": [polygon_3373["coordinates"]]}
    result = transform_geom("EPSG:3373", "EPSG:4326", geom, precision=1)
    assert all(round(x, 1) == x for x in flatten_coords(result["coordinates"]))


def test_transform_geom_array(polygon_3373):
    geom = [polygon_3373 for _ in range(10)]
    result = transform_geom("EPSG:3373", "EPSG:4326", geom, precision=1)
    assert isinstance(result, list)
    assert len(result) == 10


def test_transform_geom__geo_interface(polygon_3373):
    class GeoObj:
        @property
        def __geo_interface__(self):
            return polygon_3373

    result = transform_geom("EPSG:3373", "EPSG:4326", GeoObj(), precision=1)
    assert all(round(x, 1) == x for x in flatten_coords(result["coordinates"]))


def test_transform_geom__geo_interface__array(polygon_3373):
    class GeoObj:
        @property
        def __geo_interface__(self):
            return polygon_3373

    geom = [GeoObj() for _ in range(10)]
    results = transform_geom("EPSG:3373", "EPSG:4326", geom, precision=1)
    assert isinstance(results, list)
    assert len(results) == 10
    for result in results:
        assert all(round(x, 1) == x for x in flatten_coords(result["coordinates"]))


@pytest.mark.parametrize("method", SUPPORTED_RESAMPLING)
def test_reproject_resampling(path_rgb_byte_tif, method):
    # Expected count of nonzero pixels for each resampling method, based
    # on running rasterio with each of the following configurations
    expected = {
        Resampling.nearest: [438113],
        Resampling.bilinear: [439280],
        Resampling.cubic: [437888],
        Resampling.cubic_spline: [440475],
        Resampling.lanczos: [436001],
        Resampling.average: [439419, 439172],  # latter value for GDAL 3.1
        Resampling.mode: [437298],
        Resampling.max: [439464],
        Resampling.min: [436397],
        Resampling.med: [437194],
        Resampling.q1: [436397],
        Resampling.q3: [438948],
        Resampling.sum: [439118],
        Resampling.rms: [439385],
    }

    with rasterio.open(path_rgb_byte_tif) as src:
        source = src.read(1)

    out = np.empty(src.shape, dtype=np.uint8)
    reproject(
        source,
        out,
        src_transform=src.transform,
        src_crs=src.crs,
        dst_transform=DST_TRANSFORM,
        dst_crs="EPSG:3857",
        resampling=method,
    )

    assert np.count_nonzero(out) in expected[method]


@pytest.mark.parametrize("test3d,count_nonzero", [(True, 1309625), (False, 437686)])
def test_reproject_array_interface(test3d, count_nonzero, path_rgb_byte_tif):
    class DataArray:
        def __init__(self, data):
            self.data = data

        def __array__(self, dtype=None):
            return self.data

        @property
        def dtype(self):
            return self.data.dtype

    with rasterio.open(path_rgb_byte_tif) as src:
        if test3d:
            source = DataArray(src.read())
        else:
            source = DataArray(src.read(1))
    out = DataArray(np.empty(source.data.shape, dtype=np.uint8))
    reproject(
        source,
        out,
        src_transform=src.transform,
        src_crs=src.crs,
        src_nodata=src.nodata,
        dst_transform=DST_TRANSFORM,
        dst_crs="EPSG:3857",
        dst_nodata=99,
    )
    assert isinstance(out, DataArray)
    assert np.count_nonzero(out.data[out.data != 99]) == count_nonzero


@pytest.mark.parametrize("test3d,count_nonzero", [(True, 1309625), (False, 437686)])
def test_reproject_masked(test3d, count_nonzero, path_rgb_byte_tif):
    with rasterio.open(path_rgb_byte_tif) as src:
        if test3d:
            source = src.read(masked=True)
        else:
            source = src.read(1, masked=True)
    out = np.empty(source.shape, dtype=np.uint8)
    reproject(
        source,
        out,
        src_transform=src.transform,
        src_crs=src.crs,
        dst_transform=DST_TRANSFORM,
        dst_crs="EPSG:3857",
        dst_nodata=99,
    )
    assert np.ma.is_masked(source)
    assert np.count_nonzero(out[out != 99]) == count_nonzero


@pytest.mark.parametrize("method", SUPPORTED_RESAMPLING)
def test_reproject_resampling_alpha(method):
    """Reprojection of a source with alpha band succeeds"""
    # Expected count of nonzero pixels for each resampling method, based
    # on running rasterio with each of the following configurations
    expected = {
        Resampling.nearest: [438113],
        Resampling.bilinear: [439280],
        Resampling.cubic: [437888],
        Resampling.cubic_spline: [440475],
        Resampling.lanczos: [436001],
        Resampling.average: [439419, 439172],  # latter value for GDAL 3.1
        Resampling.mode: [437298],
        Resampling.max: [439464],
        Resampling.min: [436397],
        Resampling.med: [437194],
        Resampling.q1: [436397],
        Resampling.q3: [438948],
        Resampling.sum: [439118],
        Resampling.rms: [439385],
    }

    with rasterio.open("tests/data/RGBA.byte.tif") as src:
        source = src.read(1)

    out = np.empty(src.shape, dtype=np.uint8)
    reproject(
        source,
        out,
        src_transform=src.transform,
        src_crs=src.crs,
        dst_transform=DST_TRANSFORM,
        dst_crs="EPSG:3857",
        resampling=method,
    )

    assert np.count_nonzero(out) in expected[method]


@pytest.mark.skipif(
    gdal_version.at_least("2.0"), reason="Tests only applicable to GDAL < 2.0"
)
@pytest.mark.parametrize("method", GDAL2_RESAMPLING)
def test_reproject_not_yet_supported_resampling(method):
    """Test resampling methods not yet supported by this version of GDAL"""
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        source = src.read(1)

    dst_crs = "EPSG:32619"
    out = np.empty(src.shape, dtype=np.uint8)
    with pytest.raises(GDALVersionError):
        reproject(
            source,
            out,
            src_transform=src.transform,
            src_crs=src.crs,
            dst_transform=DST_TRANSFORM,
            dst_crs=dst_crs,
            resampling=method,
        )


def test_reproject_unsupported_resampling():
    """Values not in enums. Resampling are not supported."""
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        source = src.read(1)

    dst_crs = "EPSG:32619"
    out = np.empty(src.shape, dtype=np.uint8)
    with pytest.raises(ValueError):
        reproject(
            source,
            out,
            src_transform=src.transform,
            src_crs=src.crs,
            dst_transform=DST_TRANSFORM,
            dst_crs=dst_crs,
            resampling=99,
        )


def test_reproject_unsupported_resampling_guass():
    """Resampling.gauss is unsupported."""
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        source = src.read(1)

    dst_crs = "EPSG:32619"
    out = np.empty(src.shape, dtype=np.uint8)
    with pytest.raises(ValueError):
        reproject(
            source,
            out,
            src_transform=src.transform,
            src_crs=src.crs,
            dst_transform=DST_TRANSFORM,
            dst_crs=dst_crs,
            resampling=Resampling.gauss,
        )


@pytest.mark.parametrize("method", SUPPORTED_RESAMPLING)
def test_resample_default_invert_proj(method):
    """Nearest and bilinear should produce valid results
    with the default Env
    """

    with rasterio.open("tests/data/world.rgb.tif") as src:
        source = src.read(1)
        profile = src.profile

    dst_crs = "EPSG:32619"

    # Calculate the ideal dimensions and transformation in the new crs
    dst_affine, dst_width, dst_height = calculate_default_transform(
        src.crs, dst_crs, src.width, src.height, *src.bounds
    )

    profile["height"] = dst_height
    profile["width"] = dst_width

    out = np.empty(shape=(dst_height, dst_width), dtype=np.uint8)

    # GDAL 1.11 needs to have this config option set on to match the
    # default results in later versions.
    if gdal_version.major == 1:
        options = dict(CHECK_WITH_INVERT_PROJ=True)
    else:
        options = {}

    with rasterio.Env(**options):
        reproject(
            source,
            out,
            src_transform=src.transform,
            src_crs=src.crs,
            dst_transform=dst_affine,
            dst_crs=dst_crs,
            resampling=method,
        )

    assert out.mean() > 0


@pytest.mark.xfail(reason="Projection extents have changed with PROJ 8")
def test_target_aligned_pixels():
    """Issue 853 has been resolved"""
    with rasterio.open("tests/data/world.rgb.tif") as src:
        source = src.read(1)
        profile = src.profile

    dst_crs = "EPSG:3857"

    with rasterio.Env(CHECK_WITH_INVERT_PROJ=False):
        # Calculate the ideal dimensions and transformation in the new crs
        dst_affine, dst_width, dst_height = calculate_default_transform(
            src.crs, dst_crs, src.width, src.height, *src.bounds
        )

        dst_affine, dst_width, dst_height = aligned_target(
            dst_affine, dst_width, dst_height, 10000.0
        )

        profile["height"] = dst_height
        profile["width"] = dst_width

        out = np.empty(shape=(dst_height, dst_width), dtype=np.uint8)

        reproject(
            source,
            out,
            src_transform=src.transform,
            src_crs=src.crs,
            dst_transform=dst_affine,
            dst_crs=dst_crs,
            resampling=Resampling.nearest,
        )

        # Check that there are no black borders
        assert out[:, 0].all()
        assert out[:, -1].all()
        assert out[0, :].all()
        assert out[-1, :].all()


@pytest.mark.parametrize("method", SUPPORTED_RESAMPLING)
def test_resample_no_invert_proj(method):
    """Nearest and bilinear should produce valid results with
    CHECK_WITH_INVERT_PROJ = False
    """

    if method in (
        Resampling.bilinear,
        Resampling.cubic,
        Resampling.cubic_spline,
        Resampling.lanczos,
    ):
        pytest.xfail(
            reason="Some resampling methods succeed but produce blank images. "
            "See https://github.com/rasterio/rasterio/issues/614"
        )

    with rasterio.Env(CHECK_WITH_INVERT_PROJ=False):
        with rasterio.open("tests/data/world.rgb.tif") as src:
            source = src.read(1)
            profile = src.profile.copy()

        dst_crs = "EPSG:32619"

        # Calculate the ideal dimensions and transformation in the new crs
        dst_affine, dst_width, dst_height = calculate_default_transform(
            src.crs, dst_crs, src.width, src.height, *src.bounds
        )

        profile["height"] = dst_height
        profile["width"] = dst_width

        out = np.empty(shape=(dst_height, dst_width), dtype=np.uint8)

        # see #614, some resampling methods succeed but produce blank images
        out = np.empty(src.shape, dtype=np.uint8)
        reproject(
            source,
            out,
            src_transform=src.transform,
            src_crs=src.crs,
            dst_transform=dst_affine,
            dst_crs=dst_crs,
            resampling=method,
        )

        assert out.mean() > 0


def test_reproject_crs_none():
    """Reproject with crs is None should not cause segfault"""
    src = np.random.random(25).reshape((1, 5, 5))
    srcaff = Affine(1.1, 0.0, 0.0, 0.0, 1.1, 0.0)
    srccrs = None
    dst = np.empty(shape=(1, 11, 11))
    dstaff = Affine(0.5, 0.0, 0.0, 0.0, 0.5, 0.0)
    dstcrs = None

    with pytest.raises(ValueError):
        reproject(
            src,
            dst,
            src_transform=srcaff,
            src_crs=srccrs,
            dst_transform=dstaff,
            dst_crs=dstcrs,
            resampling=Resampling.nearest,
        )


def test_reproject_identity_src():
    """Reproject with an identity like source matrices."""
    src = np.random.random(25).reshape((1, 5, 5))
    dst = np.empty(shape=(1, 10, 10))
    dstaff = Affine(0.5, 0.0, 0.0, 0.0, 0.5, 0.0)
    crs = {"init": "epsg:3857"}

    src_affines = [
        Affine(1.0, 0.0, 0.0, 0.0, 1.0, 0.0),  # Identity both positive
        Affine(1.0, 0.0, 0.0, 0.0, -1.0, 0.0),  # Identity with negative e
    ]

    for srcaff in src_affines:
        # reproject expected to not raise any error in any of the srcaff
        reproject(
            src,
            dst,
            src_transform=srcaff,
            src_crs=crs,
            dst_transform=dstaff,
            dst_crs=crs,
            resampling=Resampling.nearest,
        )


def test_reproject_identity_dst():
    """Reproject with an identity like destination matrices."""
    src = np.random.random(100).reshape((1, 10, 10))
    srcaff = Affine(0.5, 0.0, 0.0, 0.0, 0.5, 0.0)
    dst = np.empty(shape=(1, 5, 5))
    crs = {"init": "epsg:3857"}

    dst_affines = [
        Affine(1.0, 0.0, 0.0, 0.0, 1.0, 0.0),  # Identity both positive
        Affine(1.0, 0.0, 0.0, 0.0, -1.0, 0.0),  # Identity with negative e
    ]

    for dstaff in dst_affines:
        # reproject expected to not raise any error in any of the dstaff
        reproject(
            src,
            dst,
            src_transform=srcaff,
            src_crs=crs,
            dst_transform=dstaff,
            dst_crs=crs,
            resampling=Resampling.nearest,
        )


@pytest.fixture(scope="function")
def rgb_byte_profile():
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        return src.profile


def test_reproject_gcps_transform_exclusivity():
    """gcps and transform can't be used together."""
    with pytest.raises(ValueError):
        reproject(1, 1, gcps=[0], src_transform=[0])


def test_reproject_gcps(rgb_byte_profile):
    """Reproject using ground control points for the source"""
    source = np.ones((3, 800, 800), dtype=np.uint8) * 255
    out = np.zeros(
        (3, rgb_byte_profile["height"], rgb_byte_profile["height"]), dtype=np.uint8
    )
    src_gcps = [
        GroundControlPoint(row=0, col=0, x=156113, y=2818720, z=0),
        GroundControlPoint(row=0, col=800, x=338353, y=2785790, z=0),
        GroundControlPoint(row=800, col=800, x=297939, y=2618518, z=0),
        GroundControlPoint(row=800, col=0, x=115698, y=2651448, z=0),
    ]
    reproject(
        source,
        out,
        src_crs="EPSG:32618",
        gcps=src_gcps,
        dst_transform=rgb_byte_profile["transform"],
        dst_crs=rgb_byte_profile["crs"],
        resampling=Resampling.nearest,
    )

    assert not out.all()
    assert not out[:, 0, 0].any()
    assert not out[:, 0, -1].any()
    assert not out[:, -1, -1].any()
    assert not out[:, -1, 0].any()


@requires_gdal22(
    reason="GDAL 2.2.0 and newer has different antimeridian cutting behavior."
)
def test_transform_geom_gdal22():
    """Enabling `antimeridian_cutting` has no effect on GDAL 2.2.0 or newer
    where antimeridian cutting is always enabled.  This could produce
    unexpected geometries, so an exception is raised.
    """
    geom = {"type": "Point", "coordinates": [0, 0]}
    with pytest.raises(GDALVersionError):
        transform_geom("EPSG:4326", "EPSG:3857", geom, antimeridian_cutting=False)


def test_issue1056():
    """Warp sucessfully from RGB's upper bands to an array"""
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        dst_crs = "EPSG:3857"
        out = np.zeros(src.shape, dtype=np.uint8)
        reproject(
            rasterio.band(src, 2),
            out,
            src_transform=src.transform,
            src_crs=src.crs,
            dst_transform=DST_TRANSFORM,
            dst_crs=dst_crs,
            resampling=Resampling.nearest,
        )


def test_reproject_dst_nodata():
    """Affirm resolution of issue #1395"""
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        source = src.read(1)

    dst_crs = "EPSG:3857"
    out = np.empty(src.shape, dtype=np.float32)
    reproject(
        source,
        out,
        src_transform=src.transform,
        src_crs=src.crs,
        dst_transform=DST_TRANSFORM,
        dst_crs=dst_crs,
        src_nodata=0,
        dst_nodata=np.nan,
        resampling=Resampling.nearest,
    )

    assert (out[~np.isnan(out)] > 0.0).sum() == 438113
    assert out[0, 0] != 0
    assert np.isnan(out[0, 0])


def test_issue1401():
    """The warp_mem_limit keyword argument is in effect"""
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        dst_crs = "EPSG:3857"
        out = np.zeros(src.shape, dtype=np.uint8)
        reproject(
            rasterio.band(src, 2),
            out,
            src_transform=src.transform,
            src_crs=src.crs,
            dst_transform=DST_TRANSFORM,
            dst_crs=dst_crs,
            resampling=Resampling.nearest,
            warp_mem_limit=4000,
        )


def test_reproject_dst_alpha(path_rgb_msk_byte_tif):
    """Materialization of external mask succeeds"""

    with rasterio.open(path_rgb_msk_byte_tif) as src:

        nrows, ncols = src.shape

        dst_arr = np.zeros((src.count + 1, nrows, ncols), dtype=np.uint8)

        reproject(
            rasterio.band(src, src.indexes),
            dst_arr,
            src_transform=src.transform,
            src_crs=src.crs,
            dst_transform=DST_TRANSFORM,
            dst_crs="EPSG:3857",
            dst_alpha=4,
        )

        assert dst_arr[3].any()


@pytest.mark.xfail(
    rasterio.__gdal_version__ in ["2.2.0", "2.2.1", "2.2.2", "2.2.3"],
    reason=(
        "GDAL had regression in 2.2.X series, fixed in 2.2.4,"
        " reproject used dst index instead of src index when destination was single band"
    ),
)
def test_issue1350():
    """Warp bands other than 1 or All"""

    with rasterio.open("tests/data/RGB.byte.tif") as src:
        dst_crs = "EPSG:3857"

        reprojected = []

        for dtype, idx in zip(src.dtypes, src.indexes):
            out = np.zeros((1,) + src.shape, dtype=dtype)

            reproject(
                rasterio.band(src, idx),
                out,
                resampling=Resampling.nearest,
                dst_transform=DST_TRANSFORM,
                dst_crs=dst_crs,
            )

            reprojected.append(out)

        for i in range(1, len(reprojected)):
            assert not (reprojected[0] == reprojected[i]).all()


def test_issue_1446():
    """Confirm resolution of #1446"""
    g = transform_geom(
        CRS.from_epsg(4326),
        CRS.from_epsg(32610),
        {"type": "Point", "coordinates": (-122.51403808499907, 38.06106733107932)},
    )
    assert round(g["coordinates"][0], 1) == 542630.9
    assert round(g["coordinates"][1], 1) == 4212702.1


@requires_gdal_lt_3
def test_issue_1446_b():
    """Confirm that lines aren't thrown as reported in #1446"""
    src_crs = CRS.from_epsg(4326)
    dst_crs = CRS(
        {
            "proj": "sinu",
            "lon_0": 350.85607029556,
            "x_0": 0,
            "y_0": 0,
            "a": 3396190,
            "b": 3396190,
            "units": "m",
            "no_defs": True,
        }
    )
    collection = json.load(open("tests/data/issue1446.geojson"))
    geoms = {f["properties"]["fid"]: f["geometry"] for f in collection["features"]}
    transformed_geoms = {
        k: transform_geom(src_crs, dst_crs, g) for k, g in geoms.items()
    }
    # Before the fix, this geometry was thrown eastward of 0.0. It should be between -350 and -250.
    assert all([-350 < x < -150 for x, y in transformed_geoms[183519]["coordinates"]])


def test_reproject_init_dest_nodata():
    """No pixels should transfer over"""
    crs = CRS.from_epsg(4326)
    transform = Affine.identity()
    source = np.zeros((1, 100, 100))
    destination = np.ones((1, 100, 100))
    reproject(
        source, destination, src_crs=crs, src_transform=transform,
        dst_crs=crs, dst_transform=transform,
        src_nodata=0, init_dest_nodata=False
    )
    assert destination.all()


def test_empty_transform_inputs():
    """Check for fix of #1952"""
    assert ([], []) == rasterio.warp.transform(
        "EPSG:3857", "EPSG:4326", [], [], zs=None
    )


def test_empty_transform_inputs_z():
    """Check for fix of #1952"""
    assert ([], [], []) == rasterio.warp.transform(
        "EPSG:3857", "EPSG:4326", [], [], zs=[]
    )


def test_empty_transform_inputs_length():
    """Get an exception of inputs have different lengths"""
    with pytest.raises(TransformError):
        rasterio.warp.transform("EPSG:3857", "EPSG:4326", [1], [1, 2])


def test_empty_transform_inputs_length_z():
    """Get an exception of inputs have different lengths"""
    with pytest.raises(TransformError):
        rasterio.warp.transform("EPSG:3857", "EPSG:4326", [1, 2], [1, 2], zs=[0])


def test_reproject_rpcs(caplog):
    """Reproject using rational polynomial coefficients for the source"""
    with rasterio.open('tests/data/RGB.byte.rpc.vrt') as src:
        out = np.zeros(
            (3, src.profile["width"], src.profile["height"]), dtype=np.uint8
        )
        src_rpcs = src.rpcs
        reproject(
            rasterio.band(src, src.indexes),
            out,
            src_crs="EPSG:4326",
            rpcs=src_rpcs,
            dst_crs="EPSG:3857",
            resampling=Resampling.nearest,
        )

        assert not out.all()
        assert not out[:, 0, 0].any()
        assert not out[:, 0, -1].any()
        assert not out[:, -1, -1].any()
        assert not out[:, -1, 0].any()


def test_reproject_rpcs_with_transformer_options(caplog):
    """Reproject using rational polynomial coefficients and additional transformer options"""
    with rasterio.open('tests/data/RGB.byte.rpc.vrt') as src:
        with rasterio.MemoryFile(dirname='foo', filename='dem.tif') as mem:
            crs = 'COMPD_CS["WGS 84 + EGM96 height",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]],VERT_CS["EGM96 height",VERT_DATUM["EGM96 geoid",2005,AUTHORITY["EPSG","5171"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Gravity-related height",UP],AUTHORITY["EPSG","5773"]]]'
            transform = Affine(0.001953396267361111, 0.0, -124.00013888888888, 0.0, -0.001953396267361111, 50.000138888888884)
            with mem.open(
                driver="GTiff",
                width=1024,
                height=1024,
                count=1,
                transform=transform,
                dtype="int16",
                crs=crs,
            ) as dem:
                # we flush dem dataset before letting GDAL read from vsimem
                dem.write_band(1, 500 * np.ones((1024, 1024), dtype='int16'))

            out = np.zeros(
                (3, src.profile["width"], src.profile["height"]), dtype=np.uint8
            )
            out2 = out.copy()
            src_rpcs = src.rpcs
            caplog.set_level(logging.DEBUG)
            reproject(
                rasterio.band(src, src.indexes),
                out,
                src_crs="EPSG:4326",
                rpcs=src_rpcs,
                dst_crs="EPSG:3857",
                resampling=Resampling.nearest,
                RPC_DEM=dem.name,

            )
            caplog.set_level(logging.INFO)
            reproject(
                rasterio.band(src, src.indexes),
                out2,
                src_crs="EPSG:4326",
                rpcs=src_rpcs,
                dst_crs="EPSG:3857",
                resampling=Resampling.nearest,

            )

            assert not out.all()
            assert not out2.all()
            assert "RPC_DEM" in caplog.text
            assert not np.array_equal(out, out2)


def test_warp_gcps_compute_dst_transform_automatically_array():
    """Ensure we don't raise an exception when gcps passed without dst_transform, for a source array"""
    source = np.ones((3, 800, 800), dtype=np.uint8) * 255
    out = np.zeros((3, 512, 512))
    src_gcps = [
        GroundControlPoint(row=0, col=0, x=156113, y=2818720, z=0),
        GroundControlPoint(row=0, col=800, x=338353, y=2785790, z=0),
        GroundControlPoint(row=800, col=800, x=297939, y=2618518, z=0),
        GroundControlPoint(row=800, col=0, x=115698, y=2651448, z=0),
    ]
    reproject(
        source,
        out,
        src_crs="EPSG:32618",
        gcps=src_gcps,
        dst_crs="EPSG:32618",
        resampling=Resampling.nearest
    )

    assert not out.all()
    assert not out[:, 0, 0].any()
    assert not out[:, 0, -1].any()
    assert not out[:, -1, -1].any()
    assert not out[:, -1, 0].any()


def test_warp_gcps_compute_dst_transform_automatically_reader(tmpdir):
    """Ensure we don't raise an exception when gcps passed without dst_transform, for a source dataset"""
    tiffname = str(tmpdir.join('test.tif'))
    src_gcps = [
        GroundControlPoint(row=0, col=0, x=156113, y=2818720, z=0),
        GroundControlPoint(row=0, col=800, x=338353, y=2785790, z=0),
        GroundControlPoint(row=800, col=800, x=297939, y=2618518, z=0),
        GroundControlPoint(row=800, col=0, x=115698, y=2651448, z=0),
    ]
    out = np.zeros((3, 512, 512))

    with rasterio.open(tiffname, mode='w', height=800, width=800, count=3, dtype=np.uint8) as source:
        source.gcps = (src_gcps, CRS.from_epsg(32618))

    with rasterio.open(tiffname) as source:
        reproject(
            rasterio.band(source, source.indexes),
            out,
            dst_crs="EPSG:32618",
            resampling=Resampling.nearest
        )

    assert not out.all()
    assert not out[:, 0, 0].any()
    assert not out[:, 0, -1].any()
    assert not out[:, -1, -1].any()
    assert not out[:, -1, 0].any()


def test_reproject_rpcs_exact_transformer(caplog):
    """Reproject using rational polynomial coefficients and DEM, requiring that
       we don't try to make an approximate transformer.
    """
    with rasterio.open('tests/data/RGB.byte.rpc.vrt') as src:
        with rasterio.MemoryFile(dirname='foo', filename='dem.tif') as mem:
            crs = 'COMPD_CS["WGS 84 + EGM96 height",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]],VERT_CS["EGM96 height",VERT_DATUM["EGM96 geoid",2005,AUTHORITY["EPSG","5171"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Gravity-related height",UP],AUTHORITY["EPSG","5773"]]]'
            transform = Affine(0.001953396267361111, 0.0, -124.00013888888888, 0.0, -0.001953396267361111, 50.000138888888884)
            with mem.open(
                driver="GTiff",
                width=1024,
                height=1024,
                count=1,
                transform=transform,
                dtype="int16",
                crs=crs,
            ) as dem:
                # we flush dem dataset before letting GDAL read from vsimem
                dem.write_band(1, 500 * np.ones((1024, 1024), dtype='int16'))

            out = np.zeros(
                (3, src.profile["width"], src.profile["height"]), dtype=np.uint8
            )
            src_rpcs = src.rpcs
            caplog.set_level(logging.DEBUG)
            reproject(
                rasterio.band(src, src.indexes),
                out,
                src_crs="EPSG:4326",
                rpcs=src_rpcs,
                dst_crs="EPSG:3857",
                resampling=Resampling.nearest,
                RPC_DEM=dem.name,
            )

            assert "Created exact transformer" in caplog.text


def test_reproject_rpcs_approx_transformer(caplog):
    """Reproject using rational polynomial coefficients without a DEM, for which it's
       ok to use an approximate transformer.
    """
    with rasterio.open('tests/data/RGB.byte.rpc.vrt') as src:
        out = np.zeros(
                (3, src.profile["width"], src.profile["height"]), dtype=np.uint8
            )
        src_rpcs = src.rpcs
        caplog.set_level(logging.DEBUG)
        reproject(
            rasterio.band(src, src.indexes),
            out,
            src_crs="EPSG:4326",
            rpcs=src_rpcs,
            dst_crs="EPSG:3857",
            resampling=Resampling.nearest,
        )

        assert "Created approximate transformer" in caplog.text


@pytest.fixture
def http_error_server(data):
    """Serves files from the test data directory, poorly."""
    import functools
    import multiprocessing
    import http.server
    import os


    PORT = 8000
    Handler = functools.partial(RangeRequestErrorHandler, directory=str(data))
    httpd = http.server.HTTPServer(("", PORT), Handler)
    p = multiprocessing.Process(target=httpd.serve_forever)
    p.start()
    yield
    p.terminate()
    p.join()


@requires_gdal3
@pytest.mark.skipif(
    sys.version_info < (3, 7) and sys.platform != "linux",
    reason="Python 3.7 required to serve the data fixture directory and the server fixture requires Linux",
)
def test_reproject_error_propagation(http_error_server, caplog):
    """Propagate errors up from ChunkAndWarpMulti and check for a retry."""

    with rasterio.open(
        "/vsicurl?max_retry=1&retry_delay=.1&url=http://localhost:8000/RGB.byte.tif"
    ) as src:
        out = np.zeros((src.count, src.height, src.width), dtype="uint8")

        with pytest.raises(WarpOperationError):
            reproject(
                rasterio.band(src, (1, 2, 3)),
                out,
                dst_crs=src.crs,
                dst_transform=src.transform,
            )

    assert len([rec for rec in caplog.records if "Retrying again" in rec.message]) == 2