Discrepancy between GTiff and MEM drivers

Question

I was expecting that both the drivers would behave the same way, except that 'MEM' will not create the file.

However, in this test 'MEM' pass, but 'GTiff' does not. Can someone explain why?

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

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

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

Answer 1

I have found the problem with the GTiff driver.

The problem is that I was not assigning a noDataValue to band2.

There is a difference between the GTiff and the MEM drivers.

When left unassigned, in the Gtiff driver assigns a nodatavalue for band2 to be equal to band1 noDataValue. However, the MEM driver does not.

This inconsistency should be fixed.

The following unitest will pass:

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

    print src_ds.ReadAsArray()
    print src_ds.GetRasterBand(2).GetNoDataValue()

    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()
    print band2
    np.testing.assert_array_equal(band2, np.array([[1., 0., 0.5],
                                                   [0.25, 0., 0.75],
                                                   [0., 0., 0.]]))


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

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