#   This Python module is part of the PyRate software package.
#
#   Copyright 2017 Geoscience Australia
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.
"""
This Python module contains tests for the gdal_python.py PyRate module.
"""
from subprocess import check_call

from pyrate.core.prepifg_helper import _resample_ifg

import glob
import os
import shutil
import subprocess
import tempfile
import unittest
from numpy import where, nan

import numpy as np
from osgeo import gdal, gdalconst

from pyrate.core import gdal_python, config as cf
from pyrate.core.shared import Ifg
from tests import common



class TestCrop(unittest.TestCase):
    def test_small_data_cropping(self):
        small_test_ifgs = common.small_data_setup()
        # minX, minY, maxX, maxY = extents
        extents = [150.91, -34.229999976, 150.949166651, -34.17]
        extents_str = [str(e) for e in extents]
        cmd = ['gdalwarp', '-overwrite', '-srcnodata', 'None', '-q', '-te'] \
              + extents_str

        for s in small_test_ifgs:
            temp_tif = tempfile.mktemp(suffix='.tif')
            t_cmd = cmd + [s.data_path,  temp_tif]
            subprocess.check_call(t_cmd)
            clipped_ref = gdal.Open(temp_tif).ReadAsArray()
            clipped = gdal_python.crop(s.data_path, extents)[0]
            np.testing.assert_array_almost_equal(clipped_ref, clipped)
            try:
                os.remove(temp_tif)
            except PermissionError:
                print("File opened by another process.")


class TestResample(unittest.TestCase):

    def test_small_data_resampling(self):
        small_test_ifgs = common.small_data_setup()
        # minX, minY, maxX, maxY = extents
        extents = [150.91, -34.229999976, 150.949166651, -34.17]
        extents_str = [str(e) for e in extents]
        resolutions = [0.001666666, .001, 0.002, 0.0025, .01]

        for res in resolutions:
            res = [res, -res]
            self.check_same_resampled_output(extents, extents_str, res,
                                             small_test_ifgs)

    def check_same_resampled_output(self, extents, extents_str, res,
                                    small_test_ifgs):
        cmd = ['gdalwarp', '-overwrite', '-srcnodata', 'None',
               '-q', '-r', 'near', '-te'] \
              + extents_str

        if res[0]:
            new_res_str = [str(r) for r in res]
            cmd += ['-tr'] + new_res_str
        for s in small_test_ifgs:
            temp_tif = tempfile.mktemp(suffix='.tif')
            t_cmd = cmd + [s.data_path, temp_tif]
            subprocess.check_call(t_cmd)
            resampled_ds = gdal.Open(temp_tif)
            resampled_ref = resampled_ds.ReadAsArray()

            resampled_temp_tif = tempfile.mktemp(suffix='.tif',
                                                 prefix='resampled_')
            resampled = gdal_python.resample_nearest_neighbour(s.data_path,
                                                               extents, res,
                                                               resampled_temp_tif)
            np.testing.assert_array_almost_equal(resampled_ref, resampled[0, :, :])
            try:
                os.remove(temp_tif)
            except PermissionError:
                print("File opened by another process.")

            try:
                os.remove(resampled_temp_tif)  # also proves file was written
            except PermissionError:
                print("File opened by another process.")

    def test_none_resolution_output(self):
        small_test_ifgs = common.small_data_setup()
        # minX, minY, maxX, maxY = extents
        extents = [150.91, -34.229999976, 150.949166651, -34.17]
        extents_str = [str(e) for e in extents]

        self.check_same_resampled_output(extents, extents_str, [None, None],
                                         small_test_ifgs)

    def test_output_file_written(self):
        small_test_ifgs = common.small_data_setup()
        extents = [150.91, -34.229999976, 150.949166651, -34.17]
        resolutions = [0.001666666, .001, 0.002, 0.0025, .01]
        for res in resolutions:
            for s in small_test_ifgs:
                resampled_temp_tif = tempfile.mktemp(suffix='.tif',
                                                    prefix='resampled_')
                gdal_python.resample_nearest_neighbour(s.data_path, extents,
                                                       [res, -res],
                                                       resampled_temp_tif)
                self.assertTrue(os.path.exists(resampled_temp_tif))
                os.remove(resampled_temp_tif)

    def test_small_data_crop_vs_resample(self):
        small_test_ifgs = common.small_data_setup()
        # minX, minY, maxX, maxY = extents
        extents = [150.91, -34.229999976, 150.949166651, -34.17]
        for s in small_test_ifgs:
            clipped = gdal_python.crop(s.data_path, extents)[0]
            resampled_temp_tif = tempfile.mktemp(suffix='.tif',
                                                prefix='resampled_')
            resampled = gdal_python.resample_nearest_neighbour(
                s.data_path, extents, [None, None], resampled_temp_tif)
            self.assertTrue(os.path.exists(resampled_temp_tif))
            np.testing.assert_array_almost_equal(resampled[0, :, :], clipped)
            os.remove(resampled_temp_tif)

    def test_resampled_tif_has_metadata(self):
        small_test_ifgs = common.small_data_setup()

        # minX, minY, maxX, maxY = extents
        extents = [150.91, -34.229999976, 150.949166651, -34.17]
        for s in small_test_ifgs:

            resampled_temp_tif = tempfile.mktemp(suffix='.tif',
                                                prefix='resampled_')
            gdal_python.resample_nearest_neighbour(
                s.data_path, extents, [None, None], resampled_temp_tif)
            dst_ds = gdal.Open(resampled_temp_tif)
            md = dst_ds.GetMetadata()
            self.assertDictEqual(md, s.meta_data)
            try:
                os.remove(resampled_temp_tif)
            except PermissionError:
                print("File opened by another process.")


class BasicReampleTests(unittest.TestCase):

    def test_reproject_with_no_data(self):

        data = np.array([[2, 7],
                         [2, 7]])
        src_ds = gdal.GetDriverByName('MEM').Create('', 2, 2)
        src_ds.GetRasterBand(1).WriteArray(data)
        src_ds.GetRasterBand(1).SetNoDataValue(2)
        src_ds.SetGeoTransform([10, 1, 0, 10, 0, -1])

        dst_ds = gdal.GetDriverByName('MEM').Create('', 1, 1)
        dst_ds.GetRasterBand(1).SetNoDataValue(3)
        dst_ds.GetRasterBand(1).Fill(3)
        dst_ds.SetGeoTransform([10, 2, 0, 10, 0, -2])

        gdal.ReprojectImage(src_ds, dst_ds, '', '', gdal.GRA_NearestNeighbour)
        got_data = dst_ds.GetRasterBand(1).ReadAsArray()
        expected_data = np.array([[7]])
        np.testing.assert_array_equal(got_data, expected_data)

    def test_reproject_with_no_data_2(self):

        data = np.array([[2, 7, 7, 7],
                         [2, 7, 7, 2]])
        height, width = data.shape
        src_ds = gdal.GetDriverByName('MEM').Create('', width, height)
        src_ds.GetRasterBand(1).WriteArray(data)
        src_ds.GetRasterBand(1).SetNoDataValue(2)
        src_ds.SetGeoTransform([10, 1, 0, 10, 0, -1])

        dst_ds = gdal.GetDriverByName('MEM').Create('', 2, 1)
        dst_ds.GetRasterBand(1).SetNoDataValue(3)
        dst_ds.GetRasterBand(1).Fill(3)
        dst_ds.SetGeoTransform([10, 2, 0, 10, 0, -2])

        gdal.ReprojectImage(src_ds, dst_ds, '', '', gdal.GRA_NearestNeighbour)
        got_data = dst_ds.GetRasterBand(1).ReadAsArray()
        expected_data = np.array([[7, 3]])
        np.testing.assert_array_equal(got_data, expected_data)

    def test_reproject_with_no_data_3(self):

        data = np.array([[2, 7, 7, 7],
                         [2, 7, 7, 7],
                         [2, 7, 7, 7],
                         [2, 7, 7, 2],
                         [2, 7, 7, 2]])
        src_ds = gdal.GetDriverByName('MEM').Create('', 4, 5)
        src_ds.GetRasterBand(1).WriteArray(data)
        src_ds.GetRasterBand(1).SetNoDataValue(2)
        src_ds.SetGeoTransform([10, 1, 0, 10, 0, -1])

        dst_ds = gdal.GetDriverByName('MEM').Create('', 2, 2)
        dst_ds.GetRasterBand(1).SetNoDataValue(3)
        dst_ds.GetRasterBand(1).Fill(3)
        dst_ds.SetGeoTransform([10, 2, 0, 10, 0, -2])

        gdal.ReprojectImage(src_ds, dst_ds, '', '', gdal.GRA_NearestNeighbour)
        got_data = dst_ds.GetRasterBand(1).ReadAsArray()
        expected_data = np.array([[7, 7],
                                  [7, 3]])
        np.testing.assert_array_equal(got_data, expected_data)

    def test_reproject_with_no_data_4(self):

        data = np.array([[2, 7, 7, 7, 2],
                         [2, 7, 7, 7, 2],
                         [2, 7, 7, 7, 2],
                         [2, 7, 7, 2, 2],
                         [2, 7, 7, 2, 2]])
        src_ds = gdal.GetDriverByName('MEM').Create('', 5, 5)
        src_ds.GetRasterBand(1).WriteArray(data)
        src_ds.GetRasterBand(1).SetNoDataValue(2)
        src_ds.SetGeoTransform([10, 1, 0, 10, 0, -1])

        dst_ds = gdal.GetDriverByName('MEM').Create('', 2, 2)
        dst_ds.GetRasterBand(1).SetNoDataValue(3)
        dst_ds.GetRasterBand(1).Fill(3)
        dst_ds.SetGeoTransform([10, 2, 0, 10, 0, -2])

        gdal.ReprojectImage(src_ds, dst_ds, '', '', gdal.GRA_NearestNeighbour)
        got_data = dst_ds.GetRasterBand(1).ReadAsArray()
        expected_data = np.array([[7, 7],
                                  [7, 3]])
        np.testing.assert_array_equal(got_data, expected_data)


    def test_reproject_with_no_data_5(self):

        data = np.array([[2, 7, 7, 7, 2],
                         [2, 7, 7, 7, 2],
                         [2, 7, 7, 7, 2],
                         [2, 7, 7, 2, 2],
                         [2, 7, 7, 2, 2],
                         [2, 7, 7, 2, 2]])
        src_ds = gdal.GetDriverByName('MEM').Create('', 5, 6)
        src_ds.GetRasterBand(1).WriteArray(data)
        src_ds.GetRasterBand(1).SetNoDataValue(2)
        src_ds.SetGeoTransform([10, 1, 0, 10, 0, -1])

        dst_ds = gdal.GetDriverByName('MEM').Create('', 2, 3)
        dst_ds.GetRasterBand(1).SetNoDataValue(3)
        dst_ds.GetRasterBand(1).Fill(3)
        dst_ds.SetGeoTransform([10, 2, 0, 10, 0, -2])

        gdal.ReprojectImage(src_ds, dst_ds, '', '', gdal.GRA_NearestNeighbour)
        got_data = dst_ds.GetRasterBand(1).ReadAsArray()
        expected_data = np.array([[7, 7],
                                  [7, 3],
                                  [7, 3]])
        np.testing.assert_array_equal(got_data, expected_data)


    def test_reproject_average_resampling(self):

        data = np.array([[4, 7, 7, 7, 2, 7.],
                         [4, 7, 7, 7, 2, 7.],
                         [4, 7, 7, 7, 2, 7.],
                         [4, 7, 7, 2, 2, 7.],
                         [4, 7, 7, 2, 2, 7.],
                         [4, 7, 7, 10, 2, 7.]], dtype=np.float32)
        src_ds = gdal.GetDriverByName('MEM').Create('', 6, 6, 1,
                                                    gdalconst.GDT_Float32)
        src_ds.GetRasterBand(1).WriteArray(data)
        src_ds.GetRasterBand(1).SetNoDataValue(2)
        src_ds.SetGeoTransform([10, 1, 0, 10, 0, -1])

        dst_ds = gdal.GetDriverByName('MEM').Create('', 3, 3, 1,
                                                    gdalconst.GDT_Float32)
        dst_ds.GetRasterBand(1).SetNoDataValue(3)
        dst_ds.GetRasterBand(1).Fill(3)
        dst_ds.SetGeoTransform([10, 2, 0, 10, 0, -2])

        gdal.ReprojectImage(src_ds, dst_ds, '', '', gdal.GRA_Average)
        got_data = dst_ds.GetRasterBand(1).ReadAsArray()
        expected_data = np.array([[5.5, 7, 7],
                                  [5.5, 7, 7],
                                  [5.5, 8, 7]])
        np.testing.assert_array_equal(got_data, expected_data)

    def test_reproject_average_resampling_with_2bands(self):

        data = np.array([[[4, 7, 7, 7, 2, 7.],
                         [4, 7, 7, 7, 2, 7.],
                         [4, 7, 7, 7, 2, 7.],
                         [4, 7, 7, 2, 2, 7.],
                         [4, 7, 7, 2, 2, 7.],
                         [4, 7, 7, 10, 2, 7.]],
                        [[2, 0, 0, 0, 0, 0.],
                         [2, 0, 0, 0, 0, 2.],
                         [0, 1., 0, 0, 0, 1.],
                         [0, 0, 0, 0, 0, 2],
                         [0, 0, 0, 0, 0, 0.],
                         [0, 0, 0, 0, 0, 0.]]], dtype=np.float32)
        src_ds = gdal.GetDriverByName('MEM').Create('', 6, 6, 2,
                                                    gdalconst.GDT_Float32)

        src_ds.GetRasterBand(1).WriteArray(data[0, :, :])
        src_ds.GetRasterBand(1).SetNoDataValue(2)
        src_ds.GetRasterBand(2).WriteArray(data[1, :, :])
        # src_ds.GetRasterBand(1).SetNoDataValue()
        src_ds.SetGeoTransform([10, 1, 0, 10, 0, -1])

        dst_ds = gdal.GetDriverByName('MEM').Create('', 3, 3, 2,
                                                    gdalconst.GDT_Float32)
        dst_ds.GetRasterBand(1).SetNoDataValue(3)
        dst_ds.GetRasterBand(1).Fill(3)
        dst_ds.SetGeoTransform([10, 2, 0, 10, 0, -2])

        gdal.ReprojectImage(src_ds, dst_ds, '', '', gdal.GRA_Average)
        got_data = dst_ds.GetRasterBand(1).ReadAsArray()
        expected_data = np.array([[5.5, 7, 7],
                                  [5.5, 7, 7],
                                  [5.5, 8, 7]])
        np.testing.assert_array_equal(got_data, expected_data)
        band2 = dst_ds.GetRasterBand(2).ReadAsArray()
        np.testing.assert_array_equal(band2, np.array([[1., 0., 0.5],
                                                       [0.25, 0., 0.75],
                                                       [0., 0., 0.]]))


class TestOldPrepifgVsGdalPython(unittest.TestCase):

    def setUp(self):
        self.test_dir = tempfile.mktemp()
        common.copytree(common.SML_TEST_TIF, self.test_dir)
        self.ifgs = common.small_data_setup(datafiles=glob.glob(os.path.join(self.test_dir, "*.tif")))
        self.ref_gtif = gdal.Open(self.ifgs[0].data_path, gdalconst.GA_ReadOnly)
        self.ref_proj = self.ref_gtif.GetProjection()
        self.ref_gt = self.ref_gtif.GetGeoTransform()
        self.data = self.ref_gtif.ReadAsArray()
        self.md = self.ref_gtif.GetMetadata()
        self.temp_tif = tempfile.mktemp(suffix='.tif')
        self.out_tif = tempfile.mktemp(suffix='.tif')
        del self.ref_gtif

    def tearDown(self):
        if os.path.exists(self.temp_tif):
            try:
                os.remove(self.temp_tif)
            except PermissionError:
                print("File opened by another process.")
        try:
            shutil.rmtree(self.test_dir)
        except PermissionError:
            print("File opened by another process.")

    def test_gdal_python_vs_old_prepifg_prep(self):

        for i in range(10):
            data = np.array(np.random.randint(0, 3, size=(10, 10)), dtype=np.float32)

            src_ds = gdal.GetDriverByName('GTiff').Create(self.temp_tif, 10, 10, 2, gdalconst.GDT_Float32)

            src_ds.GetRasterBand(1).WriteArray(data)
            src_ds.GetRasterBand(1).SetNoDataValue(0)
            nan_matrix = where(data == 0, nan, data)

            src_ds.GetRasterBand(2).WriteArray(np.isnan(nan_matrix))
            src_ds.GetRasterBand(2).SetNoDataValue(0)
            src_ds.SetGeoTransform([10, 1, 0, 10, 0, -1])

            dst_ds = gdal.GetDriverByName('MEM').Create('', 5, 5, 2, gdalconst.GDT_Float32)
            dst_ds.SetGeoTransform([10, 2, 0, 10, 0, -2])

            for k, v in self.md.items():
                src_ds.SetMetadataItem(k, v)
                dst_ds.SetMetadataItem(k, v)

            src_ds.FlushCache()
            gdal.ReprojectImage(src_ds, dst_ds, '', '', gdal.GRA_Average)
            nan_frac = dst_ds.GetRasterBand(2).ReadAsArray()
            avg = dst_ds.GetRasterBand(1).ReadAsArray()
            thresh = 0.1
            avg[nan_frac >= thresh] = np.nan

            ifg = Ifg(self.temp_tif)
            x_looks = y_looks = 2
            res = 2
            resolution = [res, -res]
            # minX, minY, maxX, maxY = extents
            # adfGeoTransform[0] /* top left x */
            # adfGeoTransform[1] /* w-e pixel resolution */
            # adfGeoTransform[2] /* 0 */
            # adfGeoTransform[3] /* top left y */
            # adfGeoTransform[4] /* 0 */
            # adfGeoTransform[5] /* n-s pixel resolution (negative value) */
            extents = [str(e) for e in [10, 0, 20, 10]]

            # only band 1 is resapled in warp_old
            data, self.old_prepifg_path = warp_old(ifg, x_looks, y_looks, extents, resolution,
                                                   thresh=thresh, crop_out=4, verbose=False)

            np.testing.assert_array_equal(data, avg)

            del ifg
            del src_ds

    def test_gdal_python_vs_old_prepifg_prep2(self):

        for i in range(10):
            thresh = 0.5
            x_looks = y_looks = 2
            res = 2
            extents = [10, 0, 20, 10]
            extents_str = [str(e) for e in extents]
            data = np.array(np.random.randint(0, 3, size=(10, 10)), dtype=np.float32)

            # create the  self.temp_tiff
            self.manipulation(data, self.temp_tif, self.md)

            # first file already opened by some process and never closed
            self.out_tif = self.out_tif.split(".")[0]+"_2.tif"

            # only band 1 is resapled in warp_old
            averaged_and_resampled, _ = gdal_python.crop_resample_average(self.temp_tif, extents, [res, -res], self.out_tif, thresh, match_pyrate=True)
            ifg = Ifg(self.temp_tif)
            # only band 1 is resampled in warp_old
            data, self.old_prepifg_path = warp_old(ifg, x_looks, y_looks, extents_str, [res, -res],
                                                   thresh=thresh, crop_out=4, verbose=False)

            np.testing.assert_array_almost_equal(data, averaged_and_resampled, decimal=4)
            del averaged_and_resampled

    @staticmethod
    def manipulation(data, tiff, md):
        src_ds = gdal.GetDriverByName('GTiff').Create(tiff, 10, 10, 2, gdalconst.GDT_Float32)
        src_ds.GetRasterBand(1).WriteArray(data)
        src_ds.GetRasterBand(1).SetNoDataValue(0)
        nan_matrix = where(data == 0, nan, data)
        src_ds.GetRasterBand(2).WriteArray(np.isnan(nan_matrix))
        src_ds.GetRasterBand(2).SetNoDataValue(-100)
        src_ds.SetGeoTransform([10, 1, 0, 10, 0, -1])
        dst_ds = gdal.GetDriverByName('MEM').Create('', 5, 5, 2, gdalconst.GDT_Float32)
        dst_ds.SetGeoTransform([10, 2, 0, 10, 0, -2])
        for k, v in md.items():
            src_ds.SetMetadataItem(k, v)
            dst_ds.SetMetadataItem(k, v)
        src_ds.FlushCache()

        del dst_ds
        del src_ds

    def test_gdal_python_vs_old_prepifg_no_match_pyrate(self):

        for ifg in self.ifgs:
            extents = [150.91, -34.229999976, 150.949166651, -34.17]
            extents_str = [str(e) for e in extents]
            orig_res = 0.000833333
            thresh = 0.5

            for looks in range(10):
                x_looks = y_looks = looks
                res = orig_res*x_looks
                averaged_and_resapled, out_ds = gdal_python.crop_resample_average(
                    ifg.data_path, extents, new_res=[res, -res],
                    output_file=self.temp_tif, thresh=thresh,
                    match_pyrate=False)

                # only band 1 is resampled in warp_old
                data, self.old_prepifg_path = warp_old(
                    ifg, x_looks, y_looks, extents_str, [res, -res],
                    thresh=thresh, crop_out=4, verbose=False)
                yres, xres = data.shape

                # old_prepifg warp resample method loses one
                # row at the bottom if nrows % 2 == 1
                averaged_and_resapled = averaged_and_resapled[:yres, :xres]
                nrows, ncols = averaged_and_resapled.shape
                np.testing.assert_array_almost_equal(data, averaged_and_resapled, decimal=4)

                # make sure they are the same after they are opened again
                # Last [yres:nrows, xres:ncols] won't match due to PyRate
                # dropping last few rows/columns depending on resolution/looks
                data_from_file = gdal.Open(self.old_prepifg_path).ReadAsArray()
                new_from_file = out_ds.ReadAsArray()
                np.testing.assert_array_almost_equal(data_from_file[:yres, :xres], new_from_file[:yres, :xres])

                del out_ds  # manual close
                del data_from_file
                del averaged_and_resapled
                del data


    def test_gdal_python_vs_old_prepifg(self):

        self.ifgs = common.small_data_setup(
            datafiles=glob.glob(os.path.join(self.test_dir, "*.tif")))

        for ifg in self.ifgs:
            extents = [150.91, -34.229999976, 150.949166651, -34.17]
            extents_str = [str(e) for e in extents]
            orig_res = 0.000833333
            thresh = 0.5
            for looks in range(1, 10):
                x_looks = y_looks = looks
                res = orig_res*x_looks
                averaged_and_resampled = gdal_python.crop_resample_average(
                    ifg.data_path, extents, new_res=[res, -res],
                    output_file=self.temp_tif, thresh=thresh,
                    match_pyrate=True)[0]

                # only band 1 is resampled in warp_old
                data, self.old_prepifg_path = warp_old(
                    ifg, x_looks, y_looks, extents_str, [res, -res],
                    thresh=thresh, crop_out=4, verbose=False)
                yres, xres = data.shape

                # old_prepifg warp resample method loses one row
                # at the bottom if nrows % 2 == 1
                np.testing.assert_array_almost_equal(data, averaged_and_resampled[:yres, :xres], decimal=4)
                nrows, ncols = averaged_and_resampled.shape
                # make sure they are the same after they are opened again
                data_from_file = gdal.Open(self.old_prepifg_path).ReadAsArray()
                self.assertTrue(os.path.exists(self.temp_tif))
                new_from_file = gdal.Open(self.temp_tif).ReadAsArray()
                np.testing.assert_array_almost_equal(data_from_file, new_from_file)

                del averaged_and_resampled
                del new_from_file
                del data_from_file

    def test_no_out_file_when_driver_type_is_mem(self):
        for ifg in self.ifgs:
            extents = [150.91, -34.229999976, 150.949166651, -34.17]
            extents_str = [str(e) for e in extents]
            orig_res = 0.000833333
            thresh = 0.5
            x_looks = y_looks = 6
            res = orig_res*x_looks
            averaged_and_resampled, out_ds = gdal_python.crop_resample_average(
                ifg.data_path, extents, new_res=[res, -res],
                output_file=self.temp_tif, thresh=thresh,
                out_driver_type='MEM', match_pyrate=True)

            # only band 1 is resampled in warp_old
            data, self.old_prepifg_path = warp_old(
                ifg, x_looks, y_looks, extents_str, [res, -res],
                thresh=thresh, crop_out=4, verbose=False)
            yres, xres = data.shape

            # old_prepifg warp resample method loses one row at the bottom if
            # nrows % 2 == 1
            averaged_and_resampled = averaged_and_resampled[:yres, :xres]
            np.testing.assert_array_almost_equal(data, averaged_and_resampled,
                                                 decimal=4)

            # make sure they are the same after they are opened again
            data_from_file = gdal.Open(self.old_prepifg_path).ReadAsArray()
            new_from_file = out_ds.ReadAsArray()
            np.testing.assert_array_almost_equal(data_from_file,
                                                 new_from_file)
            self.assertFalse(os.path.exists(self.temp_tif))
            out_ds = None  # manual close


class TestMEMVsGTiff(unittest.TestCase):

    @staticmethod
    def check(driver_type):

        temp_tif = tempfile.mktemp(suffix='.tif')

        data = np.array([[[4, 7, 7, 7, 2, 7.],
                         [4, 7, 7, 7, 2, 7.],
                         [4, 7, 7, 7, 2, 7.],
                         [4, 7, 7, 2, 2, 7.],
                         [4, 7, 7, 2, 2, 7.],
                         [4, 7, 7, 10, 2, 7.]],
                        [[2, 0, 0, 0, 0, 0.],
                         [2, 0, 0, 0, 0, 2.],
                         [0, 1., 0, 0, 0, 1.],
                         [0, 0, 0, 0, 0, 2],
                         [0, 0, 0, 0, 0, 0.],
                         [0, 0, 0, 0, 0, 0.]]], dtype=np.float32)
        src_ds = gdal.GetDriverByName(driver_type).Create(temp_tif, 6, 6, 2,
                                                    gdalconst.GDT_Float32)

        src_ds.GetRasterBand(1).WriteArray(data[0, :, :])
        src_ds.GetRasterBand(1).SetNoDataValue(2)
        src_ds.GetRasterBand(2).WriteArray(data[1, :, :])
        src_ds.GetRasterBand(2).SetNoDataValue(3)
        src_ds.SetGeoTransform([10, 1, 0, 10, 0, -1])
        src_ds.FlushCache()

        dst_ds = gdal.GetDriverByName('MEM').Create('', 3, 3, 2,
                                                    gdalconst.GDT_Float32)
        dst_ds.GetRasterBand(1).SetNoDataValue(3)
        dst_ds.GetRasterBand(1).Fill(3)
        dst_ds.SetGeoTransform([10, 2, 0, 10, 0, -2])

        gdal.ReprojectImage(src_ds, dst_ds, '', '', gdal.GRA_Average)
        band1 = dst_ds.GetRasterBand(1).ReadAsArray()
        np.testing.assert_array_equal(band1, np.array([[5.5, 7, 7],
                                                       [5.5, 7, 7],
                                                       [5.5, 8, 7]]))
        band2 = dst_ds.GetRasterBand(2).ReadAsArray()
        np.testing.assert_array_equal(band2, np.array([[1., 0., 0.5],
                                                       [0.25, 0., 0.75],
                                                       [0., 0., 0.]]))
        if os.path.exists(temp_tif):
            try:
                os.remove(temp_tif)
            except PermissionError:
                print("File opened by another process.")

    def test_mem(self):
        self.check('MEM')

    def test_gtiff(self):
        self.check('GTiff')


# class TestGDalAverageResampleing(unittest.TestCase):
#     pass


def warp_old(ifg, x_looks, y_looks, extents, resolution, thresh,
             crop_out, verbose, ret_ifg=True):
    """
    Resamples 'ifg' and returns a new Ifg obj.

    :param xlooks: integer factor to scale X axis by, 5 is 5x smaller,
        1 is no change.
    :param ylooks: as xlooks, but for Y axis
    :param extents: georeferenced extents for new file: (xfirst, yfirst, xlast, ylast)
    :param resolution: [xres, yres] or None. Sets resolution output Ifg metadata.
         Use *None* if raster size is not being changed.
    :param thresh: see thresh in prepare_ifgs().
    :param verbose: True to print gdalwarp output to stdout
    """
    # pylint: disable=too-many-locals
    if x_looks != y_looks:
        raise ValueError('X and Y looks mismatch')

    # dynamically build command for call to gdalwarp
    cmd = ["gdalwarp", "-overwrite", "-srcnodata", "None", "-te"] + extents
    if not verbose:
        cmd.append("-q")

    # It appears that before vrions 1.10 gdal-warp did not copy meta-data
    # ... so we need to. Get the keys from the input here
    if (sum((int(v)*i for v, i
             in zip(gdal.__version__.split('.')[:2], [10, 1]))) < 20):
        fl = ifg.data_path
        dat = gdal.Open(fl)
        md = {k:v for k, v in dat.GetMetadata().items()}
    else:
        md = None

    # HACK: if resampling, cut segment with gdalwarp & manually average tiles
    data = None
    if resolution[0]:
        data = _resample_ifg(ifg, cmd, x_looks, y_looks, thresh, md)
        cmd += ["-tr"] + [str(r) for r in resolution] # change res of final output

    # use GDAL to cut (and resample) the final output layers
    looks_path = cf.mlooked_path(ifg.data_path, y_looks, crop_out)
    cmd += [ifg.data_path, looks_path]

    check_call(cmd)
    # now write the metadata from the input to the output
    if md is not None:
        new_lyr = gdal.Open(looks_path)
        for k, v in md.iteritems():
            new_lyr.SetMetadataItem(k, v)

    # Add missing/updated metadata to resampled ifg/DEM
    new_lyr = type(ifg)(looks_path)
    new_lyr.open(readonly=False)
    # for non-DEMs, phase bands need extra metadata & conversions
    if hasattr(new_lyr, "phase_band"):
        if data is None:  # data wasn't resampled, so flag incoherent cells
            data = new_lyr.phase_band.ReadAsArray()
            data = np.where(np.isclose(data, 0.0, atol=1e-6), np.nan, data)

        # TODO: LOS conversion to vertical/horizontal (projection)
        # TODO: push out to workflow
        #if params.has_key(REPROJECTION_FLAG):
        #    reproject()

        # tricky: write either resampled or the basic cropped data to new layer
        new_lyr.phase_band.SetNoDataValue(nan)
        new_lyr.phase_band.WriteArray(data)
        new_lyr.nan_converted = True

    if ret_ifg:
        return data, looks_path
    else:
        return

if __name__ == '__main__':
    unittest.main()
