def test_dimensions_missing_params():
    """dst_width and dst_height must be specified together"""
    with pytest.raises(ValueError):
        calculate_default_transform(
            'epsg:4326', 'epsg:3857', width=1, height=1, gcps=[1],
            resolution=1, dst_width=1, dst_height=None)

    with pytest.raises(ValueError):
        calculate_default_transform(
            'epsg:4326', 'epsg:3857', width=1, height=1, gcps=[1],
            resolution=1, dst_width=None, dst_height=1)

def create_image_source(origin_uri, source_uri, image_folder, order, tile_dim):
    with rasterio.drivers():
        with rasterio.open(source_uri) as src:
            shape = src.shape
            res = src.res
            bounds = src.bounds
            (ll_transform, ll_cols, ll_rows) = calculate_default_transform(
                src.crs,
                "EPSG:4326",
                src.shape[0],
                src.shape[1],
                src.bounds.left,
                src.bounds.bottom,
                src.bounds.right,
                src.bounds.top,
            )
            w, n = ll_transform.xoff, ll_transform.yoff
            e, s = ll_transform * (ll_cols, ll_rows)
            ll_bounds = [w, s, e, n]

            (wm_transform, _, _) = calculate_default_transform(
                src.crs,
                "EPSG:3857",
                src.shape[0],
                src.shape[1],
                src.bounds.left,
                src.bounds.bottom,
                src.bounds.right,
                src.bounds.top,
            )

            resolution = max(abs(wm_transform[0]), abs(wm_transform[4]))
            zoom = get_zoom(resolution, tile_dim)
            min_tile = mercantile.tile(ll_bounds[0], ll_bounds[3], zoom)
            max_tile = mercantile.tile(ll_bounds[2], ll_bounds[1], zoom)

            return ImageSource(
                origin_uri=origin_uri,
                source_uri=source_uri,
                src_bounds=src.bounds,
                src_shape=src.shape,
                src_crs=src.crs,
                zoom=zoom,
                ll_bounds=ll_bounds,
                tile_bounds=[min_tile.x, min_tile.y, max_tile.x, max_tile.y],
                image_folder=image_folder,
                order=order,
            )

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.copy()

    dst_crs = {"init": "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)

    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_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,
            {"init": "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_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.copy()

    dst_crs = {'init': '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, 100000.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 is no black borders
        assert out[:, 0].all()
        assert out[:, -1].all()
        assert out[0, :].all()
        assert out[-1, :].all()

def test_resample_no_invert_proj(method):
    """Nearest and bilinear should produce valid results with
    CHECK_WITH_INVERT_PROJ = False
    """
    if not supported_resampling(method):
        pytest.skip()

    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 = {'init': '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 resamplin 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 warp_tif(combined_tif_path, warped_tif_path, dst_crs={
        'init': 'EPSG:3857'
}):
    logger.info('Warping tif to web mercator: %s', combined_tif_path)
    with rasterio.open(combined_tif_path) as src:
        meta = src.meta
        new_meta = meta.copy()
        transform, width, height = calculate_default_transform(
            src.crs, dst_crs, src.width, src.height, *src.bounds)
        new_meta.update({
            'crs': dst_crs,
            'transform': transform,
            'width': width,
            'height': height,
            'nodata': -28762
        })
        with rasterio.open(
                warped_tif_path, 'w', compress='DEFLATE', tiled=True,
                **new_meta) as dst:
            for i in range(1, src.count):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=transform,
                    dst_crs=dst_crs,
                    resampling=Resampling.nearest,
                    src_nodata=-28762
                )

def project_raster(src_raster, dst_raster, dst_crs,
                   resampling=1, resolution=None, num_threads=2):
    """Reproject a raster from one coordinate system to another using Rasterio
    code from: https://github.com/mapbox/rasterio/blob/master/docs/reproject.rst

    Parameters
    ----------
    src_raster : str
        Filename of source raster.
    dst_raster : str
        Filename of reprojected (destination) raster.
    dst_crs : str
        Coordinate system of reprojected raster.
        Examples:
            'EPSG:26715'
    resampling : int (see rasterio source code: https://github.com/mapbox/rasterio/blob/master/rasterio/enums.py)
        nearest = 0
        bilinear = 1
        cubic = 2
        cubic_spline = 3
        lanczos = 4
        average = 5
        mode = 6
        gauss = 7
        max = 8
        min = 9
        med = 10
        q1 = 11
        q3 = 12
    resolution : tuple of floats (len 2)
        cell size of the output raster
        (x resolution, y resolution)
    """
    rasterio = import_rasterio() # check for rasterio
    from rasterio.warp import calculate_default_transform, reproject

    with rasterio.open(src_raster) as src:
        affine, width, height = calculate_default_transform(
            src.crs, dst_crs, src.width, src.height, *src.bounds, resolution=resolution)
        kwargs = src.meta.copy()
        kwargs.update({
            'crs': dst_crs,
            'transform': affine,
            'affine': affine,
            'width': width,
            'height': height
        })
        with rasterio.open(dst_raster, 'w', **kwargs) as dst:
            for i in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.affine,
                    src_crs=src.crs,
                    dst_transform=affine,
                    dst_crs=dst_crs,
                    resampling=resampling,
                    num_threads=num_threads)

    def reproject(self, destination_file, source_file=None, resampling=RESAMPLING.nearest, **kwargs):
        """
        Reprojects the pixels of a source raster map to a destination raster, with a different reference coordinate
        system and Affine transform. It uses `Rasterio <https://github.com/mapbox/rasterio/blob/master/docs/reproject.rst>`_
        calculate_default_transform() to calculate parameters such as the resolution (if not provided), and the destination
        transform and dimensions.

        param string source_file: Full path to the source file containing a raster map

        param string destination_file: Full path to the destination file containing a raster map

        :param int resampling: Resampling method to use. Can be one of the following: ``Resampling.nearest``, ``Resampling.bilinear``, \
        ``Resampling.cubic``, ``Resampling.cubic_spline``, ``Resampling.lanczos``, ``Resampling.average``, ``Resampling.mode``.

        :param dict kwargs: Optional additional arguments passed to the method, to parametrize the reprojection. \
        For example: :attr:`dst_crs` for the target coordinate reference system, :attr:`resolution` for the target resolution, \
        in units of target coordinate reference system.

        """
        if not source_file:
            if not self.file_path:
                raise AttributeError("Please provide a source_file to load the data from.")
            else:
                source_file = self.file_path

        with rasterio.open(source_file) as src:
            affine, width, height = calculate_default_transform(src_crs=src.crs,
                                                                dst_crs=kwargs.get('dst_crs', src.crs),
                                                                width=kwargs.get('width', src.width),
                                                                height=kwargs.get('height', src.height),
                                                                left=kwargs.get('left', src.bounds.left),
                                                                bottom=kwargs.get('bottom', src.bounds.bottom),
                                                                right=kwargs.get('right', src.bounds.right),
                                                                top=kwargs.get('top', src.bounds.top),
                                                                resolution=kwargs.get('resolution', src.res)
                                                                )
            logger.info("Calculated default transformation:")
            logger.info("Affine:\n{0} \n width={1}, height={2}".format(affine, width, height))

            kwargs = src.meta.copy()
            kwargs.update({'transform': affine,
                           'affine': affine,
                           'width': width,
                           'height': height
                           })

            with rasterio.open(destination_file, 'w', **kwargs) as dst:
                for i in range(1, src.count + 1):
                    rasterio.warp.reproject(source=rasterio.band(src, i),
                                            destination=rasterio.band(dst, i),
                                            src_transform=src.affine,
                                            src_crs=src.crs,
                                            dst_transform=affine,
                                            dst_crs=kwargs.get('dst_crs', src.crs),
                                            resampling=resampling
                                            )
            logger.info("Reprojected data in %s " % destination_file)

def test_gcps_calculate_transform():
    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)]
    _, width, height = calculate_default_transform(
        'epsg:3857', 'epsg:4326', width=800, height=800, gcps=src_gcps)
    assert width == 1087
    assert height == 895

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 = {'init': '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 __compute_transform(self, in_raster, new_crs):
     affine, width, height = calculate_default_transform(
         in_raster.crs, new_crs, in_raster.width, in_raster.height, *in_raster.bounds
     )
     kwargs = in_raster.meta.copy()
     kwargs.update({
         'driver':'GTiff',
         'crs': new_crs,
         'transform': affine,
         'affine': affine,
         'width': width,
         'height': height
     })
     return affine, height, width, kwargs

    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 test_calculate_default_transform():
    target_transform = Affine(
        0.0028956983577810586, 0.0, -78.95864996545055,
        0.0, -0.0028956983577810586, 25.550873767433984
    )
    with rasterio.drivers():
        with rasterio.open('tests/data/RGB.byte.tif') as src:
            l, b, r, t = src.bounds
            wgs84_crs = {'init': 'EPSG:4326'}
            dst_transform, width, height = calculate_default_transform(
                l, b, r, t, src.width, src.height, src.crs, wgs84_crs)

            assert dst_transform.almost_equals(target_transform)
            assert width == 824
            assert height == 686

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, {'init': '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 get_zoom(input, dst_crs="EPSG:3857"):
    input = input.replace("s3://", "/vsicurl/http://s3.amazonaws.com/")
    with rasterio.drivers():
        with rasterio.open(input) as src:
            # Compute the geographic bounding box of the dataset.
            (west, east), (south, north) = transform(
                src.crs, "EPSG:4326", src.bounds[::2], src.bounds[1::2])

            affine, _, _ = calculate_default_transform(src.crs, dst_crs,
                src.width, src.height, *src.bounds, resolution=None)

            # grab the lowest resolution dimension
            resolution = max(abs(affine[0]), abs(affine[4]))

            return int(round(math.log((2 * math.pi * 6378137) /
                                      (resolution * CHUNK_SIZE)) / math.log(2)))

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/mapbox/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 = {"init": "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 to_srs(self, srs, resolution=None, src_nodata=None, dst_nodata=None,
               resampling=Resampling.nearest):
        affine, width, height = calculate_default_transform(self.crs,
                                                            srs,
                                                            self.shape[1],
                                                            self.shape[0],
                                                            *self.bounds,
                                                            resolution=resolution)
        if dst_nodata is None:
            dst_nodata = self.fill_value

        destination = RectifiedGrid(np.zeros([height, width], self.dtype),
                                    srs,
                                    affine,
                                    fill_value=dst_nodata)
        return self.to_srs_like(destination, src_nodata,
                                dst_nodata, resampling)

def test_calculate_default_transform_multiple_resolutions():
    with rasterio.drivers():
        with rasterio.open('tests/data/RGB.byte.tif') as src:
            l, b, r, t = src.bounds
            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(
                l, b, r, t, src.width, src.height, src.crs,
                {'init': 'EPSG:4326'}, resolution=target_resolution
            )

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

def test_project():
    data = np.array([[0,0,1,0,0],
                     [0,0,1,0,0],
                     [1,1,1,1,1],
                     [0,0,1,0,0],
                     [0,0,1,0,0]],dtype=np.int32)
    geodict = {'xmin':50,'xmax':50.4,'ymin':50,'ymax':50.4,'dx':0.1,'dy':0.1,'nx':5,'ny':5}
    gd = GeoDict(geodict)
    grid = GDALGrid(data,gd)
    projstr = "+proj=utm +zone=40 +north +ellps=WGS84 +datum=WGS84 +units=m +no_defs "
    newgrid = grid.project(projstr,method='nearest')

    try:
        tdir = tempfile.mkdtemp()
        outfile = os.path.join(tdir,'output.bil')
        grid.save(outfile)
        with rasterio.open(outfile) as src:
            aff = src.transform
            data = src.read(1)
            src_crs = CRS().from_string(GeoDict.DEFAULT_PROJ4).to_dict()
            dst_crs = CRS().from_string(projstr).to_dict()
            nrows,ncols = data.shape
            left = aff.xoff
            top = aff.yoff
            right,bottom = aff * (ncols-1, nrows-1)
            dst_transform,width,height = calculate_default_transform(src_crs,dst_crs,
                                                                     ncols,nrows,
                                                                     left,bottom,
                                                                     right,top)
            destination = np.zeros((height,width))
            reproject(data,
                      destination,
                      src_transform=aff,
                      src_crs=src_crs,
                      dst_transform=dst_transform,
                      dst_crs=dst_crs,
                      src_nodata=src.nodata,
                      dst_nodata=np.nan,
                      resampling=Resampling.nearest)
            x = 1
    except:
        pass
    finally:
        shutil.rmtree(tdir)