def get_resized_8chan_image_test(image_id, datapath, bs_rgb, bs_mul):
    im = []

    fn = get_test_image_path_from_imageid(image_id, datapath)
    with rasterio.open(fn, 'r') as f:
        values = f.read().astype(np.float32)
        for chan_i in range(3):
            min_val = bs_rgb[chan_i]['min']
            max_val = bs_rgb[chan_i]['max']
            values[chan_i] = np.clip(values[chan_i], min_val, max_val)
            values[chan_i] = (values[chan_i] - min_val) / (max_val - min_val)
            im.append(skimage.transform.resize(
                values[chan_i],
                (INPUT_SIZE, INPUT_SIZE)))

    fn = get_test_image_path_from_imageid(image_id, datapath, mul=True)
    with rasterio.open(fn, 'r') as f:
        values = f.read().astype(np.float32)
        usechannels = [1, 2, 5, 6, 7]
        for chan_i in usechannels:
            min_val = bs_mul[chan_i]['min']
            max_val = bs_mul[chan_i]['max']
            values[chan_i] = np.clip(values[chan_i], min_val, max_val)
            values[chan_i] = (values[chan_i] - min_val) / (max_val - min_val)
            im.append(skimage.transform.resize(
                values[chan_i],
                (INPUT_SIZE, INPUT_SIZE)))

    im = np.array(im)  # (ch, w, h)
    im = np.swapaxes(im, 0, 2)  # -> (h, w, ch)
    im = np.swapaxes(im, 0, 1)  # -> (w, h, ch)
    return im

def tiffs(tmpdir):
    with rasterio.open('tests/data/RGB.byte.tif') as src:
        profile = src.profile

        shadowed_profile = profile.copy()
        shadowed_profile['count'] = 4
        with rasterio.open(
                str(tmpdir.join('shadowed.tif')), 'w',
                **shadowed_profile) as dst:

            for i, band in enumerate(src.read(masked=False), 1):
                dst.write(band, i)
            dst.write(band, 4)

        del profile['nodata']
        with rasterio.open(
                str(tmpdir.join('no-nodata.tif')), 'w',
                **profile) as dst:
            dst.write(src.read(masked=False))

        with rasterio.open(
                str(tmpdir.join('sidecar-masked.tif')), 'w',
                **profile) as dst:
            dst.write(src.read(masked=False))
            mask = np.zeros(src.shape, dtype='uint8')
            dst.write_mask(mask)

    return tmpdir

def test_warp_reproject_like(runner, tmpdir):
    likename = str(tmpdir.join('like.tif'))
    kwargs = {
        "crs": {'init': 'epsg:4326'},
        "transform": affine.Affine(0.001, 0, -106.523,
                                   0, -0.001, 39.6395),
        "count": 1,
        "dtype": rasterio.uint8,
        "driver": "GTiff",
        "width": 10,
        "height": 10,
        "nodata": 0
    }

    with rasterio.open(likename, 'w', **kwargs) as dst:
        data = np.zeros((10, 10), dtype=rasterio.uint8)
        dst.write(data, indexes=1)

    srcname = 'tests/data/shade.tif'
    outputname = str(tmpdir.join('test.tif'))
    result = runner.invoke(main_group, [
        'warp', srcname, outputname, '--like', likename])
    assert result.exit_code == 0
    assert os.path.exists(outputname)

    with rasterio.open(outputname) as output:
        assert output.crs == {'init': 'epsg:4326'}
        assert np.allclose(
            [0.001, 0.001], [output.transform.a, -output.transform.e])
        assert output.width == 10
        assert output.height == 10

def test_mask_crop(runner, tmpdir):
    output = str(tmpdir.join('test.tif'))

    with rasterio.open('tests/data/shade.tif') as src:

        result = runner.invoke(
            features.mask,
            [
                'tests/data/shade.tif', output,
                '--crop',
                '--geojson-mask', '-'
            ],
            input=TEST_MERC_FEATURES
        )
        assert result.exit_code == 0
        assert os.path.exists(output)
        with rasterio.open(output) as out:
            assert out.shape[1] == src.shape[1]
            assert out.shape[0] < src.shape[0]
            assert out.shape[0] == 824

    # Adding invert option after crop should be ignored
    result = runner.invoke(
        features.mask,
        [
            'tests/data/shade.tif', output,
            '--crop',
            '--invert',
            '--geojson-mask', '-'
        ],
        input=TEST_MERC_FEATURES
    )
    assert result.exit_code == 0
    assert 'Invert option ignored' in result.output

def test_mask_crop(runner, tmpdir, basic_feature, pixelated_image):
    """
    In order to test --crop option, we need to use a transform more similar to
    a normal raster, with a negative y pixel size.
    """

    image = pixelated_image
    outfilename = str(tmpdir.join('pixelated_image.tif'))
    kwargs = {
        "crs": CRS({'init': 'epsg:4326'}),
        "transform": Affine(1, 0, 0, 0, -1, 0),
        "count": 1,
        "dtype": rasterio.uint8,
        "driver": "GTiff",
        "width": image.shape[1],
        "height": image.shape[0],
        "nodata": 255}
    with rasterio.open(outfilename, 'w', **kwargs) as out:
        out.write(image, indexes=1)

    output = str(tmpdir.join('test.tif'))

    truth = np.zeros((3, 3))
    truth[0:2, 0:2] = 1

    result = runner.invoke(
        main_group,
        ['mask', outfilename, output, '--crop', '--geojson-mask', '-'],
        input=json.dumps(basic_feature))
    assert result.exit_code == 0
    assert os.path.exists(output)
    with rasterio.open(output) as out:
        assert np.array_equal(
            truth,
            out.read(1, masked=True).filled(0))

def test_update_band(tmpdir):
    tiffname = str(tmpdir.join('foo.tif'))
    shutil.copy('rasterio/tests/data/RGB.byte.tif', tiffname)
    with rasterio.open(tiffname, 'r+') as f:
        f.write_band(1, numpy.zeros(f.shape, dtype=f.dtypes[0]))
    with rasterio.open(tiffname) as f:
        assert not f.read_band(1).any()

def test_tags_update(tmpdir):
    tiffname = str(tmpdir.join('foo.tif'))
    with rasterio.open(
            tiffname, 
            'w', 
            driver='GTiff', 
            count=1, 
            dtype=rasterio.uint8, 
            width=10, 
            height=10) as dst:

        dst.update_tags(a='1', b='2')
        dst.update_tags(1, c=3)
        with pytest.raises(IndexError):
            dst.update_tags(4, d=4)

        assert dst.tags() == {'a': '1', 'b': '2'}
        assert dst.tags(1) == {'c': '3' }
        
        # Assert that unicode tags work.
        # Russian text appropriated from pytest issue #319
        # https://bitbucket.org/hpk42/pytest/issue/319/utf-8-output-in-assertion-error-converted
        dst.update_tags(ns='rasterio_testing', rus=u'другая строка')
        assert dst.tags(ns='rasterio_testing') == {'rus': u'другая строка'}

    with rasterio.open(tiffname) as src:
        assert src.tags() == {'a': '1', 'b': '2'}
        assert src.tags(1) == {'c': '3'}
        assert src.tags(ns='rasterio_testing') == {'rus': u'другая строка'}

def test_write_read_gcps(tmpdir):
    tiffname = str(tmpdir.join('test.tif'))
    gcps = [GroundControlPoint(1, 1, 100.0, 1000.0, z=0.0)]

    with rasterio.open(tiffname, 'w', driver='GTiff', dtype='uint8', count=1,
                       width=10, height=10, crs='epsg:4326', gcps=gcps) as dst:
        pass

    with rasterio.open(tiffname, 'r+') as dst:
        gcps, crs = dst.gcps
        assert crs['init'] == 'epsg:4326'
        assert len(gcps) == 1
        point = gcps[0]
        assert (1, 1) == (point.row, point.col)
        assert (100.0, 1000.0, 0.0) == (point.x, point.y, point.z)

        dst.gcps = [
            GroundControlPoint(1, 1, 100.0, 1000.0, z=0.0),
            GroundControlPoint(2, 2, 200.0, 2000.0, z=0.0)], crs

        gcps, crs = dst.gcps

        assert crs['init'] == 'epsg:4326'
        assert len(gcps) == 2
        point = gcps[1]
        assert (2, 2) == (point.row, point.col)
        assert (200.0, 2000.0, 0.0) == (point.x, point.y, point.z)

    def __next__(self):
        length = len(self)
        if length == 0:
            raise StopIteration()
        if self._curent_position == length:
            if self._loop:
                self.reset()
            else:
                raise StopIteration()

        entry = self._datasets[self._curent_position]
        env = getattr(self, 'rasterio_env', {})
        self._curent_position += 1
        entry_name, entry_components = entry
        new_components = {}
        cache_data = self._cache_data
        use_tensorflow_io = False
        for component_name, component_path in entry_components.items():
            if isinstance(component_path, DatasetReader):
                component_path = component_path.name
            local_component_path = component_path
            url_components = urlparse(component_path)
            if not url_components.scheme:
                cache_data = False
                if url_components.path.startswith('/vsigs/'):
                    cache_data = True  # We should check if we run inside GCP ML Engine
                    use_tensorflow_io = True
                    component_path = url_components.path[6:]
                    component_path = "gs:/" + component_path
            else:
                if url_components.scheme == 'file':
                    local_component_path = url_components.path
                    use_tensorflow_io = False
                    cache_data = False

            with rasterio.Env(**env):
                if use_tensorflow_io:
                    real_path = component_path
                    data = IOUtils.open_file(real_path, "rb").read()
                    if cache_data:
                        hash = sha224(component_path.encode("utf8")).hexdigest()
                        hash_part = "/".join(list(hash)[:3])
                        dataset_path = os.path.join(self._temp_dir, hash_part)
                        if not IOUtils.file_exists(dataset_path):
                            IOUtils.recursive_create_dir(dataset_path)
                        dataset_path = os.path.join(dataset_path, os.path.basename(component_path))
                        if not IOUtils.file_exists(dataset_path):
                            f = IOUtils.open_file(dataset_path, "wb")
                            f.write(data)
                            f.close()
                        component_src = rasterio.open(dataset_path)
                    else:
                        with NamedTemporaryFile() as tmpfile:
                            tmpfile.write(data)
                            tmpfile.flush()
                            component_src = rasterio.open(tmpfile.name)
                else:
                    component_src = rasterio.open(local_component_path)
                new_components[component_name] = component_src
        return (entry_name, new_components)

def tiffs(tmpdir):

    data = numpy.ones((1, 1, 1), 'uint8')

    kwargs = {'count': '1',
              'driver': 'GTiff',
              'dtype': 'uint8',
              'height': 1,
              'width': 1}

    kwargs['transform'] = Affine( 1, 0, 1,
                                  0,-1, 1)
    with rasterio.open(str(tmpdir.join('a-sw.tif')), 'w', **kwargs) as r:
        r.write(data * 40)

    kwargs['transform'] = Affine( 1, 0, 2,
                                  0,-1, 2)
    with rasterio.open(str(tmpdir.join('b-ct.tif')), 'w', **kwargs) as r:
        r.write(data * 60)

    kwargs['transform'] = Affine( 2, 0, 3,
                                  0,-2, 4)
    with rasterio.open(str(tmpdir.join('c-ne.tif')), 'w', **kwargs) as r:
        r.write(data * 90)

    kwargs['transform'] = Affine( 2, 0, 2,
                                  0,-2, 4)
    with rasterio.open(str(tmpdir.join('d-ne.tif')), 'w', **kwargs) as r:
        r.write(data * 120)

    return tmpdir

def test_data_dir_2(tmpdir):
    kwargs = {
        "crs": {'init': 'epsg:4326'},
        "transform": (-114, 0.2, 0, 46, 0, -0.1),
        "count": 1,
        "dtype": rasterio.uint8,
        "driver": "GTiff",
        "width": 10,
        "height": 10
        # these files have undefined nodata.
    }

    with rasterio.drivers():

        with rasterio.open(str(tmpdir.join('b.tif')), 'w', **kwargs) as dst:
            data = numpy.zeros((10, 10), dtype=rasterio.uint8)
            data[0:6, 0:6] = 255
            dst.write_band(1, data)

        with rasterio.open(str(tmpdir.join('a.tif')), 'w', **kwargs) as dst:
            data = numpy.zeros((10, 10), dtype=rasterio.uint8)
            data[4:8, 4:8] = 254
            dst.write_band(1, data)

    return tmpdir

    def run(self):
        out_fname = pjoin(self.out_dir,
                          CONFIG.get('outputs', 'rasterise_filename'))
        ds_list_fname = pjoin(self.out_dir,
                              CONFIG.get('outputs', 'query_filename'))

        with open(ds_list_fname, 'r') as infile:
            ds_list = pickle.load(infile)

        vector_fname = CONFIG.get('work', 'vector_filename')

        img_fname = ds_list[0].datasets[DatasetType.FC25].path
        with rasterio.open(img_fname) as src:
            crs = src.crs
            transform = src.affine
            height = src.height
            width = src.width

        res = rasterise_vector(vector_fname, shape=(height, width),
                               transform=transform, crs=crs)

        kwargs = {'count': 1,
                  'width': width,
                  'height': height,
                  'crs': crs,
                  'transform': transform,
                  'dtype': res.dtype.name,
                  'nodata': 0}

        with rasterio.open(out_fname, 'w', **kwargs) as src:
            src.write(1, res)

        # We could just set the image as the Luigi completion target...
        with self.output().open('w') as outf:
            outf.write('Complete')

def calc_monthly_average(basepath, filename, layers_by_month, epsg="3857"):
    print "-----Averages"

    if not os.path.exists(basepath):
        os.makedirs(basepath)

    for month in layers_by_month:
        output_path = basepath + "/" + filename + "_" + month + "_" + epsg + ".tif"

        data = None
        kwargs = None

        print "Processing: ", str(month)
        for f in layers_by_month[month]:
            print "Reading: ",  f
            r = rasterio.open(f)

            if data is None:
                data, kwargs = initialize_rasterio_raster(r, rasterio.float32)

            data = data + r.read_band(1).astype(float)

        data = data / len(layers_by_month[month])

        # writing
        print "Writing: ", output_path
        with rasterio.open(output_path, 'w', **kwargs) as dst:
            dst.write_band(1, data.astype(rasterio.float32))

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 test_warp_reproject_dst_bounds(runner, tmpdir):
    """--x-dst-bounds option works."""
    srcname = 'tests/data/shade.tif'
    outputname = str(tmpdir.join('test.tif'))
    out_bounds = [-106.45036, 39.6138, -106.44136, 39.6278]
    result = runner.invoke(
        warp.warp, [srcname, outputname, '--dst-crs', 'EPSG:4326',
                    '--res', 0.001, '--x-dst-bounds'] + out_bounds)
    assert result.exit_code == 0
    assert os.path.exists(outputname)

    with rasterio.open(srcname) as src:
        with rasterio.open(outputname) as output:
            assert output.crs == {'init': 'epsg:4326'}
            assert np.allclose(output.bounds[0::3],
                                  [-106.45036, 39.6278])
            assert np.allclose([0.001, 0.001],
                                  [output.affine.a, -output.affine.e])

            # XXX: an extra row and column is produced in the dataset
            # because we're using ceil instead of floor internally.
            # Not necessarily a bug, but may change in the future.
            assert np.allclose([output.bounds[2]-0.001, output.bounds[1]+0.001],
                                  [-106.44136, 39.6138])
            assert output.width == 10
            assert output.height == 15

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 test_read_basic(self):
     with rasterio.open('tests/data/shade.tif') as s:
         a = s.read(masked=True)  # Gray
         self.assertEqual(a.ndim, 3)
         self.assertEqual(a.shape, (1, 1024, 1024))
         self.assertTrue(hasattr(a, 'mask'))
         self.assertEqual(a.fill_value, 255)
         self.assertEqual(list(set(s.nodatavals)), [255])
         self.assertEqual(a.dtype, rasterio.ubyte)
         self.assertEqual(a.sum((1, 2)).tolist(), [0])
     with rasterio.open('tests/data/RGB.byte.tif') as s:
         a = s.read(masked=True)  # RGB
         self.assertEqual(a.ndim, 3)
         self.assertEqual(a.shape, (3, 718, 791))
         self.assertTrue(hasattr(a, 'mask'))
         self.assertEqual(a.fill_value, 0)
         self.assertEqual(list(set(s.nodatavals)), [0])
         self.assertEqual(a.dtype, rasterio.ubyte)
         a = s.read(masked=False)  # no mask
         self.assertFalse(hasattr(a, 'mask'))
         self.assertEqual(list(set(s.nodatavals)), [0])
         self.assertEqual(a.dtype, rasterio.ubyte)
     with rasterio.open('tests/data/float.tif') as s:
         a = s.read(masked=True)  # floating point values
         self.assertEqual(a.ndim, 3)
         self.assertEqual(a.shape, (1, 2, 3))
         self.assert_(hasattr(a, 'mask'))
         self.assertEqual(list(set(s.nodatavals)), [None])
         self.assertEqual(a.dtype, rasterio.float64)

def run( args ):
	# unpack some args being passed into the run wrapper
	x = args['x']
	y = args['y']
	z = args['z']
	meshgrid_10min = args['meshgrid_10min']
	output_filename_10min = args['output_filename_10min']
	meta_10min = args['meta_10min']
	meta_10min.update( compress='lzw' )

	# interpolate to a global (unmasked) 10 min output
	new_grid = interpolate_to_grid( x, y, z, meshgrid_10min, output_filename_10min, meta_10min )
	new_grid = rasterio.open( output_filename_10min ).read( 1 ) # read it back in << this is a hack!

	# unpack some args being passed in
	template_raster = args['template_raster']
	src_affine = args['meta_10min']['affine']
	dst_affine = template_raster.affine
	method = 'cubic_spline'
	output_filename = args['output_filename']
	src_crs = {'init':'epsg:4326'}
	dst_crs = {'init':'epsg:3338'}
	meta_akcan = template_raster.meta
	meta_akcan.update( compress='lzw', crs=dst_crs  )
	mask = args['mask']

	# regrid to 3338 and AKCAN extent -- currently hardwired
	akcan = regrid( new_grid, template_raster.read( 1 ), src_affine, src_crs, dst_affine, dst_crs, method=method )
	with rasterio.open( output_filename, 'w', **meta_akcan ) as out:
		if mask is not None:
			akcan[ mask == 0 ] = meta_akcan[ 'nodata' ]
		out.write( akcan.astype( np.float32 ), 1 ) # watch this hardwired type!
	return output_filename

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 test_blockysize_guard(tmpdir):
    """blockysize can't be greater than image height."""
    tiffname = str(tmpdir.join('foo.tif'))
    with pytest.raises(ValueError):
        profile = default_gtiff_profile.copy()
        profile.update(count=1, width=256, height=128)
        rasterio.open(tiffname, 'w', **profile)