'''
This module contains the unit tests. Run "python tests.py" at the command line
to run the tests.
'''
import hashlib
import os
import random
import unittest
import unmixing
from unmixing.utils import *
from unmixing.lsma import FCLSAbundanceMapper, PPI, convex_hull_graham, endmembers_by_maximum_angle, endmembers_by_maximum_area, endmembers_by_query, endmembers_by_maximum_volume, hall_rectification, point_to_pixel_geometry
from unmixing.sasma import concat_endmember_arrays
from unmixing.transform import biophysical_composition_index, tasseled_cap_tm, mnf_rotation
from unmixing.visualize import FeatureSpace
from osgeo import gdal
from pysptools.noise import MNF
# For backwards compatibility in GDAL
gdal.SetConfigOption('GDAL_ARRAY_OPEN_BY_FILENAME', 'TRUE')
TEST_DIR = os.path.join(os.path.dirname(unmixing.__file__), 'test')
# Required to get consistent results for WKT strings
os.environ['OGR_WKT_PRECISION'] = '7'
class Tester(unittest.TestCase):
@classmethod
def tearDownClass(self):
for p in [
'temp.tiff', 'temp2.tiff', 'temp3.tiff', 'FeatureSpace_selection_test_1.kml',
'rect_multi7_raster2.tiff'
]:
try:
os.unlink(os.path.join(self.test_dir, p))
except FileNotFoundError:
pass
class FCLS(Tester):
test_dir = TEST_DIR
test_data_100_110_hex_string = 'bea34345168cd1858fe5f9066c4403a806b9d6f2ffa7ddbdb3dd0ab354b513cd'
def test_fcls_unmixing_with_single_endmember_spectra(self):
'''
Should calculate abundance from a mixed image for single endmember
spectra.
'''
em_locs = [(326701, 4696895),(324978, 4699651), (328823, 4696835)]
arr, gt, wkt = as_array(os.path.join(self.test_dir,
'LT05_020030_merge_19950712_stack_clip.tiff'))
arr_mnf = mnf_rotation(arr).T
endmembers = spectra_at_xy(arr_mnf, em_locs, gt, wkt)
fcls_mapper = FCLSAbundanceMapper(arr_mnf[:,100:110,100:110],
gt, wkt, processes = 1)
result = fcls_mapper.map_abundance(endmembers)
hasher = hashlib.sha256()
hasher.update(result)
self.assertEqual(hasher.hexdigest(), self.test_data_100_110_hex_string)
def test_fcls_unmixing_with_single_endmember_spectra_multicore(self):
'''
Should calculate abundance from a mixed image for single endmember
spectra; result should be the same for 1 or 2 processes.
'''
em_locs = [(326701, 4696895),(324978, 4699651), (328823, 4696835)]
arr, gt, wkt = as_array(os.path.join(self.test_dir,
'LT05_020030_merge_19950712_stack_clip.tiff'))
endmembers = spectra_at_xy(mnf_rotation(arr).T, em_locs, gt, wkt)
arr_mnf = mnf_rotation(arr).T
endmembers = spectra_at_xy(arr_mnf, em_locs, gt, wkt)
fcls_mapper1 = FCLSAbundanceMapper(arr_mnf[:,100:110,100:110],
gt, wkt, processes = 1)
fcls_mapper2 = FCLSAbundanceMapper(arr_mnf[:,100:110,100:110],
gt, wkt, processes = 2)
result1 = fcls_mapper1.map_abundance(endmembers)
result2 = fcls_mapper2.map_abundance(endmembers)
self.assertTrue(np.all(np.equal(result1, result2)))
hasher = hashlib.sha256()
hasher.update(result2)
self.assertEqual(hasher.hexdigest(), self.test_data_100_110_hex_string)
class SASMA(Tester):
test_dir = TEST_DIR
def test_concatenation_of_endmember_arrays(self):
'''
Spectra arrays for multiple endmember types should be concatenated
correctly; this step anticipates LSMA with multiple endmember spectra.
'''
arr, gt, wkt = as_array(os.path.join(self.test_dir,
'LT05_020030_merge_19950712_stack_clip.tiff'))
vbd, gt, wkt = as_array(os.path.join(self.test_dir,
'LT05_020030_merge_19950712_VBD_endmember_PIFs.tiff'))
emv = np.where(vbd == 1, arr, 0)
emb = np.where(vbd == 2, arr, 0)
emd = np.where(vbd == 3, arr, 0)
endmembers = concat_endmember_arrays(emv, emb, emd)
self.assertTrue(np.all(np.equal(endmembers.shape, (56639, 3, 6))))
class LSMA(Tester):
test_dir = TEST_DIR
cases = {
'Vegetation': [
(341663, 4709229), (314333, 4694229), (301583, 4697919),
(320843, 4709769), (288053, 4741239)
],
'High/Bright': [
(331583, 4690839), (343313, 4707999), (351143, 4718739),
(334913, 4709109), (305603, 4676919), (300683, 4707459),
(318923, 4724259)
],
'Low/Dark': [
(325973, 4682799), (322193, 4678389), (321383, 4686279),
(320033, 4694229), (334793, 4700829), (350393, 4727439)
]
}
def test_endmember_search_by_maximum_area(self):
'''Should be able to search for endmembers by maximum area.'''
rast, gt, wkt = as_array(os.path.join(self.test_dir, 'multi7_raster.tiff'))
hsi_post_mnf = mnf_rotation(rast)
# Test that dictionary type works
result = endmembers_by_maximum_area(hsi_post_mnf.T, self.cases,
gt=gt, wkt=wkt, dd=False)
self.assertTrue(isinstance(result[0], np.ndarray))
self.assertTrue(list(map(np.ceil, result[0][0].tolist())),
[15.0, -9.0, 1.0])
self.assertTrue(list(map(np.ceil, result[0][0].tolist())),
[16.0, -8.0, 2.0])
# Test that sequence/ array type works
result = endmembers_by_maximum_area(hsi_post_mnf.T, self.cases['Vegetation'],
gt=gt, wkt=wkt, dd=False)
self.assertTrue(isinstance(result[0], np.ndarray))
self.assertTrue(list(map(np.ceil, result[0][0].tolist())),
[-12.0, 0.0, 0.0])
self.assertTrue(list(map(np.ceil, result[0][0].tolist())),
[-11.0, 1.0, 1.0])
def test_endmember_search_by_maximum_volume(self):
'''Should be able to search for endmembers by maximum volume.'''
rast, gt, wkt = as_array(os.path.join(self.test_dir, 'multi7_raster.tiff'))
hsi_post_mnf = mnf_rotation(rast)
# Test that dictionary type works
result = endmembers_by_maximum_volume(hsi_post_mnf.T, self.cases,
gt=gt, wkt=wkt, dd=False)
self.assertTrue(isinstance(result[0], np.ndarray))
# Test that sequence/ array type works
result = endmembers_by_maximum_volume(hsi_post_mnf.T, self.cases['Vegetation'],
gt=gt, wkt=wkt, dd=False)
self.assertTrue(isinstance(result[0], np.ndarray))
ref_target = (300953, 4682679)
result = endmembers_by_maximum_volume(hsi_post_mnf.T, self.cases['Vegetation'],
ref_target, gt=gt, wkt=wkt, dd=False)
self.assertTrue(isinstance(result[0], np.ndarray))
def test_endmember_search_by_maximum_angle(self):
'''Should be able to search for endmembers by maximum solid angle.'''
rast, gt, wkt = as_array(os.path.join(self.test_dir, 'multi7_raster.tiff'))
hsi = rast.transpose()
hsi[hsi==-9999] = 0
mnf = MNF()
hsi_post_mnf = mnf.apply(hsi)
# Test that sequence/ array type works
ref_target = (342683, 4703919)
result = endmembers_by_maximum_angle(hsi_post_mnf.T, self.cases['Vegetation'],
ref_target, gt=gt, wkt=wkt, dd=False)
self.assertTrue(isinstance(result[0], np.ndarray))
def test_composite(self):
'''Compositing multi-date images should produce the expected result.'''
nodata = -9999
ex1 = [
np.array([
[[-9999., 0.1], [-9999., -9999.]],
[[-9999., 0.5], [-9999., -9999.]]
]),
np.array([
[[-9999., 0.9], [-9999., -9999.]],
[[-9999., 0.2], [-9999., -9999.]]
]),
np.array([
[[-9999., 0.03], [-9999., -9999.]],
[[-9999., 0.03], [-9999., -9999.]]
]),
np.array([
[[-9999., -9999.], [-9999., 0.2]],
[[-9999., -9999.], [-9999., 0.2]]
]),
np.array([
[[-9999., 0.3], [0.4, 0.5]],
[[-9999., 0.3], [0.4, 0.5]]
])
]
# These tests use the function with masks that handles NoData correctly
self.assertTrue(np.where(np.round(composite2(('max', 'median'), *ex1, normalize=None), 2) == np.array([[[-9999., 0.9], [0.4, 0.5]], [[-9999., 0.25], [0.4, 0.35]]]), True, False).all())
self.assertTrue(np.where(np.round(composite2(('min', 'median'), *ex1, normalize=None), 2) == np.array([[[-9999., 0.03], [0.4, 0.2]], [[-9999., 0.25], [0.4, 0.35]]]), True, False).all())
self.assertTrue(np.where(np.round(composite2(('max', 'median'), *ex1, normalize=None, dtype=int), 2) == np.array([[[-9999., 0.9], [0.4, 0.5]], [[-9999., 0.25], [0.4, 0.35]]]), True, False).all()) # Setting `dtype=int` simulates the type coercion error we want to test
# These tests use the function without masks that requires star functions
self.assertTrue(np.where(np.round(composite3(('max', 'median'), *ex1, normalize=None), 2) == np.array([[[-9999., 0.9], [0.4, 0.5]], [[-9999., 0.25], [0.4, 0.35]]]), True, False).all())
self.assertTrue(np.where(np.round(composite3(('min', 'median'), *ex1, normalize=None), 2) == np.array([[[-9999., -9999.], [-9999., -9999.]], [[-9999., 0.25], [0.4, 0.35]]]), True, False).all())
self.assertTrue(np.where(np.round(composite3(('max', 'median'), *ex1, normalize=None, dtype=int), 2) == np.array([[[-9999., 0.9], [0.4, 0.5]], [[-9999., 0], [0, 0]]]), True, False).all()) # Setting `dtype=int` simulates the type coercion error we want to test
# Should be able to subtract an endmember by setting None
self.assertTrue(np.where(composite2(('max', None), *ex1, normalize=None) == np.array([[[-9999., 0.9], [0.4, 0.5]]]), True, False).all())
self.assertTrue(np.where(composite3(('max', None), *ex1, normalize=None) == np.array([[[-9999., 0.9], [0.4, 0.5]]]), True, False).all())
# When subtracting one endmember from two and normalizing, all non-null
# pixels should sum to one
self.assertTrue(np.where(composite2(('max', None), *ex1, normalize='sum') == np.array([[[-9999., 1], [1, 1]]]), True, False).all())
self.assertTrue(np.where(composite3(('max', None), *ex1, normalize='sum') == np.array([[[-9999., 1], [1, 1]]]), True, False).all())
def test_convex_hull_graham(self):
'''Should correctly calculate a convex hull.'''
result = convex_hull_graham(self.cases['High/Bright'])
self.assertEqual(result, [
(300683, 4707459), (305603, 4676919), (331583, 4690839),
(351143, 4718739),
(318923, 4724259)
])
def test_hall_rectification(self):
'''Should rectify an image in the expected way.'''
control_sets = {
'High/Bright': [(331501.45,4694346.66), (319495.39,4706820.66), (298527.006,4691417.99)],
'Low/Dark': [(322577.40,4658508.99), (361612.79,4694665.62), (378823.69,4692132.56)]
}
ref_raster = gdal.Open(os.path.join(self.test_dir, 'multi7_raster.tiff'))
sub_raster = gdal.Open(os.path.join(self.test_dir, 'multi7_raster2.tiff'))
# NOTE: Using same control sets for reference, subject
hall_rectification(ref_raster, sub_raster, self.test_dir, control_sets, control_sets)
arr, gt, wkt = as_array(os.path.join(self.test_dir, 'rect_multi7_raster2.tiff'))
self.assertTrue(np.array_equal(arr.shape, (6, 74, 81)))
self.assertTrue(np.array_equiv(arr[:,50,50].round(5), np.array([
17, 1331, 1442, 3422, 2916, 2708
]).round(5)))
def test_point_to_pixel_geometry(self):
'''Should correctly calculate pixel geometry from center points.'''
ds = gdal.Open(os.path.join(self.test_dir, 'multi3_raster.tiff'))
xy_coords = pixel_to_xy([[10, 10], [20, 20]], ds.GetGeoTransform(),
ds.GetProjection())
geometries = point_to_pixel_geometry(xy_coords, source_epsg = 32617, target_epsg=4326)
self.assertTrue(isinstance(geometries[0], ogr.Geometry))
self.assertEqual(geometries[0].ExportToWkt(), 'POLYGON ((-84.89202 42.45761 0,-84.89165 42.45763 0,-84.89164 42.45736 0,-84.892 42.45734 0,-84.89202 42.45761 0))')
self.assertEqual(geometries[1].ExportToWkt(), 'POLYGON ((-84.51282 42.19974 0,-84.51246 42.19975 0,-84.51244 42.19948 0,-84.51281 42.19947 0,-84.51282 42.19974 0))')
def test_principle(self):
'''
Demonstrate the critical error that arises in composite() if the
`dtype` argument is not set to `np.float32` when that kind of data is
present.
'''
nodata = -9999
# A 5-band array where the first pixel is -9999 in each band
ex1 = np.array([
[[ -9999, -9999],
[ -9999, -9999]],
[[ -9999, -9999],
[ -9999, -9999]],
[[ -9999, 3.18850607e-01],
[ -9999, -9999]],
[[ -9999, -9999],
[ -9999, 2.08578259e-01]],
[[ -9999, 2.99650848e-01],
[ 3.32257748e-01, 4.80180055e-01]]])
# A 5-band array where the first pixel is NOT -9999 in each band
ex2 = np.array([
[[ -9999, -9999],
[ -9999, -9999]],
[[ -9999, -9999],
[ -9999, -9999]],
[[ -9999, -9999],
[ -9999, -9999]],
[[ -9999, -9999],
[ -9999, 2.08578259e-01]],
[[ 3.18850607e-01, 2.99650848e-01],
[ 3.32257748e-01, 4.80180055e-01]]])
# Array with integer-valued -9999; causes erroneous behavior
ex1a_out = np.round(np.apply_along_axis(lambda xy: np.median([
a for a in xy if a != nodata
]) if not np.where(xy == nodata, True, False).all() else -9999,
0, ex1), 2)
self.assertTrue(np.where(ex1a_out == np.array([
[-9.999e+03, 0],
[ 0, 0]]), True, False).all())
# Array with integer-valued -9999 but doesn't cause error
ex2a_out = np.round(np.apply_along_axis(lambda xy: np.median([
a for a in xy if a != nodata
]) if not np.where(xy == nodata, True, False).all() else -9999,
0, ex2), 3)
self.assertTrue(np.where(ex2a_out == np.array([
[ 0.319, 0.3 ],
[ 0.332, 0.344]]), True, False).all())
# Float-valued -9999.0 fixes error
ex1b_out = np.round(np.apply_along_axis(lambda xy: np.median([
a for a in xy if a != nodata
]) if not np.where(xy == nodata, True, False).all() else -9999.0,
0, ex1), 2)
self.assertTrue(np.where(ex1b_out == np.array([
[-9.999e+03, 3.1e-01],
[ 3.300e-01, 3.4e-01]]), True, False).all())
ex2b_out = np.round(np.apply_along_axis(lambda xy: np.median([
a for a in xy if a != nodata
]) if not np.where(xy == nodata, True, False).all() else -9999.0,
0, ex2), 3)
self.assertTrue(np.where(ex2b_out == np.array([
[ 0.319, 0.3 ],
[ 0.332, 0.344]]), True, False).all())
def test_mnf_rotation(self):
'''MNF rotation should achieve an expected result.'''
file_path = os.path.join(self.test_dir, 'multi7_raster.tiff')
ds = gdal.Open(file_path)
nodata = ds.GetRasterBand(1).GetNoDataValue()
hsi = ds.ReadAsArray().transpose()
hsi[hsi==nodata] = 0
mnf = MNF()
hsi_post_mnf = mnf.apply(hsi)
self.assertEqual(int(hsi_post_mnf.diagonal()[0].sum()), 132)
def test_ppi_endmember_extraction(self):
'''
Extraction by PPI, following an MNF rotation, should achieve an
expected result.
'''
file_path = os.path.join(self.test_dir, 'multi7_raster.tiff')
ds = gdal.Open(file_path)
ppi = PPI()
nodata = ds.GetRasterBand(1).GetNoDataValue()
hsi = ds.ReadAsArray().transpose()
hsi[hsi==nodata] = 0
mnf = MNF()
hsi_post_mnf = mnf.apply(hsi)
members = ppi.extract(hsi_post_mnf[...,0:3], 3, numSkewers=10000)
self.assertEqual(int(members.sum()), 934)
class Transforms(Tester):
test_dir = TEST_DIR
def test_tasseled_cap_transformation(self):
'''
Should apply the tasseled cap transform and produce expected output.
'''
in_fname = os.path.join(self.test_dir, 'multi7_raster.tiff')
out_fname = os.path.join(self.test_dir, 'temp.tiff')
# Load the array; apply the transform
rast, gt, wkt = as_array(in_fname)
tm_image = tasseled_cap_tm(rast, ncomp = 6)
self.assertEqual(tm_image.shape, (6, 74, 81))
self.assertEqual(tm_image.mean().round(0), 212.0)
self.assertTrue((tm_image[:,1,1] == np.array([0, 0, 0, 0, 0, 0])).all())
self.assertTrue((tm_image[:,50,50].round(0) == np.array([ 3269., 2328., -1752., -1140., 407., 18.])).all())
def test_bci_calculation(self):
'''
Should calculate the biophysical composition index and produce
expected output.
'''
in_fname = os.path.join(self.test_dir, 'multi7_raster.tiff')
out_fname = os.path.join(self.test_dir, 'temp2.tiff')
# Load the array; apply the transform
rast, gt, wkt = as_array(in_fname)
bci_image = biophysical_composition_index(rast, tc_func=tasseled_cap_tm)
self.assertEqual(bci_image.shape, (1, 74, 81))
self.assertEqual((bci_image.mean() * 100).round(0), -65.0)
self.assertTrue((bci_image[:,1,1].round(3) == np.array([-0.574])).all())
self.assertTrue((bci_image[:,50,50].round(3) == np.array([-0.702])).all())
class Utilities(Tester):
test_dir = TEST_DIR
def test_combine_dicts(self):
'''
Combining dictionaries of lists should also be idempotent/ with no
side effects.
'''
foo = {'a': [1, 2, 3], 'b': [1]}
bar = {'a': [4], 'b': [2, 3]}
self.assertEqual(combine_dicts(foo, bar), {
'a': [1, 2, 3, 4],
'b': [1, 2, 3]
})
self.assertEqual(foo, {'a': [1, 2, 3], 'b': [1]})
self.assertEqual(bar, {'a': [4], 'b': [2, 3]})
def test_file_raster_and_array_access(self):
'''
Tests that essential file reading and raster/array conversion utilities
are working properly.
'''
from_as_array = as_array(os.path.join(self.test_dir, 'multi3_raster.tiff'))
from_as_raster = as_raster(os.path.join(self.test_dir, 'multi3_raster.tiff'))
self.assertTrue(len(from_as_array) == len(from_as_raster) == 3)
self.assertTrue(isinstance(from_as_array[0], np.ndarray))
self.assertTrue(isinstance(from_as_raster[0], gdal.Dataset))
def test_array_to_raster_interface(self):
'''
The array_to_raster() and array_to_raster_clone functions should
perform as expected.
'''
# First round
ds = gdal.Open(os.path.join(self.test_dir, 'multi3_raster.tiff'))
gt = ds.GetGeoTransform()
wkt = ds.GetProjection()
arr = ds.ReadAsArray()
ds = None
rast = array_to_raster(arr, gt, wkt)
self.assertEqual(gt, rast.GetGeoTransform())
self.assertEqual(wkt, rast.GetProjection())
# Second round
rast2 = array_to_raster_clone(arr, os.path.join(self.test_dir,
'multi7_raster.tiff'))
self.assertEqual(gt, rast2.GetGeoTransform())
self.assertEqual(wkt, rast2.GetProjection())
def test_binary_mask(self):
'''
Should successfully apply a binary mask to an image.
'''
out_fname = os.path.join(self.test_dir, 'temp.tiff')
rast, gt, wkt = as_array(os.path.join(self.test_dir, 'multi7_raster.tiff'))
masked = binary_mask(rast, density_slice(rast[4, ...]))
self.assertEqual(masked.shape, (6, 74, 81))
self.assertEqual(masked.mean().round(0), -3627.0)
self.assertTrue((masked[:,15,19] == np.array([348, 576, 374, 3159, 2033, 2386])).all())
def test_saturation_mask(self):
'''
Should successfully mask out saturated pixels.
'''
out_fname = os.path.join(self.test_dir, 'temp.tiff')
rast, gt, wkt = as_array(os.path.join(self.test_dir, 'multi7_raster.tiff'))
mask = saturation_mask(rast)
masked = binary_mask(rast, mask)
self.assertEqual(masked.shape, (6, 74, 81))
self.assertEqual(masked.mean().round(0), 881.0)
self.assertTrue((masked[:,15,19] == np.array([348, 576, 374, 3159, 2033, 2386])).all())
def test_density_slicing(self):
'''
Density slicing with the density_slice() function should perform as
expected.
'''
rast, gt, wkt = as_array(os.path.join(self.test_dir, 'multi7_raster.tiff'))
r1 = density_slice(rast[4, ...])
r2 = density_slice(rast[4, ...], threshold=500)
r3 = density_slice(rast[4, ...], rel=np.greater_equal)
self.assertEqual(r1.shape, (74, 81))
self.assertEqual(r2.shape, (74, 81))
self.assertEqual(r3.shape, (74, 81))
self.assertEqual(np.count_nonzero(r1), 2675)
self.assertEqual(np.count_nonzero(r2), 2642)
self.assertEqual(np.count_nonzero(r3), 3319)
def test_lat_lng_to_pixel_and_reverse(self):
'''
Conversions from map coordinates to pixel coordinates and vice-versa
should not introduce pixel-position errors of greater than 1 pixel.
'''
coords = [(random.randint(1, 70),
random.randint(1, 70)) for x in range(1, 10)]
file_path = os.path.join(self.test_dir, 'multi3_raster.tiff')
p2ll = pixel_to_xy(coords, path=file_path, dd=True)
ll2p = xy_to_pixel(p2ll, path=file_path, dd=True)
self.assertTrue(np.max(np.abs((np.array(ll2p) - np.array(coords)))) <= 1)
def test_lat_lng_to_pixel_and_reverse_without_reference_dataset(self):
'''
Should accurately locate the lat-long coordinates of two pixel
coordinate pairs.
'''
coords = ((-84.5983, 42.7256), (-85.0807, 41.1138))
pixels = [(18, 0), (2, 59)]
file_path = os.path.join(self.test_dir, 'multi3_raster.tiff')
ds = gdal.Open(file_path)
gt = ds.GetGeoTransform()
wkt = ds.GetProjection()
p2ll = pixel_to_xy(pixels, gt=gt, wkt=wkt, dd=True)
ll2p = xy_to_pixel(coords, gt=gt, wkt=wkt, dd=True)
self.assertEqual(ll2p, [(18, 0), (2, 59)])
self.assertTrue(np.max(np.abs((np.array(p2ll) - np.array(coords)))) <= 1)
def test_spectral_profile(self):
'''
Should correctly retrieve a spectral profile from a raster dataset.
'''
coords = ((-84.5983, 42.7256), (-85.0807, 41.1138))
pixels = [(18, 0), (2, 59)]
file_path = os.path.join(self.test_dir, 'multi3_raster.tiff')
ds = gdal.Open(file_path)
kwargs = {
'gt': ds.GetGeoTransform(),
'wkt': ds.GetProjection(),
'dd': True
}
# The true spectral profile
spectra = np.array([[237, 418, 325], [507, 616, 445]], dtype=np.int16)
sp1 = spectra_at_xy(ds, coords, **kwargs)
sp2 = spectra_at_xy(ds.ReadAsArray(), coords, **kwargs)
sp3 = spectra_at_idx(ds.ReadAsArray().transpose(), pixels)
self.assertEqual(spectra.tolist(), sp1.tolist())
self.assertEqual(spectra.tolist(), sp2.tolist())
self.assertEqual(spectra.tolist(), sp3.tolist())
def test_masking(self):
'''
Masking should go on without a hitch and the result should be just
as expected.
'''
file_path = os.path.join(self.test_dir, 'multi7_raster.tiff')
ds = gdal.Open(file_path)
raw_mask = gdal.Open(os.path.join(self.test_dir, 'multi7_mask.tiff'))
mask = cfmask(raw_mask, nodata=-9999)
masked = binary_mask(ds.ReadAsArray(), mask)
self.assertEqual(ds.ReadAsArray().diagonal()[0,0], 0)
self.assertEqual(masked.diagonal()[0,0], -9999)
class Visualize(Tester):
test_dir = TEST_DIR
def test_interactive_feature_space(self):
'''Tests that interactive plotting works as expected.'''
path = os.path.join(self.test_dir, 'multi7_raster.tiff')
vis = FeatureSpace(path, mask = None, cut_dim = 2, transform = True,
nodata = -9999, epsg = 32617, keyword = 'test')
fig = vis.plot_feature_space(c = None, interact = True, hold = True)
vis.on_reset()
vis.x0 = 161.56909597769811
vis.y0 = 140.10416666666663
vis.x1 = 171.02747909199519
vis.y1 = 133.83487654320987
vis.on_draw(output_dir = self.test_dir)
with open(os.path.join(self.test_dir, 'FeatureSpace_selection_test.kml'), 'r') as stream:
ref_contents = stream.read()
with open(os.path.join(self.test_dir, 'FeatureSpace_selection_test_1.kml'), 'r') as stream:
file_contents = stream.read()
self.assertEqual(ref_contents, file_contents)
if __name__ == '__main__':
unittest.main()
