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