test_geoarray.py 26.6 KB
Newer Older
1
2
#!/usr/bin/env python
# -*- coding: utf-8 -*-
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

# geoarray, A fast Python interface for image geodata - either on disk or in memory.
#
# Copyright (C) 2019  Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de)
#
# This software was developed within the context of the GeoMultiSens project funded
# by the German Federal Ministry of Education and Research
# (project grant code: 01 IS 14 010 A-C).
#
# This program is free software: you can redistribute it and/or modify it under
# the terms of the GNU Lesser General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
# details.
#
# You should have received a copy of the GNU Lesser General Public License along
# with this program.  If not, see <http://www.gnu.org/licenses/>.

25
26
"""
test_geoarray
27
-------------
28

29
Tests for the geoarray.baseclasses module.
30
"""
31
32
33
34
35


from unittest import TestCase
import tempfile
import os
36
from collections import OrderedDict
37
38
from typing import Iterable

39
import numpy as np
40
from shapely.geometry import Polygon
41
from osgeo.osr import SpatialReference  # noqa
42
import matplotlib
43
from parameterized import parameterized  # noqa
44
from py_tools_ds.geo.vector import geometry
45
from py_tools_ds.geo.coord_trafo import transform_any_prj
46

47
from geoarray import GeoArray, masks, __file__ as geoarray_root  # noqa
48
49
50
51
52
from geoarray.metadata import GDAL_Metadata


__author__ = ["Daniel Scheffler", "Jessica Palka"]

53

54
# Path of the tests-directory in the geoarray-package.
55
tests_path = os.path.abspath(os.path.join(geoarray_root, "..", ".."))
56
path_data = os.path.join(tests_path, "tests", "data", "L8_2bands_extract10x11.tif")
57

58
59
60
61
# disable matplotlib figure popups
matplotlib.use('Template')  # disables matplotlib figure popups


62
63
64
65
66
67
68
69
70
def _get_gA_inMem_notInMem():
    gA_inMem = GeoArray(path_data)
    gA_inMem.to_mem()
    gA_inMem.filePath = ""

    gA_notInMem = GeoArray(path_data)

    params_inMem_NotInMem = [("inMem", gA_inMem),
                             ("notInMem", gA_notInMem)]
71

72
    return params_inMem_NotInMem
73

74
75
76
77

class Test_GeoArray(TestCase):
    # expected results
    given_bandnames = ['B1', 'B2']
78
79

    @classmethod
80
    def setUpClass(cls) -> None:
81
        cls.gA = _get_gA_inMem_notInMem()[0][1]
82

83
    @parameterized.expand(_get_gA_inMem_notInMem())
84
85
86
87
    def test_bandnames(self, _, gA):
        """Test GeoArray.bandnames."""
        self.assertEqual(gA.bandnames, OrderedDict([('B1', 0), ('B2', 1)]))

88
    @parameterized.expand(_get_gA_inMem_notInMem())
89
90
91
    def test_shape(self, _, gA):
        """Test GeoArray.shape, .ndim, .rows, .columns, and .bands"""
        shape_attrs = (gA.ndim, gA.rows, gA.columns, gA.bands)
92
        shape_property = ('DIMENSIONS', 'ROWS', 'COLUMNS', 'BANDS')
93
94
95
96
97
98
99
100
101
102
103
        R_exp, C_exp, B_exp = 10, 11, 2
        expected_result = (3, R_exp, C_exp, B_exp)  # dimensions, rows, columns, band

        self.assertEqual(gA.shape, (R_exp, C_exp, B_exp))

        for i in range(4):
            with self.subTest(i=i):
                self.assertEqual(shape_attrs[i], expected_result[i],
                                 msg='The number of {i} is different from the expected result!'.format(
                                     i=shape_property[i]))

104
    @parameterized.expand(_get_gA_inMem_notInMem())
105
106
107
108
    def test_dtype(self, _, gA):
        """Test GeoArray.dtype."""
        self.assertEqual(gA.dtype, np.dtype('float32'))

109
    @parameterized.expand(_get_gA_inMem_notInMem())
110
111
112
113
114
115
116
    def test_geotransform_gsd_grid(self, _, gA):
        """Test GeoArray.geotransform, .xgsd, .ygsd, and .xygrid_specs."""
        self.assertEqual(gA.geotransform, [365985.0, 30.0, 0.0, 5916615.0, 0.0, -30.0])
        self.assertEqual((gA.xgsd, gA.ygsd), (30, 30))
        self.assertEqual(gA.xygrid_specs, [[365985.0, 366015.0], [5916615.0, 5916585.0]],
                         msg='The [[xOrigin, xGSD], [yOrigin, yGSD]]-grid is not as expected!')

117
    @parameterized.expand(_get_gA_inMem_notInMem())
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
    def test_projection_epsg(self, _, gA):
        """Test GeoArray.projection and .epsg."""

        # Convert WKT-string of the projection to a Proj4_string
        srs = SpatialReference()
        srs.ImportFromWkt(gA.projection)
        proj4 = srs.ExportToProj4()

        self.assertEqual(proj4.strip(' /t/n/r'), '+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs')
        self.assertEqual(gA.epsg, 32633)

    def test_nodata_not_given(self):
        """Test GeoArray.nodata in case no nodata value is given (i.e., auto-detected)."""
        gA = GeoArray(np.random.randint(1, 10, (5, 10)))
        self.assertIsNone(gA._nodata)
        self.assertEqual(self.gA.nodata, -9999.0)

    def test_nodata_given(self):
        """Test GeoArray.nodata in case nodata value is given."""
        gA = GeoArray(np.array([1, 2]), nodata=-9999)
        self.assertEqual(gA.nodata, -9999.0)

140
141
142
    def test_nodata_out_of_range(self):
        """Test if given nodata value is valid with respect to the array data type."""
        # test nodata value in range
Daniel Scheffler's avatar
Daniel Scheffler committed
143
        GeoArray(np.array([1, 2]).astype(np.uint8), nodata=255)
144
145
146

        # test nodata value out of range in range
        with self.assertRaises(ValueError):
Daniel Scheffler's avatar
Daniel Scheffler committed
147
            GeoArray(np.array([1, 2]).astype(np.uint8), nodata=256)
148

149
        # test float nodata value together with integer array data type
Daniel Scheffler's avatar
Daniel Scheffler committed
150
        GeoArray(np.array([1, 2]).astype(np.uint8), nodata=0.)
151

152
    @parameterized.expand(_get_gA_inMem_notInMem())
153
154
    def test_calc_mask_nodata(self, _, gA):
        def test_with_flag_all(_gA):
155
            """Test mask output with flag='all' (flag when nodata in ALL bands)."""
156
157
            mnd_all = _gA.calc_mask_nodata(overwrite=True, flag='all')
            mnd_all_exp = np.ones((10, 11)).astype(bool)
158
159
160
            mnd_all_exp[0, 0] = False  # False represents nodata
            self.assertTrue(np.array_equal(mnd_all, mnd_all_exp))

161
        def test_with_flag_any(_gA):
162
            """Test mask output with flag='any' (flag when nodata in ANY band)."""
163
164
            mnd_any = _gA.calc_mask_nodata(overwrite=True, flag='any')
            mnd_any_exp = np.ones((10, 11)).astype(bool)
165
166
167
            mnd_any_exp[0, 0:3] = False  # False represents nodata
            self.assertTrue(np.array_equal(mnd_any, mnd_any_exp))

168
169
        if gA.is_inmem:
            _real_nodata = gA.nodata
170
171
            try:
                for nodata in [_real_nodata, np.nan]:
172
                    gA.nodata = nodata
173

174
175
176
                    gA[0, 0, :] = nodata  # set pixel [0, 0] to nodata in all bands
                    gA[0, 1, 0] = nodata  # set pixel [0, 1] to nodata in first band only
                    gA[0, 2, 1] = nodata  # set pixel [0, 2] to nodata in second band only
177

178
179
                    test_with_flag_all(gA)
                    test_with_flag_any(gA)
180
181

            finally:
182
                gA.nodata = _real_nodata
183
184

            with self.assertRaises(ValueError):
185
                gA.calc_mask_nodata(overwrite=True, flag='bad_string')
186

187
    @parameterized.expand(_get_gA_inMem_notInMem())
188
    def test___getitem__(self, _, gA):
189
190
191
192
        def validate(array, exp_shape):
            self.assertIsInstance(array, np.ndarray)
            self.assertEqual(array.shape, exp_shape)

193
        R, C, B = gA.shape  # (10, 11, 2)
194
195

        # test row/col subset
196
197
198
199
        validate(gA[:1, :3, :], (1, 3, B))  # only one row is requested, given as a slice
        validate(gA[0, :3, :], (3, B))  # only one row is requested, given as an int
        validate(gA[2:5, :3], (3, 3, B))  # third dimension is not given
        validate(gA[2:5, :3, :], (3, 3, B))
200
201

        # test band subset
202
203
204
205
        validate(gA[:, :, 0:1], (R, C, 1))  # band slice  # returns 3D array
        validate(gA[:, :, 0], (R, C))  # band indexing  # returns 2D array
        validate(gA[1], (R, C))  # only band is given  # returns 2D
        validate(gA['B1'], (R, C))  # only bandname is given
206

207
208
209
210
211
        # test read point data in case only data of a single band is needed
        validate(gA[[1, 5], [2, 4], 1], (2, ))
        if not gA.is_inmem:
            self.assertNotEqual(gA._arr_cache['arr_cached'].shape, gA.shape)

212
        # test wrong inputs
213
        self.assertRaises(ValueError, gA.__getitem__, 'B01')
214

215
216
        # test full array  # NOTE: This sets self.gA.arr!
        validate(gA[:], (R, C, B))
217

218
219
        # TODO: add further tests

220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
    def test___getitem__consistency(self):
        testarr = np.zeros((2, 2, 2))
        testarr[0, :, :] = [[11, 12], [13, 14]]
        testarr[1, :, :] = [[21, 22], [23, 24]]

        gA_inmem = GeoArray(testarr)
        inmem_res = gA_inmem[0, :, :]

        with tempfile.TemporaryDirectory() as tf:
            gA_inmem.save(os.path.join(tf, 'test.bsq'))

            gA_notinmem = GeoArray(os.path.join(tf, 'test.bsq'))
            notinmem_res = gA_notinmem[0, :, :]

        self.assertEqual(inmem_res.ndim, notinmem_res.ndim)
        self.assertEqual(inmem_res.shape, notinmem_res.shape)

Daniel Scheffler's avatar
Bugfix.    
Daniel Scheffler committed
237
    def test___getitem__consistency_3d_array_1_column(self):
238
        testarr = np.array([[1, 2], [3, 4]]).reshape((2, 1, 2))
Daniel Scheffler's avatar
Bugfix.    
Daniel Scheffler committed
239
240
241
242
243
244
245
246
247
248
249
250
251

        gA_inmem = GeoArray(testarr)
        inmem_res = gA_inmem[:]

        with tempfile.TemporaryDirectory() as tf:
            gA_inmem.save(os.path.join(tf, 'test.bsq'))

            gA_notinmem = GeoArray(os.path.join(tf, 'test.bsq'))
            notinmem_res = gA_notinmem[:]

        self.assertEqual(inmem_res.ndim, notinmem_res.ndim)
        self.assertEqual(inmem_res.shape, notinmem_res.shape)

252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
    def test___getitem__consistency_2d_array(self):
        testarr = np.zeros((2, 2))
        testarr[:, :] = [[11, 12], [13, 14]]

        gA_inmem = GeoArray(testarr)
        inmem_res = gA_inmem[0, 0]

        with tempfile.TemporaryDirectory() as tf:
            gA_inmem.save(os.path.join(tf, 'test.bsq'))

            gA_notinmem = GeoArray(os.path.join(tf, 'test.bsq'))
            notinmem_res = gA_notinmem[0, 0]

        self.assertEqual(inmem_res.ndim, notinmem_res.ndim)
        self.assertEqual(inmem_res.shape, notinmem_res.shape)

268
    def test_numpy_array(self):
269
        arr = np.array(self.gA)
270
        self.assertIsInstance(arr, np.ndarray)
271
        self.assertEqual(arr.shape, self.gA.shape)
272

273
    def test_box(self):
274
275
        """Test GeoArray.box."""
        self.assertIsInstance(self.gA.box, geometry.boxObj)
276

277
    def test_mask_nodata(self):
278
279
280
        # TODO: Consider the dependency of mask_nodata on calc_mask_nodata().
        """Test GeoArray.mask_nodata."""
        self.assertIsInstance(self.gA.mask_nodata, masks.NoDataMask)
281

282
    def test_mask_baddata(self):
283
284
285
286
        """Test GeoArray.mask_baddata."""
        bdm = masks.BadDataMask(os.path.join(tests_path, "tests", "data", "L8_BadDataMask10x11.tif"))

        for i in range(2):
287
288
            with self.subTest(i=1):
                if i == 0:
289
290
                    # mask_baddata-function is "None", when it is not actively set
                    self.assertIsNone(self.gA.mask_baddata)
291
                if i == 1:
292
293
                    self.gA.mask_baddata = bdm
                    self.assertIsInstance(self.gA.mask_baddata, masks.BadDataMask)
294

295
    def test_footprint_poly(self):
296
297
        # TODO: Test the validation of the footprint_poly-function.
        # TODO: Consider the dependencies of the footprint_poly-function on mask_nodata, boxObj.
298
299
        """Test GeoArray.footprint_poly."""
        self.assertIsInstance(self.gA.footprint_poly, Polygon)
300

301
    def test_metadata(self):
302
        # TODO: Create a metadata-file for the tested TIFF-Image.
303
        # TODO: Test, if the metadata-function gives an output
304
305
        """Test GeoArray.metadata."""
        self.assertIsInstance(self.gA.metadata, GDAL_Metadata)
306

307
    def test_tiles(self):
308
309
        test_gAs = [self.gA,  # 3D
                    GeoArray(self.gA[:, :, 0], geotransform=self.gA.gt, projection=self.gA.prj)]  # 2D
310

311
312
313
314
315
        for gA in test_gAs:
            tiles = gA.tiles(tilesize=(50, 50))
            self.assertIsInstance(tiles, Iterable)

            for ((rS, rE), (cS, cE)), tile in tiles:
316
                self.assertTrue(np.array_equal(tile, gA[rS: rE + 1, cS: cE + 1]))
317

318
    def test_get_subset_3D_geoarray(self):
319
        # test without resetting band names
320
        sub_gA = self.gA.get_subset(xslice=slice(2, 5), yslice=slice(None, 3), zslice=slice(1, 2))
321
        self.assertIsInstance(sub_gA, GeoArray)
322
        self.assertTrue(list(sub_gA.bandnames), list(self.gA.bandnames)[1])
323

324
        # test with providing only xslice
325
        sub_gA = self.gA.get_subset(xslice=slice(2, 5))
326
327
328
        self.assertIsInstance(sub_gA, GeoArray)

        # test with providing only yslice
329
        sub_gA = self.gA.get_subset(yslice=slice(None, 3))
330
331
        self.assertIsInstance(sub_gA, GeoArray)

332
        # test with zslice provided as list
333
        sub_gA = self.gA.get_subset(xslice=slice(2, 5), yslice=slice(None, 3), zslice=[0, 1])
334
335
        self.assertIsInstance(sub_gA, GeoArray)

336
        # test without providing zslice
337
        sub_gA = self.gA.get_subset(xslice=slice(2, 5), yslice=slice(None, 3))
338
339
        self.assertIsInstance(sub_gA, GeoArray)

340
        # test requesting only one column
341
        sub_gA = self.gA.get_subset(xslice=slice(0, 1), yslice=slice(None, 3))
342
343
        self.assertIsInstance(sub_gA, GeoArray)

344
        # test with resetting band names
345
346
        sub_gA = self.gA.get_subset(xslice=slice(2, 5), yslice=slice(None, 3), zslice=slice(1, 2),
                                    reset_bandnames=True)
347
348
        self.assertTrue(list(sub_gA.bandnames), ['B1'])

349
        # test arrays are equal
350
351
        sub_gA = self.gA.get_subset(xslice=slice(2, 5), yslice=slice(None, 3), zslice=slice(1, 2))
        sub_arr = self.gA[:3, 2:5, 1:2]
352
        self.assertTrue(np.array_equal(sub_gA[:], sub_arr))
353

354
        # test deepcopied arrays (modification of sub_gA.arr must not affect self.gA.arr)
355
        sub_gA[:2, :2] = 99
356
357
358
        self.assertTrue(np.array_equal(sub_gA[:2, :2], np.full((2, 2, 1), 99, self.gA.dtype)))
        self.assertNotEqual(np.mean(sub_arr[:2, :2]), 99)
        self.assertNotEqual(np.std(sub_arr[:2, :2]), 0)
359
360
361
362
363

        # test metadata
        self.assertEqual(sub_gA.meta.bands, 1)
        self.assertEqual(len(list(sub_gA.meta.band_meta.values())[0]), 1)
        self.assertEqual(len(list(sub_gA.bandnames.keys())), 1)
364
365
        self.assertNotEqual(sub_gA.gt, self.gA.gt)
        self.assertEqual(sub_gA.prj, self.gA.prj)
366

367
        # test not to return GeoArray
368
369
        out = self.gA.get_subset(xslice=slice(2, 5), yslice=slice(None, 3), zslice=slice(1, 2),
                                 return_GeoArray=False)
370

371
372
373
374
375
376
        self.assertIsInstance(out, tuple)
        self.assertTrue(len(out), 3)
        self.assertIsInstance(out[0], np.ndarray)
        self.assertIsInstance(out[1], tuple)
        self.assertIsInstance(out[2], str)

377
    def test_get_subset_2D_geoarray(self):
378
        gA_2D = GeoArray(self.gA[0])
379
380

        # test without resetting band names
381
        sub_gA = self.gA.get_subset(xslice=slice(2, 5), yslice=slice(None, 3))
382
        self.assertIsInstance(sub_gA, GeoArray)
383
        self.assertTrue(list(sub_gA.bandnames), list(self.gA.bandnames)[1])
384
385
386
387
388
389
390
391
392
393
394
395
396

        # test with providing only xslice
        sub_gA = gA_2D.get_subset(xslice=slice(2, 5))
        self.assertIsInstance(sub_gA, GeoArray)

        # test with providing only yslice
        sub_gA = gA_2D.get_subset(yslice=slice(None, 3))
        self.assertIsInstance(sub_gA, GeoArray)

        # test without providing zslice
        sub_gA = gA_2D.get_subset(xslice=slice(2, 5), yslice=slice(None, 3))
        self.assertIsInstance(sub_gA, GeoArray)

397
        # test requesting only one column
398
        sub_gA = self.gA.get_subset(xslice=slice(0, 1), yslice=slice(None, 3))
399
400
        self.assertIsInstance(sub_gA, GeoArray)

401
402
403
404
        # test with resetting band names
        sub_gA = gA_2D.get_subset(xslice=slice(2, 5), yslice=slice(None, 3), reset_bandnames=True)
        self.assertTrue(list(sub_gA.bandnames), ['B1'])

405
406
407
408
409
        # test arrays are equal
        sub_gA = gA_2D.get_subset(xslice=slice(2, 5), yslice=slice(None, 3))
        sub_gA_2D = gA_2D[:3, 2:5]
        self.assertTrue(np.array_equal(sub_gA[:], sub_gA_2D))

410
        # test deepcopied arrays (modification of sub_gA.arr must not affect self.gA.arr)
411
412
413
414
415
416
417
418
419
420
421
422
        sub_gA[:2, :2] = 99
        self.assertTrue(np.array_equal(sub_gA[:2, :2], np.full((2, 2), 99, gA_2D.dtype)))
        self.assertNotEqual(np.mean(sub_gA_2D[:2, :2]), 99)
        self.assertNotEqual(np.std(sub_gA_2D[:2, :2]), 0)

        # test metadata
        self.assertEqual(sub_gA.meta.bands, 1)
        self.assertEqual(len(list(sub_gA.meta.band_meta.values())[0]), 1)
        self.assertEqual(len(list(sub_gA.bandnames.keys())), 1)
        self.assertNotEqual(sub_gA.gt, gA_2D.gt)
        self.assertEqual(sub_gA.prj, gA_2D.prj)

423
424
425
426
427
428
        # test not to return GeoArray
        out = gA_2D.get_subset(xslice=slice(2, 5), yslice=slice(None, 3), return_GeoArray=False)

        # test with provided zslice
        with self.assertRaises(ValueError):
            gA_2D.get_subset(xslice=slice(2, 5), yslice=slice(None, 3), zslice=slice(1, 2))
429
430
        with self.assertRaises(ValueError):
            gA_2D.get_subset(xslice=slice(2, 5), yslice=slice(None, 3), zslice=[1, 2])
431
432
433
434
435
436
437

        self.assertIsInstance(out, tuple)
        self.assertTrue(len(out), 3)
        self.assertIsInstance(out[0], np.ndarray)
        self.assertIsInstance(out[1], tuple)
        self.assertIsInstance(out[2], str)

438
    @parameterized.expand(_get_gA_inMem_notInMem())
439
440
    def test_save(self, _, gA):
        """Test GeoArray.save."""
441
442
        with tempfile.TemporaryDirectory(dir=os.path.join(tests_path, "tests", "data", "output"),
                                         prefix='test_save_') as td:
443
444
445

            # save GeoArray and validate output
            for fmt in ['ENVI', 'GTiff']:
446
447
448
                outpath = os.path.join(td, "TestGeoArray_10x11_copy.tif")
                gA.save(outpath, fmt=fmt)
                self.assertTrue(os.path.exists(outpath))
449

450
                gA_copy = GeoArray(outpath)
451
452
453
                self.assertIsInstance(gA_copy, GeoArray)
                self.assertTrue(np.array_equal(gA[:], gA_copy[:]))
                self.assertEqual(gA._nodata, gA_copy._nodata)
454

455
456
            # TODO validate written metadata (inMem, notInMem, in case a notInMem dataset already has custom metadata)

457
458
    def test_show(self):
        # test 3D case
459
        self.gA.show()
460
        self.gA.show(interactive=True)  # only works if test is started with ipython.
461

462
463
464
465
466
467
468
        # test 2D case
        gA = GeoArray(self.gA[:, :, 0])
        gA.show()
        with self.assertWarnsRegex(UserWarning, 'Currently there is no interactive mode for single-band arrays.'):
            gA.show(interactive=True)

    def test_show_map(self):
469
        self.gA.show_map()
470

471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
    def test_show_map_noepsg(self):
        wkt_noepsg = \
            """
            PROJCRS["BU MEaSUREs Lambert Azimuthal Equal Area - SA - V01",
                BASEGEOGCRS["WGS 84",
                    DATUM["World Geodetic System 1984",
                        ELLIPSOID["WGS 84",6378137,298.257223563,
                            LENGTHUNIT["metre",1]]],
                    PRIMEM["Greenwich",0,
                        ANGLEUNIT["degree",0.0174532925199433]],
                    ID["EPSG",4326]],
                CONVERSION["unnamed",
                    METHOD["Lambert Azimuthal Equal Area",
                        ID["EPSG",9820]],
                    PARAMETER["Latitude of natural origin",-15,
                        ANGLEUNIT["degree",0.0174532925199433],
                        ID["EPSG",8801]],
                    PARAMETER["Longitude of natural origin",-60,
                        ANGLEUNIT["degree",0.0174532925199433],
                        ID["EPSG",8802]],
                    PARAMETER["False easting",0,
                        LENGTHUNIT["metre",1],
                        ID["EPSG",8806]],
                    PARAMETER["False northing",0,
                        LENGTHUNIT["metre",1],
                        ID["EPSG",8807]]],
                CS[Cartesian,2],
                    AXIS["easting",east,
                        ORDER[1],
                        LENGTHUNIT["metre",1]],
                    AXIS["northing",north,
                        ORDER[2],
                        LENGTHUNIT["metre",1]]]
            """
        wkt_noepsg = ' '.join(wkt_noepsg.split())

        gA = GeoArray(self.gA[:], self.gA.gt, wkt_noepsg)
        gA.show_map()

510
    def test_show_histogram(self):
511
        self.gA.show_histogram()
512

513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
    def test_show_xprofile(self):
        nd_old = self.gA._nodata
        self.gA._nodata = self.gA[0, 0, 0]
        self.gA.show_xprofile(0, 0)
        self.gA.show_xprofile(0, 0, show_nodata=False)
        self.gA._nodata = nd_old

    def test_show_yprofile(self):
        self.gA.show_yprofile(0, 0)
        self.gA.show_yprofile(0, 0, show_nodata=False)

    def test_show_zprofile(self):
        self.gA.show_zprofile(0, 0)
        self.gA.show_zprofile(0, 0, show_nodata=False)
        self.gA.show_zprofile(0, 0, ylim=(0, 5))

        with self.assertRaises(RuntimeError):
            GeoArray(np.random.randint(1, 10, [10, 10])).show_zprofile(0, 0)

        if 'wavelength' not in self.gA.meta.band_meta:
            self.gA.meta.band_meta['wavelength'] = range(self.gA.bands)
            self.gA.show_zprofile(0, 0, ylim=(0, 5))
            del self.gA.meta.band_meta['wavelength']

537
538
    def test_get_mapPos(self):
        # TODO: validate results (pixel dimensions, values, metadata subset, ...)
539
540
541
542
543
        xmin, xmax, ymin, ymax = self.gA.box.boundsMap
        xgsd = self.gA.xgsd
        ygsd = self.gA.ygsd
        arr_sub, gt, prj = self.gA.get_mapPos(mapBounds=(xmin + xgsd, ymin, xmax + xgsd, ymax - ygsd / 2),
                                              mapBounds_prj=self.gA.prj)
544
545
546
        self.assertIsInstance(arr_sub, np.ndarray)
        self.assertIsInstance(gt, (tuple, list))
        self.assertIsInstance(prj, str)
547

548
    def test_read_pointData__singleCoord_oneband(self):
549
550
        vals_z = self.gA.read_pointData(mapXY_points=np.array([[366015, 5916585]]),
                                        band=1)
551

552
        self.assertIsInstance(vals_z, self.gA.dtype)
553
554
555
        self.assertEqual(vals_z, 5766)

    def test_read_pointData__singleCoord_allbands(self):
556
        vals_z = self.gA.read_pointData(mapXY_points=np.array([[366015, 5916585]]))
557
558

        self.assertIsInstance(vals_z, np.ndarray)
559
        self.assertEqual(vals_z.size, self.gA.bands)
560
561
562
        self.assertTrue(np.array_equal(vals_z, np.array([[5352, 5766]])))

    def test_read_pointData__singleCoord_allbands_otherPrj(self):
563
564
565
        X, Y = transform_any_prj(self.gA.epsg, 4326, 366015.1, 5916584.9)
        vals_z = self.gA.read_pointData(mapXY_points=np.array([[X, Y]]),
                                        mapXY_points_prj=4326)
566
567

        self.assertIsInstance(vals_z, np.ndarray)
568
        self.assertEqual(vals_z.size, self.gA.bands)
569
570
571
        self.assertTrue(np.array_equal(vals_z, np.array([[5352, 5766]])))

    def test_read_pointData__multiCoords_allbands(self):
572
573
        vals_z = self.gA.read_pointData(mapXY_points=np.array([[366015, 5916585],
                                                               [366165, 5916400]]))
574
575

        self.assertIsInstance(vals_z, np.ndarray)
576
        shape_exp = (2, 1, self.gA.bands)
577
578
579
580
581
        self.assertEqual(vals_z.shape, shape_exp)
        self.assertTrue(np.array_equal(vals_z, np.array([[5352, 5766],
                                                         [5080, 5436]]).reshape(shape_exp)))

    def test_read_pointData__multiCoords_allbands_otherPrj(self):
582
        X, Y = transform_any_prj(self.gA.epsg, 4326,
583
584
                                 np.array([366015.1, 366165]),
                                 np.array([5916584.9, 5916399.9]))
585
586
        vals_z = self.gA.read_pointData(mapXY_points=np.array([[X, Y]]),
                                        mapXY_points_prj=4326)
587
588

        self.assertIsInstance(vals_z, np.ndarray)
589
        shape_exp = (2, 1, self.gA.bands)
590
591
        self.assertEqual(vals_z.shape, shape_exp)
        self.assertTrue(np.array_equal(vals_z, np.array([[5352, 5766],
Daniel Scheffler's avatar
Daniel Scheffler committed
592
593
                                                         [5080, 5436]]).reshape(shape_exp)),
                        'vals_z: %s' % repr(vals_z))
594
595

    def test_read_pointData__singleCoord_offside_oneband(self):
596
597
598
        vals_z = self.gA.read_pointData(mapXY_points=np.array([[365900, 5916585]]),
                                        band=1,
                                        offside_val=-9999)
599
600
601
602
603

        self.assertIsInstance(vals_z, int)
        self.assertEqual(vals_z, -9999)

    def test_read_pointData__multiCoords_offside_oneband(self):
604
605
606
607
608
609
610
611
        xtoolarge = self.gA.box.boundsMap[1] + self.gA.xgsd
        vals_z = self.gA.read_pointData(mapXY_points=np.array([[365900, 5916585],  # x too small
                                                               [365900, 5917000],  # x too small, y to large
                                                               [xtoolarge, 5916585],  # x too large
                                                               [366015, 5916585],
                                                               [366165, 5916400]]),
                                        offside_val=-9999,
                                        band=1)
612
613
614
615
616
617
618
619
620
621
622

        self.assertIsInstance(vals_z, np.ndarray)
        shape_exp = (5, 1)
        self.assertEqual(vals_z.shape, shape_exp)
        self.assertTrue(np.array_equal(vals_z, np.array([[-9999],
                                                         [-9999],
                                                         [-9999],
                                                         [5766],
                                                         [5436]])))

    def test_read_pointData__multiCoords_offside_allbands(self):
623
624
625
626
627
628
629
        xtoolarge = self.gA.box.boundsMap[1] + self.gA.xgsd
        vals_z = self.gA.read_pointData(mapXY_points=np.array([[365900, 5916585],  # x too small
                                                               [365900, 5917000],  # x too small, y to large
                                                               [xtoolarge, 5916585],  # x too large
                                                               [366015, 5916585],
                                                               [366165, 5916400]]),
                                        offside_val=-9999)
630
631

        self.assertIsInstance(vals_z, np.ndarray)
632
        shape_exp = (5, 1, self.gA.bands)
633
634
635
636
637
638
        self.assertEqual(vals_z.shape, shape_exp)
        self.assertTrue(np.array_equal(vals_z, np.array([[-9999, -9999],
                                                         [-9999, -9999],
                                                         [-9999, -9999],
                                                         [5352, 5766],
                                                         [5080, 5436]]).reshape(shape_exp)))