def test_warp_from_to_file(tmpdir):
    """File to file"""
    tiffname = str(tmpdir.join('foo.tif'))
    with rasterio.open('rasterio/tests/data/RGB.byte.tif') as src:
        dst_transform = [-8789636.708, 300.0, 0.0, 2943560.235, 0.0, -300.0]
        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 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 reproject_dataset(geotiff_path):
    """Project a GeoTIFF to the WGS84 coordinate reference system.
    See https://mapbox.github.io/rasterio/topics/reproject.html"""

    # We want to project the GeoTIFF coordinate reference system (crs)
    # to WGS84 (e.g. into the familiar Lat/Lon pairs). WGS84 is analogous
    # to EPSG:4326
    dst_crs = 'EPSG:4326'

    with rasterio.open(geotiff_path) as src:
        transform, width, height = rasterio.warp.calculate_default_transform(
            src.crs, dst_crs, src.width, src.height, *src.bounds)
        kwargs = src.meta.copy()
        kwargs.update({
            'crs': dst_crs,
            'transform': transform,
            'width': width,
            'height': height
        })

        satellite_img_name = get_file_name(geotiff_path)
        out_file_name = "{}_wgs84.tif".format(satellite_img_name)
        out_path = os.path.join(WGS84_DIR, out_file_name)
        with rasterio.open(out_path, '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.transform,
                    src_crs=src.crs,
                    dst_transform=transform,
                    dst_crs=dst_crs,
                    resampling=rasterio.warp.Resampling.nearest)

        return rasterio.open(out_path), out_path

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)

def get_bands(inputs, d, i=None):
    """Get a rasterio.Band object from calc's inputs"""
    path = inputs[d] if d in dict(inputs) else inputs[int(d)-1][1]
    if i:
        return rasterio.band(rasterio.open(path), i)
    else:
        src = rasterio.open(path)
        return [rasterio.band(src, i) for i in src.indexes]

    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 __reproject(self, in_raster, out_raster, affine, new_crs):
     for k in range(1, in_raster.count + 1):
         reproject(
             source=rasterio.band(in_raster, k),
             destination=rasterio.band(out_raster, k),
             src_transform=affine,
             src_crs=in_raster.crs,
             dst_transform=affine,
             dst_crs=new_crs,
             resampling=RESAMPLING.nearest)
     return out_raster

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_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_band():
    with rasterio.open('tests/data/RGB.byte.tif') as src:
        b = rasterio.band(src, 1)
        assert b.ds == src
        assert b.bidx == 1
        assert b.dtype in src.dtypes
        assert b.shape == src.shape

def test_warp_from_file():
    """File to ndarray"""
    with rasterio.open('rasterio/tests/data/RGB.byte.tif') as src:
        dst_transform = [-8789636.708, 300.0, 0.0, 2943560.235, 0.0, -300.0]
        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 = numpy.empty(src.shape, dtype=numpy.uint8)
        reproject(
            rasterio.band(src, 1), 
            destin, 
            dst_transform=dst_transform, 
            dst_crs=dst_crs)
    assert destin.any()
    try:
        import matplotlib.pyplot as plt
        plt.imshow(destin)
        plt.gray()
        plt.savefig('test_warp_from_filereproject.png')
    except:
        pass

def process_tile(tile):
    """Process a single MBTiles tile."""
    global base_kwds, src
    # Get the bounds of the tile.
    ulx, uly = mercantile.xy(*mercantile.ul(tile.x, tile.y, tile.z))
    lrx, lry = mercantile.xy(*mercantile.ul(tile.x + 1, tile.y + 1, tile.z))

    kwds = base_kwds.copy()
    kwds["transform"] = from_bounds(ulx, lry, lrx, uly, 256, 256)

    with rasterio.open("/vsimem/tileimg", "w", **kwds) as tmp:
        # Reproject the src dataset into image tile.
        for bidx in tmp.indexes:
            reproject(rasterio.band(src, bidx), rasterio.band(tmp, bidx))

    # Get contents of the virtual file.
    contents = bytearray(virtual_file_to_buffer("/vsimem/tileimg"))
    return tile, contents

def test_shapes_band_shortcut():
    """Access to shapes of labeled features"""
    with rasterio.drivers():
        with rasterio.open('rasterio/tests/data/shade.tif') as src:
            shapes = ftrz.shapes(rasterio.band(src, 1))
            shape, val = next(shapes)
            assert shape['type'] == 'Polygon'
            assert len(shape['coordinates']) == 1
            assert val == 255

def process_chunk(tile, input, creation_options, resampling):
    """Process a single tile."""

    from rasterio.warp import RESAMPLING

    input = input.replace("s3://", "/vsicurl/http://s3.amazonaws.com/")

    print("Chunking initial image for", tile)

    # Get the bounds of the tile.
    ulx, uly = mercantile.xy(
        *mercantile.ul(tile.x, tile.y, tile.z))
    lrx, lry = mercantile.xy(
        *mercantile.ul(tile.x + 1, tile.y + 1, tile.z))

    tmp_path = "/vsimem/tile"

    with rasterio.drivers():
        with rasterio.open(input, "r") as src:
            meta = src.meta.copy()
            meta.update(creation_options)
            meta["height"] = CHUNK_SIZE
            meta["width"] = CHUNK_SIZE
            meta["transform"] = from_bounds(ulx, lry, lrx, uly, CHUNK_SIZE, CHUNK_SIZE)

            # write to a tmp file to allow GDAL to handle the transform
            with rasterio.open(tmp_path, "w", **meta) as tmp:
                # Reproject the src dataset into image tile.
                for bidx in src.indexes:
                    reproject(
                        source=rasterio.band(src, bidx),
                        destination=rasterio.band(tmp, bidx),
                        resampling=getattr(RESAMPLING, resampling),
                        num_threads=multiprocessing.cpu_count(),
                    )

                # check for chunks containing only NODATA
                data = tmp.read(masked=True)

            if data.mask.all():
                return

            # TODO hard-coded for the first band
            return (tile, data[0])

def test_shapes_band(pixelated_image, pixelated_image_file):
    """Shapes from a band should match shapes from an array."""
    truth = list(shapes(pixelated_image))

    with rasterio.open(pixelated_image_file) as src:
        band = rasterio.band(src, 1)
        assert truth == list(shapes(band))

        # Mask band should function, but will mask out some results
        assert truth[0] == list(shapes(band, mask=band))[0]

def process_chunk_task(task):
    """
    Chunks the image into tile_dim x tile_dim tiles,
    and saves them to the target folder (s3 or local)

    Returns the extent of the output raster.
    """

    creation_options = {
        "driver": "GTiff",
        "crs": "EPSG:3857",
        "tiled": True,
        "compress": "deflate",
        "predictor":   2, # 3 for floats, 2 otherwise
        "sparse_ok": True
    }

    with rasterio.open(task.source_uri, "r") as src:
        meta = src.meta.copy()
        meta.update(creation_options)
        meta.update(task.target_meta)
        
        cols = meta["width"]
        rows = meta["height"]

        tmp_path = "/vsimem/" + get_filename(task.target)

        with rasterio.open(tmp_path, "w", **meta) as tmp:
            # Reproject the src dataset into image tile.
            warped = []
            for bidx in src.indexes:
                source = rasterio.band(src, bidx)
                warped.append(numpy.zeros((cols, rows), dtype=meta['dtype']))

                warp.reproject(
                    source=source,
                    src_nodata=0,
                    destination=warped[bidx - 1],
                    dst_transform=meta["transform"],
                    dst_crs=meta["crs"],
                    # resampling=RESAMPLING.bilinear
                )

            # check for chunks containing only zero values
            if not any(map(lambda b: b.any(), warped)):
                return False

            # write out our warped data to the vsimem raster
            for bidx in src.indexes:
                tmp.write_band(bidx, warped[bidx - 1])

    contents = bytearray(virtual_file_to_buffer(tmp_path))

    write_bytes_to_target(task.target, contents)
    return True

def project_raster(src_raster, dst_raster, dst_crs):
    """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' 
    """
    try:
        import rasterio
        from rasterio.warp import calculate_default_transform, reproject, RESAMPLING
    except:
        print('This function requires rasterio.')

    with rasterio.open(src_raster) as src:
        affine, width, height = calculate_default_transform(
            src.crs, dst_crs, src.width, src.height, *src.bounds)
        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.nearest)

def project2wgs(gtiff, output):
    with rio.Env():
        with rio.open(gtiff) as src:
            out_kwargs = src.meta.copy()
            out_kwargs['driver'] = 'GTiff'

            print(out_kwargs)

            res = (0.01, 0.01)
            dst_crs = crs.from_string('+units=m +init=epsg:4326')

            #dst_width, dst_height = src.width, src.height
            xmin, ymin, xmax, ymax = [-127.8294048826629989,5.1830409679864857,
                                      -59.0561278820333229,49.9999999955067977]
            dst_transform = Affine(res[0], 0, xmin, 0, -res[1], ymax)
            dst_width = max(int(ceil((xmax - xmin) / res[0])), 1)
            dst_height = max(int(ceil((ymax - ymin) / res[1])), 1)
            print(dst_transform)

            out_kwargs.update({
                    'crs': dst_crs,
                    'transform': dst_transform,
                    'affine': dst_transform,
                    'width': dst_width,
                    'height': dst_height
                })

            print(out_kwargs)

            with rio.open(output, 'w', **out_kwargs) as dst:
                reproject(
                    source=rio.band(src, 1),
                    destination=rio.band(dst, 1),
                    src_transform=src.affine,
                    src_crs=src.crs,
                    src_nodata=src.nodata,
                    dst_transform=out_kwargs['transform'],
                    dst_crs=out_kwargs['crs'],
                    dst_nodata=src.nodata,
                    resampling=0,
                    num_threads=2)

def ConvertRaster2LatLong(InputRasterFile,OutputRasterFile):

    """
    Convert a raster to lat long WGS1984 EPSG:4326 coordinates for global plotting

    MDH

    """

    # import modules
    import rasterio
    from rasterio.warp import reproject, calculate_default_transform as cdt, Resampling

    # read the source raster
    with rasterio.open(InputRasterFile) as src:
        #get input coordinate system
        Input_CRS = src.crs
        # define the output coordinate system
        Output_CRS = {'init': "epsg:4326"}
        # set up the transform
        Affine, Width, Height = cdt(Input_CRS,Output_CRS,src.width,src.height,*src.bounds)
        kwargs = src.meta.copy()
        kwargs.update({
            'crs': Output_CRS,
            'transform': Affine,
            'affine': Affine,
            'width': Width,
            'height': Height
        })

        with rasterio.open(OutputRasterFile, '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=Output_CRS,
                    resampling=Resampling.bilinear)

def run(input_file, out_file, dst_crs):
    try:
        print("Starting proccess...")

        createOutFolder(out_file)

        with rasterio.open(input_file) as src:

            affine, width, height = calculate_default_transform(
                src.crs, dst_crs, src.width, src.height,*src.bounds)
            kwargs = src.meta.copy()
            kwargs.update({
                'crs': dst_crs,
                'transform': affine,
                'affine': affine,
                'width': width,
                'height': height,
                'compress': 'lzw',
                'nodata': 0
            })

            with rasterio.open(out_file, 'w', **kwargs) as dst:

                reproject(
                    source=rasterio.band(src, 1),
                    destination=rasterio.band(dst, 1),
                    src_transform=src.affine,
                    src_crs=src.crs,
                    dst_transform=affine,
                    dst_crs=dst_crs,
                    resampling=RESAMPLING.nearest)

                for i in dst.indexes:
                    band_dst = dst.read(i) / 100
                    dst.write_band(i, band_dst)

        print("Successfully finished proccess!")

    except Exception as error:
        print("Error creating out folder: {}".format(error))