def test_transform_bounds_densify():
    # This transform is non-linear along the edges, so densification produces
    # a different result than otherwise
    src_crs = {'init': 'EPSG:4326'}
    dst_crs = {'init': 'EPSG:32610'}
    assert numpy.allclose(
        transform_bounds(
            src_crs,
            dst_crs,
            -120, 40, -80, 64,
            densify_pts=0
        ),
        (
            646695.227266598, 4432069.056898901,
            4201818.984205882, 7807592.187464975
        )
    )

    assert numpy.allclose(
        transform_bounds(
            src_crs,
            dst_crs,
            -120, 40, -80, 64,
            densify_pts=100
        ),
        (
            646695.2272665979, 4432069.056898901,
            4201818.984205882, 7807592.187464977
        )
    )

def test_transform_bounds_densify_out_of_bounds():
    with pytest.raises(ValueError):
        transform_bounds(
            {'init': 'EPSG:4326'},
            {'init': 'EPSG:32610'},
            -120, 40, -80, 64,
            densify_pts=-10
        )

def test_transform_bounds():
    """CRSError is raised."""
    left, bottom, right, top = (
        -11740727.544603072, 4852834.0517692715, -11584184.510675032,
        5009377.085697309)
    src_crs = 'EPSG:3857'
    dst_crs = {'proj': 'foobar'}
    with pytest.raises(CRSError):
        transform_bounds(src_crs, dst_crs, left, bottom, right, top)

def test_transform_bounds_densify_out_of_bounds():
    with pytest.raises(ValueError):
        transform_bounds(
            {"init": "epsg:4326"},
            {"init": "epsg:32610"},
            -120,
            40,
            -80,
            64,
            densify_pts=-10,
        )

def test_transform_bounds_densify():
    # This transform is non-linear along the edges, so densification produces
    # a different result than otherwise
    src_crs = {"init": "epsg:4326"}
    dst_crs = {"init": "epsg:2163"}
    assert np.allclose(
        transform_bounds(src_crs, dst_crs, -120, 40, -80, 64, densify_pts=0),
        (-1684649.41338, -350356.81377, 1684649.41338, 2234551.18559),
    )

    assert np.allclose(
        transform_bounds(src_crs, dst_crs, -120, 40, -80, 64, densify_pts=100),
        (-1684649.41338, -555777.79210, 1684649.41338, 2234551.18559),
    )

def generate_chunk_tasks(image_source, tile_dim):
    tasks = []
    zoom = image_source.zoom
    (min_col, max_col) = (image_source.tile_bounds[0], image_source.tile_bounds[2])
    (min_row, max_row) = (image_source.tile_bounds[1], image_source.tile_bounds[3])

    for tile_col in range(min_col, min(max_col + 1, 2**zoom)):
        for tile_row in range(min_row, min(max_row + 1, 2**zoom)):
            tile_bounds = mercantile.bounds(tile_col, tile_row, zoom)
            (wm_left, wm_bottom, wm_right, wm_top)  = warp.transform_bounds("EPSG:4326",
                                                                           "EPSG:3857",
                                                                            tile_bounds.west,
                                                                            tile_bounds.south,
                                                                            tile_bounds.east,
                                                                            tile_bounds.north)
            affine = transform.from_bounds(wm_left, wm_bottom, wm_right, wm_top, tile_dim, tile_dim)
            target_meta = { 
                "transform": affine[:6],
                "width": tile_dim,
                "height": tile_dim 
            }

            target = os.path.join(image_source.image_folder, "%d/%d/%d.tif" % (zoom, tile_col, tile_row))
            task = ChunkTask(source_uri=image_source.source_uri,
                             target_meta=target_meta,
                             target=target)

            tasks.append(task)

    return tasks

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({"init": "epsg:4326"},
                         dst_projection_string,
                         left,
                         bottom,
                         right,
                         top),
        (
            1721263.7931814701,
            219684.49332178483,
            2002926.56696663,
            479360.16562217404),
    )

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():
    with rasterio.open('tests/data/RGB.byte.tif') as src:
        l, b, r, t = src.bounds
        assert np.allclose(
            transform_bounds(src.crs, {'init': 'EPSG:4326'}, l, b, r, t),
            (
                -78.95864996545055, 23.564991210854686,
                -76.57492370013823, 25.550873767433984
            )
        )

def get_bounds(path):
    """ Retrun bounds in WGS84 system """

    with rasterio.drivers():
        src = rasterio.open(path)

        return transform_bounds(
                    src.crs,
                    {'init': 'EPSG:4326'},
                    *src.bounds)

    def clip(self):
        """ Clip images based on bounds provided
        Implementation is borrowed from
        https://github.com/brendan-ward/rasterio/blob/e3687ce0ccf8ad92844c16d913a6482d5142cf48/rasterio/rio/convert.py
        """

        self.output("Clipping", normal=True)

        # create new folder for clipped images
        path = check_create_folder(join(self.scene_path, 'clipped'))

        try:
            temp_bands = copy(self.bands)
            temp_bands.append('QA')
            for i, band in enumerate(temp_bands):
                band_name = self._get_full_filename(band)
                band_path = join(self.scene_path, band_name)

                self.output("Band %s" % band, normal=True, color='green', indent=1)
                with rasterio.open(band_path) as src:
                    bounds = transform_bounds(
                        {
                            'proj': 'longlat',
                            'ellps': 'WGS84',
                            'datum': 'WGS84',
                            'no_defs': True
                        },
                        src.crs,
                        *self.bounds
                    )

                    if disjoint_bounds(bounds, src.bounds):
                        bounds = adjust_bounding_box(src.bounds, bounds)

                    window = src.window(*bounds)

                    out_kwargs = src.meta.copy()
                    out_kwargs.update({
                        'driver': 'GTiff',
                        'height': window[0][1] - window[0][0],
                        'width': window[1][1] - window[1][0],
                        'transform': src.window_transform(window)
                    })

                    with rasterio.open(join(path, band_name), 'w', **out_kwargs) as out:
                        out.write(src.read(window=window))

            # Copy MTL to the clipped folder
            copyfile(join(self.scene_path, self.scene + '_MTL.txt'), join(path, self.scene + '_MTL.txt'))

            return path

        except IOError as e:
            exit(e.message, 1)

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, {"init": "epsg:4326"}, l, b, r, t),
            (
                -78.95864996545055,
                23.564991210854686,
                -76.57492370013823,
                25.550873767433984,
            ),
        )

def _grid_datasets(datasets, bounds_override, grid_proj, grid_size):
    tiles = defaultdict(list)
    for dataset in datasets:
        dataset_proj = dataset.crs
        dataset_bounds = dataset.bounds
        bounds = bounds_override or BoundingBox(*transform_bounds(dataset_proj, grid_proj, *dataset_bounds))

        for y in range(int(bounds.bottom // grid_size[1]), int(bounds.top // grid_size[1]) + 1):
            for x in range(int(bounds.left // grid_size[0]), int(bounds.right // grid_size[0]) + 1):
                tile_index = (x, y)
                if _check_intersect(tile_index, grid_size, grid_proj, dataset_bounds, dataset_proj):
                    tiles[tile_index].append(dataset)

    return tiles

        def __call__(self):
            for i, path in enumerate(input):
                with rasterio.open(path) as src:
                    bounds = src.bounds
                    if dst_crs:
                        bbox = transform_bounds(src.crs,
                                                dst_crs, *bounds)
                    elif projection == 'mercator':
                        bbox = transform_bounds(src.crs,
                                                {'init': 'epsg:3857'}, *bounds)
                    elif projection == 'geographic':
                        bbox = transform_bounds(src.crs,
                                                {'init': 'epsg:4326'}, *bounds)
                    else:
                        bbox = bounds

                if precision >= 0:
                    bbox = [round(b, precision) for b in bbox]

                yield {
                    'type': 'Feature',
                    'bbox': bbox,
                    'geometry': {
                        'type': 'Polygon',
                        'coordinates': [[
                            [bbox[0], bbox[1]],
                            [bbox[2], bbox[1]],
                            [bbox[2], bbox[3]],
                            [bbox[0], bbox[3]],
                            [bbox[0], bbox[1]]]]},
                    'properties': {
                        'id': str(i),
                        'title': path,
                        'filename': os.path.basename(path)}}

                self._xs.extend(bbox[::2])
                self._ys.extend(bbox[1::2])

def clip_from_wgs(input_path, output_path, bounds):
    """
        Clips a raster based on WGS coordinates.
        Heavily borrowed from: https://github.com/mapbox/rasterio/blob/master/rasterio/rio/clip.py
        :param input_path The path to the raster to clip out of
        :param output_path The path to place the clipped raster
        :param bounds (left, bottom, right, top) bounds in WGS 84 projection
    """

    wgs_crs = CRS({'init': 'epsg:4326'})  # WGS 84 Projection

    with rasterio.open(input_path, 'r') as src:

        window = src.window(*transform_bounds(wgs_crs, src.crs, *bounds))
        height = window[0][1] - window[0][0]
        width = window[1][1] - window[1][0]
        t = 2048    # Threshold, for if the user selects an area greater than 2048 pixels in size
                    # If so, select the center of the area at size of 2048 by 2048
        if width > t:
            if width % 2 != 0:
                width -= 1
            trim = int((width - t) / 2)
            width = t
            window = (window[0], (window[1][0] + trim, window[1][0] + t + trim))

        if height > t:
            if height % 2 != 0:
                height -= 1
            trim = int((height - t) / 2)
            height = t
            window = ((window[0][0] + trim, window[0][0] + t + trim), window[1])

        out_kwargs = src.meta.copy()
        out_kwargs.update({
            'height': height,
            'width': width,
            'transform': src.window_transform(window),
            'dtype': 'int32',   # SyncroSim needs a signed int32
            'nodata': 0
            })

        with rasterio.open(output_path, 'w', **out_kwargs) as dst:
            dst.write(src.read(1, window=window).astype('int32'), 1)

def test_crs_should_be_set(path_rgb_byte_tif, tmpdir, complex):

    """When ``dst_height``, ``dst_width``, and ``dst_transform`` are set
    :py:class:`rasterio.warp.WarpedVRT` calls ``GDALCreateWarpedVRT()``,
    which requires the caller to then set a projection with
    ``GDALSetProjection()``.

    Permalink to ``GDALCreateWarpedVRT()`` call:

        https://github.com/mapbox/rasterio/blob/1f759e5f67628f163ea2550d8926b91545245712/rasterio/_warp.pyx#L753

    """

    vrt_path = str(tmpdir.join("test_crs_should_be_set.vrt"))

    with rasterio.open(path_rgb_byte_tif) as src:

        dst_crs = "EPSG:4326"
        dst_height = dst_width = 10
        dst_bounds = transform_bounds(src.crs, dst_crs, *src.bounds)

        # Destination transform
        left, bottom, right, top = dst_bounds
        xres = (right - left) / dst_width
        yres = (top - bottom) / dst_height
        dst_transform = affine.Affine(xres, 0.0, left, 0.0, -yres, top)

        # The 'complex' test case hits the affected code path
        vrt_options = {"dst_crs": dst_crs}
        if complex:
            vrt_options.update(
                dst_crs=dst_crs,
                dst_height=dst_height,
                dst_width=dst_width,
                dst_transform=dst_transform,
                resampling=Resampling.nearest,
            )

        with WarpedVRT(src, **vrt_options) as vrt:
            rio_shutil.copy(vrt, vrt_path, driver="VRT")
        with rasterio.open(vrt_path) as src:
            assert src.crs

def test_windows(runner, path_rgb_byte_tif):
    result = runner.invoke(main_group, [
        'blocks', path_rgb_byte_tif])
    assert result.exit_code == 0

    fc = json.loads(result.output)

    with rasterio.open(path_rgb_byte_tif) as src:
        block_windows = tuple(src.block_windows())

        # Check the coordinates of the first output feature
        actual_first = fc['features'][0]
        expected_first = block_windows[0]
        bounds = src.window_bounds(expected_first[1])
        xmin, ymin, xmax, ymax = transform_bounds(src.crs, 'EPSG:4326', *bounds)
        coordinates = [[xmin, ymin], [xmin, ymax], [xmax, ymax], [xmax, ymin]]
        assert np.array_equal(
            actual_first['geometry']['coordinates'][0], coordinates)

        assert len(fc['features']) == len(block_windows)
        assert check_features_block_windows(fc['features'], src, bidx=1)

def test_windows_sequence(runner, path_rgb_byte_tif):
    result = runner.invoke(main_group, [
        'blocks',
        path_rgb_byte_tif,
        '--sequence'])
    assert result.exit_code == 0

    features = tuple(map(json.loads, result.output.splitlines()))

    with rasterio.open(path_rgb_byte_tif) as src:

        # Check the coordinates of the first output feature
        actual_first = features[0]
        expected_first = next(src.block_windows())
        bounds = src.window_bounds(expected_first[1])
        xmin, ymin, xmax, ymax = transform_bounds(src.crs, 'EPSG:4326', *bounds)
        coordinates = [[xmin, ymin], [xmin, ymax], [xmax, ymax], [xmax, ymin]]
        assert np.array_equal(
            actual_first['geometry']['coordinates'][0], coordinates)

        assert check_features_block_windows(features, src, bidx=1)

def _make_tiles(bbox, src_crs, minz, maxz):
    '''
    Given a bounding box, zoom range, and source crs,
    find all tiles that would intersect

    Parameters
    -----------
    bbox: list
        [w, s, e, n] bounds
    src_crs: str
        the source crs of the input bbox
    minz: int
        minumum zoom to find tiles for
    maxz: int
        maximum zoom to find tiles for

    Returns
    --------
    tiles: generator
        generator of [x, y, z] tiles that intersect
        the provided bounding box
    '''
    w, s, e, n = transform_bounds(*[src_crs, 'epsg:4326'] + bbox, densify_pts=0)

    EPSILON = 1.0e-10

    w += EPSILON
    s += EPSILON
    e -= EPSILON
    n -= EPSILON

    for z in range(minz, maxz + 1):
        for x, y in _tile_range(
            mercantile.tile(w, n, z),
            mercantile.tile(e, s, z)):

            yield [x, y, z]

def clip(
        ctx,
        files,
        output,
        bounds,
        like,
        driver,
        creation_options):
    """Clips a raster using bounds input directly or from a template raster.

    \b
      $ rio clip input.tif output.tif --bounds xmin ymin xmax ymax
      $ rio clip input.tif output.tif --like template.tif

    If using --bounds, values must be in coordinate reference system of input.
    If using --like, bounds will automatically be transformed to match the
    coordinate reference system of the input.

    It can also be combined to read bounds of a feature dataset using Fiona:

    \b
      $ rio clip input.tif output.tif --bounds $(fio info features.shp --bounds)

    """

    from rasterio.warp import transform_bounds

    verbosity = (ctx.obj and ctx.obj.get('verbosity')) or 1

    with rasterio.drivers(CPL_DEBUG=verbosity > 2):

        output, files = resolve_inout(files=files, output=output)
        input = files[0]

        with rasterio.open(input) as src:
            if bounds:
                if disjoint_bounds(bounds, src.bounds):
                    raise click.BadParameter('must overlap the extent of '
                                             'the input raster',
                                             param='--bounds',
                                             param_hint='--bounds')
            elif like:
                with rasterio.open(like) as template_ds:
                    bounds = template_ds.bounds
                    if template_ds.crs != src.crs:
                        bounds = transform_bounds(template_ds.crs, src.crs,
                                                  *bounds)

                    if disjoint_bounds(bounds, src.bounds):
                        raise click.BadParameter('must overlap the extent of '
                                                 'the input raster',
                                                 param='--like',
                                                 param_hint='--like')

            else:
                raise click.UsageError('--bounds or --like required')

            window = src.window(*bounds)

            out_kwargs = src.meta.copy()
            out_kwargs.update({
                'driver': driver,
                'height': window[0][1] - window[0][0],
                'width': window[1][1] - window[1][0],
                'transform': src.window_transform(window)
            })
            out_kwargs.update(**creation_options)

            with rasterio.open(output, 'w', **out_kwargs) as out:
                out.write(src.read(window=window))