def interp_ds( anom, base, src_crs, src_nodata, dst_nodata, src_transform, resample_type='bilinear',*args, **kwargs ):
		'''	
		anom = [numpy.ndarray] 2-d array representing a single monthly timestep of the data to be downscaled. 
								Must also be representative of anomalies.
		base = [str] filename of the corresponding baseline monthly file to use as template and downscale 
								baseline for combining with anomalies.
		src_transform = [affine.affine] 6 element affine transform of the input anomalies. [should be greenwich-centered]
		resample_type = [str] one of ['bilinear', 'count', 'nearest', 'mode', 'cubic', 'index', 'average', 'lanczos', 'cubic_spline']
		'''	
		import rasterio
		from rasterio.warp import reproject, RESAMPLING

		resampling = {'average':RESAMPLING.average,
					'cubic':RESAMPLING.cubic,
					'lanczos':RESAMPLING.lanczos,
					'bilinear':RESAMPLING.bilinear,
					'cubic_spline':RESAMPLING.cubic_spline,
					'mode':RESAMPLING.mode,
					'count':RESAMPLING.count,
					'index':RESAMPLING.index,
					'nearest':RESAMPLING.nearest }
		
		base = rasterio.open( base )
		baseline_arr = base.read( 1 )
		baseline_meta = base.meta
		baseline_meta.update( compress='lzw' )
		output_arr = np.empty_like( baseline_arr )
		
		reproject( anom, output_arr, src_transform=src_transform, src_crs=src_crs, src_nodata=src_nodata, \
				dst_transform=baseline_meta['affine'], dst_crs=baseline_meta['crs'],\
				dst_nodata=dst_nodata, resampling=resampling[ resample_type ], SOURCE_EXTRA=1000 )
		return output_arr

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_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 downscale_wrapper( arr, affine, crs, baseline, output_filename, downscaling_operation, post_downscale_function ):
	# rotate
	if ( self.data.ds.lon > 200.0 ).any() == True:
		dat, lons = utils.shiftgrid( 180., self.data.anomalies, self.historical.ds.lon, start=False )
		a,b,c,d,e,f,g,h,i = affine #flip it to the greenwich-centering
		src_transform = affine.Affine( a, b, -180.0, d, e, 180.0 )
	else:
		dat, lons = ( self.data.ds, self.historical.ds.lon )
		src_transform = affine
	
	# reproject / resample
	src_crs = {'init':'epsg:4326'}
	src_nodata = None # DangerTown™
	baseline_meta = baseline.meta
	baseline_meta.update( compress='lzw' )
	output_arr = np.empty_like( baeline.read( 1 ) )

	# TODO: make this function available for manipulation if used for different needs
	reproject( arr, output_arr, src_transform=src_transform, src_crs=src_crs, src_nodata=src_nodata, \
			dst_transform=baseline_meta['affine'], dst_crs=baseline_meta['crs'],\
			dst_nodata=None, resampling=RESAMPLING.cubic_spline, SOURCE_EXTRA=1000 )
	
	# downscale
	return utils.downscale( arr, output_arr, output_filename, downscaling_operation, \
					baseline_meta, post_downscale_function, mask=None, mask_value=0 )

def test_reproject_dst_nodata_default():
    """
    If nodata is not provided, destination will be filled with 0
    instead of nodata
    """

    params = default_reproject_params()

    with rasterio.drivers():
        source = numpy.ones((params.width, params.height), dtype=numpy.uint8)
        out = numpy.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() == 4461
        assert (out == 0).sum() == (params.dst_width *
                                    params.dst_height - 4461)

def test_warp_from_file():
    """File to ndarray"""
    with rasterio.open('tests/data/RGB.byte.tif') as src:
        dst_transform = affine.Affine.from_gdal(-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 warp_raster_to_template(input_raster_file, template_raster_file, output_file):
    print "reprojecting raster"
    with rasterio.open(input_raster_file) as input_r:
        with rasterio.open(template_raster_file) as template_r:
            input_array = input_r.read(1)
            dest_array = np.zeros((template_r.width, template_r.height), dtype=input_array.dtype)

    reproject(
        input_array,
        dest_array,
        src_transform=input_r.affine,
        src_crs=input_r.crs,
        dst_transform=template_r.affine,
        dst_crs=template_r.crs,
        resampling=Resampling.nearest)

    print "now saving raster"

    with rasterio.open(output_file, 'w',
                  driver = "GTiff",
                  width=template_r.width,
                  height=template_r.height,
                  count=1,
                  dtype=dest_array.dtype,
                  crs=template_r.crs,
                  transform=template_r.affine,
                  nodata=0
                 ) as output:
        output.write(dest_array, indexes=1)

def test_reproject_multi():
    """Ndarry to ndarray"""
    with rasterio.drivers():
        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 = numpy.empty(source.shape, dtype=numpy.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_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,
        Resampling.mode: 437298,
        Resampling.max: 439464,
        Resampling.min: 436397,
        Resampling.med: 437194,
        Resampling.q1: 436397,
        Resampling.q3: 438948
    }

    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={'init': 'EPSG:3857'},
        resampling=method)

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

	def interp_ds( anom, base, output_filename, src_transform, downscaling_operation, post_downscale_function=None, mask=None, mask_value=0 ):
		'''	
		anom = [numpy.ndarray] 2-d array representing a single monthly timestep of the data to be downscaled. Must also be representative of anomalies.
		base = [str] filename of the corresponding baseline monthly file to use as template and downscale baseline for combining with anomalies.
		output_filename = [str] path to the output file to be created following downscaling
		src_transform = [affine.affine] 6 element affine transform of the input anomalies. [should be greenwich-centered]
		downscaling_operation = [str] one of 'add' or 'mult' depending on absolute or relative delta downscaling.
		post_downscale_function = [function] function that takes as input a single 2-d array and returns a 2-d array in the same shape as the input.  
		mask = [numpy.ndarray] 2-d array showing what values should be masked. Masked=0, unmasked=1. must be same shape as base.
		mask_value = [int] what value to use as the masked values. default=0.

		'''		
		from rasterio.warp import reproject, RESAMPLING
		# reproject / resample
		src_crs = {'init':'epsg:4326'}
		src_nodata = None # DangerTown™
		base = rasterio.open( base )
		baseline_arr = base.read( 1 )
		baseline_meta = base.meta
		baseline_meta.update( compress='lzw' )
		output_arr = np.empty_like( baseline_arr )

		# TODO: make this function available for manipulation if used for different needs
		reproject( anom, output_arr, src_transform=src_transform, src_crs=src_crs, src_nodata=src_nodata, \
				dst_transform=baseline_meta['affine'], dst_crs=baseline_meta['crs'],\
				dst_nodata=None, resampling=RESAMPLING.cubic_spline, SOURCE_EXTRA=1000 )
		# downscale
		return utils.downscale( output_arr, baseline_arr, output_filename, downscaling_operation, \
				baseline_meta, post_downscale_function, mask, mask_value )

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_reproject_identity():
    """Reproject with an identity matrix."""
    # note the affines are both positive e, src is identity
    src = np.random.random(25).reshape((1, 5, 5))
    srcaff = Affine(1.0, 0.0, 0.0, 0.0, 1.0, 0.0)  # Identity
    srccrs = {'init': 'epsg:3857'}

    dst = np.empty(shape=(1, 10, 10))
    dstaff = Affine(0.5, 0.0, 0.0, 0.0, 0.5, 0.0)
    dstcrs = {'init': 'epsg:3857'}

    reproject(
        src, dst,
        src_transform=srcaff,
        src_crs=srccrs,
        dst_transform=dstaff,
        dst_crs=dstcrs,
        resampling=Resampling.nearest)

    # note the affines are both positive e, dst is identity
    src = np.random.random(100).reshape((1, 10, 10))
    srcaff = Affine(0.5, 0.0, 0.0, 0.0, 0.5, 0.0)
    srccrs = {'init': 'epsg:3857'}

    dst = np.empty(shape=(1, 5, 5))
    dstaff = Affine(1.0, 0.0, 0.0, 0.0, 1.0, 0.0)  # Identity
    dstcrs = {'init': 'epsg:3857'}

    reproject(
        src, dst,
        src_transform=srcaff,
        src_crs=srccrs,
        dst_transform=dstaff,
        dst_crs=dstcrs,
        resampling=Resampling.nearest)

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 to_srs_like(self, rgrid, src_nodata=None, dst_nodata=None,
                    resampling=Resampling.bilinear):
        if src_nodata is None:
            src_nodata = self.fill_value
        if dst_nodata is None:
            dst_nodata = src_nodata
        source = self.copy()
        source.fill_underlying_data(src_nodata)

        # dst_shape = rgrid.shape
        dst_transform = rgrid.gtransform
        dst_crs = rgrid.crs
        destination = rgrid.astype(self.dtype).copy()
        reproject(
            source,
            destination=destination,
            src_transform=self.gtransform,
            src_nodata=src_nodata,
            src_crs=self.crs,
            dst_transform=dst_transform,
            dst_crs=dst_crs,
            resampling=resampling,
            dst_nodata=dst_nodata)

        destination.masked_equal(dst_nodata)
        return destination

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_transform = affine.Affine.from_gdal(-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_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 run( df, meshgrid_tuple, lons_pcll, template_raster_fn, src_transform, src_crs, src_nodata, output_filename ):
	'''
	run the interpolation to a grid, and reprojection / resampling to the Alaska / Canada rasters
	extent, resolution, origin (template_raster).

	This function is intended to be used to run a pathos.multiprocessing Pool's map function
	across a list of pre-computed arguments.
			
	RETURNS:

	[str] path to the output filename generated

	'''
	template_raster = rasterio.open( template_raster_fn )
	interp_arr = xyz_to_grid( np.array(df['lon'].tolist()), \
					np.array(df['lat'].tolist()), \
					np.array(df['anom'].tolist()), grid=meshgrid_tuple, method='cubic' ) 

	src_nodata = -9999.0 # nodata
	interp_arr[ np.isnan( interp_arr ) ] = src_nodata
	dat, lons = shiftgrid( 180., interp_arr, lons_pcll, start=False )
	output_arr = np.empty_like( template_raster.read( 1 ) )
	# mask it with the internal mask in the template raster, where 0 is oob.
	output_arr = np.ma.masked_where( template_raster.read_masks( 1 ) == 0, output_arr )
	template_meta = template_raster.meta

	if 'transform' in template_meta.keys():
		template_meta.pop( 'transform' )

	reproject( dat, output_arr, src_transform=src_transform, src_crs=src_crs, src_nodata=src_nodata, \
				dst_transform=template_meta['affine'], dst_crs=template_meta['crs'],\
				dst_nodata=None, resampling=RESAMPLING.nearest, num_threads=1, SOURCE_EXTRA=1000 )	
	return write_gtiff( output_arr, template_meta, output_filename, compress=True )

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

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,
        Resampling.mode: 437298,
        Resampling.max: 439464,
        Resampling.min: 436397,
        Resampling.med: 437194,
        Resampling.q1: 436397,
        Resampling.q3: 438948,
    }

    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={"init": "epsg:3857"},
        resampling=method,
    )

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

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