test_sensormapgeo.py 18.4 KB
Newer Older
Daniel Scheffler's avatar
Daniel Scheffler committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
#!/usr/bin/env python
# -*- coding: utf-8 -*-

# sensormapgeo, A package for transforming remote sensing images between sensor and map geometry.
#
# Copyright (C) 2020  Daniel Scheffler (GFZ Potsdam, danschef@gfz-potsdam.de)
#
# This software was developed within the context of the EnMAP project supported
# by the DLR Space Administration with funds of the German Federal Ministry of
# Economic Affairs and Energy (on the basis of a decision by the German Bundestag:
# 50 EE 1529) and contributions from DLR, GFZ and OHB System AG.
#
# 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/>.

"""Tests for `sensormapgeo` package."""


import os
from unittest import TestCase

import numpy as np
from gdalnumeric import LoadFile
from py_tools_ds.geo.coord_calc import corner_coord_to_minmax, get_corner_coordinates

from sensormapgeo import __path__
from sensormapgeo import SensorMapGeometryTransformer, SensorMapGeometryTransformer3D


tests_path = os.path.abspath(os.path.join(__path__[0], "..", "tests"))
rsp_algs = ['nearest', 'bilinear', 'gauss']
42
mp_algs = ['bands', 'tiles']
Daniel Scheffler's avatar
Daniel Scheffler committed
43
44
45
46
47
48
49
50


class Test_SensorMapGeometryTransformer(TestCase):
    def setUp(self):
        self.dem_map_geo = LoadFile(os.path.join(tests_path, 'data', 'dem_map_geo.bsq'))
        self.dem_sensor_geo = LoadFile(os.path.join(tests_path, 'data', 'dem_sensor_geo.bsq'))
        self.lons = LoadFile(os.path.join(tests_path, 'data', 'lons_full_vnir.bsq'))
        self.lats = LoadFile(os.path.join(tests_path, 'data', 'lats_full_vnir.bsq'))
51
        self.dem_area_extent_coarse_subset_utm = (622613.864409047,  # LL_x
Daniel Scheffler's avatar
Daniel Scheffler committed
52
53
                                                  5254111.40255343,  # LL_x
                                                  660473.864409047,  # LL_x
54
                                                  5269351.40255343)  # UR_y
Daniel Scheffler's avatar
Daniel Scheffler committed
55

56
        self.expected_dem_area_extent_lonlat = (10.685733901515151,  # LL_x
Daniel Scheffler's avatar
Daniel Scheffler committed
57
58
                                                47.44113415492957,  # LL_y
                                                11.073066098484848,  # UR_x
59
                                                47.54576584507042)  # UR_y
Daniel Scheffler's avatar
Daniel Scheffler committed
60

61
        self.expected_dem_area_extent_utm = (626938.928052,  # LL_x
Daniel Scheffler's avatar
Daniel Scheffler committed
62
63
                                             5256253.56579,  # LL_y
                                             656188.928052,  # UR_x
64
                                             5267203.56579)  # UR_y
Daniel Scheffler's avatar
Daniel Scheffler committed
65

Daniel Scheffler's avatar
Daniel Scheffler committed
66
67
68
69
70
        self.expected_dem_area_extent_utm_ongrid = (626910,  # LL_x
                                                    5256240,  # LL_y
                                                    656190,  # UR_x
                                                    5267220)  # UR_y

Daniel Scheffler's avatar
Daniel Scheffler committed
71
72
73
74
75
76
    def test_to_sensor_geometry(self):
        for rsp_alg in rsp_algs:
            SMGT = SensorMapGeometryTransformer(lons=self.lons,
                                                lats=self.lats,
                                                resamp_alg=rsp_alg,
                                                radius_of_influence=30 if rsp_alg != 'bilinear' else 45)
77
78
79
80
81
82
            dem_sensor_geo = SMGT.to_sensor_geometry(self.dem_map_geo,
                                                     src_prj=32632, src_extent=self.dem_area_extent_coarse_subset_utm)
            self.assertIsInstance(dem_sensor_geo, np.ndarray)
            self.assertFalse(np.array_equal(np.unique(dem_sensor_geo), np.array([0])))
            self.assertEqual(dem_sensor_geo.shape, (150, 1000))
            self.assertEqual(self.dem_map_geo.dtype, dem_sensor_geo.dtype)
Daniel Scheffler's avatar
Daniel Scheffler committed
83
84
85
86
87
88

    def test_to_sensor_geometry_3DInput(self):
        for rsp_alg in rsp_algs:
            SMGT = SensorMapGeometryTransformer(lons=self.lons,
                                                lats=self.lats,
                                                resamp_alg=rsp_alg)
89
90
91
92
93
94
95
            dem_sensor_geo = SMGT.to_sensor_geometry(np.dstack([self.dem_map_geo] * 2),
                                                     src_prj=32632, src_extent=self.dem_area_extent_coarse_subset_utm)
            self.assertIsInstance(dem_sensor_geo, np.ndarray)
            self.assertFalse(np.array_equal(np.unique(dem_sensor_geo), np.array([0])))
            self.assertEqual(dem_sensor_geo.shape, (150, 1000, 2))
            self.assertTrue(np.array_equal(dem_sensor_geo[:, :, 0], dem_sensor_geo[:, :, 1]))
            self.assertEqual(self.dem_map_geo.dtype, dem_sensor_geo.dtype)
Daniel Scheffler's avatar
Daniel Scheffler committed
96
97
98
99
100
101
102
103
104

    def test_to_map_geometry_lonlat(self):
        for rsp_alg in rsp_algs:
            SMGT = SensorMapGeometryTransformer(lons=self.lons,
                                                lats=self.lats,
                                                resamp_alg=rsp_alg)

            # to Lon/Lat
            dem_map_geo, dem_gt, dem_prj = SMGT.to_map_geometry(self.dem_sensor_geo, tgt_prj=4326)
105
106
107
108
109

            # from geoarray import GeoArray
            # GeoArray(dem_map_geo, dem_gt, dem_prj)\
            #     .save(os.path.join(tests_path, 'test_output', 'resampled_pyresample_ll.bsq'))

Daniel Scheffler's avatar
Daniel Scheffler committed
110
            self.assertIsInstance(dem_map_geo, np.ndarray)
111
112
            self.assertEqual(dem_map_geo.shape, (SMGT.area_definition.height,
                                                 SMGT.area_definition.width))
Daniel Scheffler's avatar
Daniel Scheffler committed
113
114
115
116
117
118
            xmin, xmax, ymin, ymax = corner_coord_to_minmax(get_corner_coordinates(gt=dem_gt,
                                                                                   cols=dem_map_geo.shape[1],
                                                                                   rows=dem_map_geo.shape[0]))
            self.assertTrue(False not in np.isclose(np.array([xmin, ymin, xmax, ymax]),
                                                    np.array(self.expected_dem_area_extent_lonlat)))
            self.assertFalse(np.array_equal(np.unique(dem_map_geo), np.array([0])))
Daniel Scheffler's avatar
Daniel Scheffler committed
119
120
121
            self.assertTrue(np.isclose(np.mean(dem_map_geo[dem_map_geo != 0]),
                                       np.mean(self.dem_sensor_geo),
                                       rtol=0.01))
122
            self.assertEqual(self.dem_sensor_geo.dtype, dem_map_geo.dtype)
Daniel Scheffler's avatar
Daniel Scheffler committed
123
124
125
126
127
128
129
130

            with self.assertRaises(ValueError):
                SMGT.to_map_geometry(self.dem_sensor_geo[:10, :10], tgt_prj=4326)  # must have the shape of lons/lats

    def test_to_map_geometry_utm(self):
        for rsp_alg in rsp_algs:
            SMGT = SensorMapGeometryTransformer(lons=self.lons,
                                                lats=self.lats,
131
132
133
134
135
                                                resamp_alg=rsp_alg,
                                                # neighbours=8,
                                                # radius_of_influence=45,
                                                # epsilon=0
                                                )
Daniel Scheffler's avatar
Daniel Scheffler committed
136
137

            # to UTM32
Daniel Scheffler's avatar
Daniel Scheffler committed
138
139
140
141
142
143
            # dem_map_geo, dem_gt, dem_prj = SMGT.to_map_geometry(self.dem_sensor_geo, tgt_prj=32632, tgt_res=(30, 30))
            dem_map_geo, dem_gt, dem_prj = SMGT.to_map_geometry(self.dem_sensor_geo, tgt_prj=32632,
                                                                tgt_res=(30, 30),
                                                                # tgt_extent=self.expected_dem_area_extent_utm,
                                                                tgt_coordgrid=((0, 30), (0, 30))
                                                                )
144
            # from geoarray import GeoArray
145
146
            # GeoArray(dem_map_geo, dem_gt, dem_prj).save(os.path.join(tests_path, 'test_output',
            #                                                          'resampled_pyresample_bilinear_16n.bsq'))
Daniel Scheffler's avatar
Daniel Scheffler committed
147

Daniel Scheffler's avatar
Daniel Scheffler committed
148
            self.assertIsInstance(dem_map_geo, np.ndarray)
Daniel Scheffler's avatar
Daniel Scheffler committed
149
            self.assertEqual(dem_map_geo.shape, (366, 976))
Daniel Scheffler's avatar
Daniel Scheffler committed
150
151
152
153
            xmin, xmax, ymin, ymax = corner_coord_to_minmax(get_corner_coordinates(gt=dem_gt,
                                                                                   cols=dem_map_geo.shape[1],
                                                                                   rows=dem_map_geo.shape[0]))
            self.assertTrue(False not in np.isclose(np.array([xmin, ymin, xmax, ymax]),
Daniel Scheffler's avatar
Daniel Scheffler committed
154
                                                    np.array(self.expected_dem_area_extent_utm_ongrid)))
Daniel Scheffler's avatar
Daniel Scheffler committed
155
            self.assertFalse(np.array_equal(np.unique(dem_map_geo), np.array([0])))
Daniel Scheffler's avatar
Daniel Scheffler committed
156
157
158
            self.assertTrue(np.isclose(np.mean(dem_map_geo[dem_map_geo != 0]),
                                       np.mean(self.dem_sensor_geo),
                                       rtol=0.01))
159
            self.assertEqual(self.dem_sensor_geo.dtype, dem_map_geo.dtype)
Daniel Scheffler's avatar
Daniel Scheffler committed
160
161
162
163
164
165
166
167
168
169
170


class Test_SensorMapGeometryTransformer3D(TestCase):
    def setUp(self):
        dem_map_geo = LoadFile(os.path.join(tests_path, 'data', 'dem_map_geo.bsq'))
        dem_sensor_geo = LoadFile(os.path.join(tests_path, 'data', 'dem_sensor_geo.bsq'))
        lons = LoadFile(os.path.join(tests_path, 'data', 'lons_full_vnir.bsq'))
        lats = LoadFile(os.path.join(tests_path, 'data', 'lats_full_vnir.bsq'))

        self.data_map_geo_3D = np.dstack([dem_map_geo, dem_map_geo])
        self.data_sensor_geo_3D = np.dstack([dem_sensor_geo, dem_sensor_geo])
Daniel Scheffler's avatar
Daniel Scheffler committed
171
172
        self.lons_3D = np.dstack([lons, lons + .002])  # assume slightly different coordinates in both bands
        self.lats_3D = np.dstack([lats, lats + .002])
Daniel Scheffler's avatar
Daniel Scheffler committed
173

174
        self.dem_area_extent_coarse_subset_utm = (622613.864409047,  # LL_x
Daniel Scheffler's avatar
Daniel Scheffler committed
175
176
                                                  5254111.40255343,  # LL_x
                                                  660473.864409047,  # LL_x
177
                                                  5269351.40255343)  # UR_y
Daniel Scheffler's avatar
Daniel Scheffler committed
178

Daniel Scheffler's avatar
Daniel Scheffler committed
179
        # this differs from the 2D version because the geolayer in the second band has slightly different coordinates
180
        self.expected_dem_area_extent_lonlat = (10.685733901515151,  # LL_x
Daniel Scheffler's avatar
Daniel Scheffler committed
181
                                                47.44113415492957,  # LL_y
Daniel Scheffler's avatar
Daniel Scheffler committed
182
183
                                                11.075064739115845,  # UR_x
                                                47.54772559829233)  # UR_y
Daniel Scheffler's avatar
Daniel Scheffler committed
184

185
        self.expected_dem_area_extent_utm = (626938.928052,  # LL_x
Daniel Scheffler's avatar
Daniel Scheffler committed
186
187
                                             5256253.56579,  # LL_y
                                             656188.928052,  # UR_x
188
                                             5267203.56579)  # UR_y
Daniel Scheffler's avatar
Daniel Scheffler committed
189

Daniel Scheffler's avatar
Daniel Scheffler committed
190
191
192
193
194
195
        # this differs from the 2D version because the geolayer in the second band has slightly different coordinates
        self.expected_dem_area_extent_utm_ongrid = (626910,  # LL_x
                                                    5256240,  # LL_y
                                                    656340,  # UR_x
                                                    5267430)  # UR_y

Daniel Scheffler's avatar
Daniel Scheffler committed
196
197
    def test_to_map_geometry_lonlat_3D_geolayer(self):
        for rsp_alg in rsp_algs:
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
            for mp_alg in mp_algs:
                SMGT = SensorMapGeometryTransformer3D(lons=self.lons_3D,
                                                      lats=self.lats_3D,
                                                      # resamp_alg='nearest',
                                                      resamp_alg=rsp_alg,
                                                      mp_alg=mp_alg
                                                      )

                # to Lon/Lat
                data_mapgeo_3D, dem_gt, dem_prj = SMGT.to_map_geometry(self.data_sensor_geo_3D, tgt_prj=4326)

                # from geoarray import GeoArray
                # GeoArray(data_mapgeo_3D, dem_gt, dem_prj)\
                #     .save(os.path.join(tests_path, 'test_output', 'resampled_3D_02_ll.bsq'))

                self.assertIsInstance(data_mapgeo_3D, np.ndarray)
                # only validate number of bands (height and width are validated in 2D version
                #   fixed numbers may fail here due to float uncertainty errors
                self.assertGreater(np.dot(*data_mapgeo_3D.shape[:2]), np.dot(*self.data_sensor_geo_3D.shape[:2]))
                self.assertEqual(data_mapgeo_3D.shape[2], 2)
                xmin, xmax, ymin, ymax = corner_coord_to_minmax(get_corner_coordinates(gt=dem_gt,
                                                                                       cols=data_mapgeo_3D.shape[1],
                                                                                       rows=data_mapgeo_3D.shape[0]))
                self.assertTrue(False not in np.isclose(np.array([xmin, ymin, xmax, ymax]),
                                                        np.array(self.expected_dem_area_extent_lonlat)))

                self.assertTrue(np.isclose(np.mean(data_mapgeo_3D[data_mapgeo_3D != 0]),
                                           np.mean(self.data_sensor_geo_3D),
                                           rtol=0.01))
                self.assertEqual(self.data_sensor_geo_3D.dtype, data_mapgeo_3D.dtype)

    # def test_to_map_geometry_lonlat_3D_geolayer_bands(self):
    #     SMGT = SensorMapGeometryTransformer3D(lons=self.lons_3D,
    #                                           lats=self.lats_3D,
    #                                           # resamp_alg='nearest',
    #                                           resamp_alg='gauss',
    #                                           mp_alg='bands'
    #                                           )
    #
    #     # to Lon/Lat
    #     data_mapgeo_3D, dem_gt, dem_prj = SMGT.to_map_geometry(self.data_sensor_geo_3D, tgt_prj=4326)
    #
    #     # from geoarray import GeoArray
    #     # GeoArray(data_mapgeo_3D, dem_gt, dem_prj)\
    #     #     .save(os.path.join(tests_path, 'test_output', 'resampled_3D_02_ll.bsq'))
    #
    #     self.assertIsInstance(data_mapgeo_3D, np.ndarray)
    #     # only validate number of bands (height and width are validated in 2D version
    #     #   fixed numbers may fail here due to float uncertainty errors
    #     self.assertGreater(np.dot(*data_mapgeo_3D.shape[:2]), np.dot(*self.data_sensor_geo_3D.shape[:2]))
    #     self.assertEqual(data_mapgeo_3D.shape[2], 2)
    #     xmin, xmax, ymin, ymax = corner_coord_to_minmax(get_corner_coordinates(gt=dem_gt,
    #                                                                            cols=data_mapgeo_3D.shape[1],
    #                                                                            rows=data_mapgeo_3D.shape[0]))
    #     self.assertTrue(False not in np.isclose(np.array([xmin, ymin, xmax, ymax]),
    #                                             np.array(self.expected_dem_area_extent_lonlat)))
    #
    #     self.assertTrue(np.isclose(np.mean(data_mapgeo_3D[data_mapgeo_3D != 0]),
    #                                np.mean(self.data_sensor_geo_3D),
    #                                rtol=0.01))
    #     self.assertEqual(self.data_sensor_geo_3D.dtype, data_mapgeo_3D.dtype)
Daniel Scheffler's avatar
Daniel Scheffler committed
259
260
261

    def test_to_map_geometry_utm_3D_geolayer(self):
        for rsp_alg in rsp_algs:
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
            for mp_alg in mp_algs:
                SMGT = SensorMapGeometryTransformer3D(lons=self.lons_3D,
                                                      lats=self.lats_3D,
                                                      # resamp_alg='nearest',
                                                      resamp_alg=rsp_alg,
                                                      mp_alg=mp_alg
                                                      )

                # to Lon/Lat
                data_mapgeo_3D, dem_gt, dem_prj = SMGT.to_map_geometry(self.data_sensor_geo_3D,
                                                                       tgt_prj=32632,
                                                                       tgt_res=(30, 30),
                                                                       # tgt_extent=self.expected_dem_area_extent_utm,
                                                                       tgt_coordgrid=((0, 30), (0, 30))
                                                                       )
                # from geoarray import GeoArray
                # GeoArray(data_mapgeo_3D, dem_gt, dem_prj)\
                #     .save(os.path.join(tests_path, 'test_output', 'resampled_3D_02_pyresample.bsq'))

                self.assertIsInstance(data_mapgeo_3D, np.ndarray)
                # only validate number of bands (height and width are validated in 2D version
                #   fixed numbers may fail here due to float uncertainty errors
                self.assertGreater(np.dot(*data_mapgeo_3D.shape[:2]), np.dot(*self.data_sensor_geo_3D.shape[:2]))
                self.assertEqual(data_mapgeo_3D.shape[2], 2)
                xmin, xmax, ymin, ymax = corner_coord_to_minmax(get_corner_coordinates(gt=dem_gt,
                                                                                       cols=data_mapgeo_3D.shape[1],
                                                                                       rows=data_mapgeo_3D.shape[0]))
                self.assertTrue(False not in np.isclose(np.array([xmin, ymin, xmax, ymax]),
                                                        np.array(self.expected_dem_area_extent_utm_ongrid)))

                self.assertTrue(np.isclose(np.mean(data_mapgeo_3D[data_mapgeo_3D != 0]),
                                           np.mean(self.data_sensor_geo_3D),
                                           rtol=0.01))
                self.assertEqual(self.data_sensor_geo_3D.dtype, data_mapgeo_3D.dtype)
Daniel Scheffler's avatar
Daniel Scheffler committed
296

Daniel Scheffler's avatar
Daniel Scheffler committed
297
298
    def test_to_sensor_geometry(self):
        for rsp_alg in rsp_algs:
299
300
301
302
303
304
305
306
307
308
309
            for mp_alg in mp_algs:
                SMGT = SensorMapGeometryTransformer3D(lons=self.lons_3D,
                                                      lats=self.lats_3D,
                                                      resamp_alg=rsp_alg,
                                                      mp_alg=mp_alg
                                                      )
                dem_sensors_geo = SMGT.to_sensor_geometry(self.data_map_geo_3D,
                                                          src_prj=32632, src_extent=self.dem_area_extent_coarse_subset_utm)
                self.assertIsInstance(dem_sensors_geo, np.ndarray)
                self.assertEqual(dem_sensors_geo.shape, (150, 1000,  2))
                self.assertEqual(self.data_map_geo_3D.dtype, dem_sensors_geo.dtype)