Target align option for Image.warp

[MWieland]

In many applications it is extremely useful to be able to force align two raster grids to match each others crs, extent, dimensions and grid alignment. With a combination of warp and mask we can come close to align two rasters (this actually works in many cases). It may, however, happen that the raster grids are shifted with subpixel resolution, which prevents them from being exactly on the same grid and therefore they cannot be compared at array level with for example numpy.

I propose the following solution to add a target_align option to warp, which aligns one Image instance to another one:

    def warp(self, dst_crs, resampling_method=0, num_threads=4, resolution=None, nodata=None, target_align=None):
        """Reproject a source raster to a destination raster.
        :param dst_crs: CRS or dict, Target coordinate reference system.
        :param resampling_method: Resampling algorithm, int, defaults to 0 (Nearest)
            numbers: https://github.com/mapbox/rasterio/blob/master/rasterio/enums.py#L28
        :param num_threads: int, number of workers, optional (default: 4)
        :param resolution: tuple (x resolution, y resolution) or float, optional.
            Target resolution, in units of target coordinate reference system.
        :param target_align: raster to which to align resolution, extent and gridspacing, optional (Image).
        :param nodata: nodata value of source, int or float, optional.
        """
        if target_align:
            transform = target_align.dataset.transform
            width = target_align.dataset.width
            height = target_align.dataset.height

        else:
            if resolution:
                transform, width, height = calculate_default_transform(
                    self.dataset.crs,
                    dst_crs,
                    self.dataset.width,
                    self.dataset.height,
                    *self.dataset.bounds,
                    resolution=resolution,
                )

            else:
                transform, width, height = calculate_default_transform(
                    self.dataset.crs, dst_crs, self.dataset.width, self.dataset.height, *self.dataset.bounds,
                )

        destination = np.zeros((self.dataset.count, height, width), self.__arr.dtype)

        self.__arr, transform = reproject(
            source=self.__arr,
            destination=destination,
            src_transform=self.dataset.transform,
            src_crs=self.dataset.crs,
            src_nodata=nodata,
            dst_transform=transform,
            dst_crs=dst_crs,
            dst_nodata=nodata,
            resampling=resampling_method,
            num_threads=num_threads,
        )

        self.dataset.close()
        self.dataset = self.__update_dataset(dst_crs, transform, nodata=nodata)