test_sensormapgeo.py 11.6 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
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
#!/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']


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'))
        self.dem_area_extent_coarse_subset_utm = [622613.864409047,  # LL_x
                                                  5254111.40255343,  # LL_x
                                                  660473.864409047,  # LL_x
                                                  5269351.40255343]  # UR_y

        self.expected_dem_area_extent_lonlat = [10.685733901515151,  # LL_x
                                                47.44113415492957,  # LL_y
                                                11.073066098484848,  # UR_x
                                                47.54576584507042]  # UR_y

        self.expected_dem_area_extent_utm = [626938.928052,  # LL_x
                                             5256253.56579,  # LL_y
                                             656188.928052,  # UR_x
                                             5267203.56579]  # UR_y

    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)
71
72
73
74
75
76
            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
77
78
79
80
81
82

    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)
83
84
85
86
87
88
89
            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
90
91
92
93
94
95
96
97
98
99

    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)
            self.assertIsInstance(dem_map_geo, np.ndarray)
100
101
            self.assertEqual(dem_map_geo.shape, (SMGT.area_definition.height,
                                                 SMGT.area_definition.width))
Daniel Scheffler's avatar
Daniel Scheffler committed
102
103
104
105
106
107
            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
108
109
110
            self.assertTrue(np.isclose(np.mean(dem_map_geo[dem_map_geo != 0]),
                                       np.mean(self.dem_sensor_geo),
                                       rtol=0.01))
111
            self.assertEqual(self.dem_sensor_geo.dtype, dem_map_geo.dtype)
Daniel Scheffler's avatar
Daniel Scheffler committed
112
113
114
115
116
117
118
119
120
121
122
123
124

            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,
                                                resamp_alg=rsp_alg)

            # to UTM32
            dem_map_geo, dem_gt, dem_prj = SMGT.to_map_geometry(self.dem_sensor_geo, tgt_prj=32632, tgt_res=(30, 30))
            self.assertIsInstance(dem_map_geo, np.ndarray)
125
            self.assertEqual(dem_map_geo.shape, (365, 975))
Daniel Scheffler's avatar
Daniel Scheffler committed
126
127
128
129
130
131
            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_utm)))
            self.assertFalse(np.array_equal(np.unique(dem_map_geo), np.array([0])))
Daniel Scheffler's avatar
Daniel Scheffler committed
132
133
134
            self.assertTrue(np.isclose(np.mean(dem_map_geo[dem_map_geo != 0]),
                                       np.mean(self.dem_sensor_geo),
                                       rtol=0.01))
135
            self.assertEqual(self.dem_sensor_geo.dtype, dem_map_geo.dtype)
Daniel Scheffler's avatar
Daniel Scheffler committed
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177


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])
        self.lons_3D = np.dstack([lons, lons])  # TODO use different lons per band here
        self.lats_3D = np.dstack([lats, lats])  # TODO use different lats per band here

        self.dem_area_extent_coarse_subset_utm = [622613.864409047,  # LL_x
                                                  5254111.40255343,  # LL_x
                                                  660473.864409047,  # LL_x
                                                  5269351.40255343]  # UR_y

        self.expected_dem_area_extent_lonlat = [10.685733901515151,  # LL_x
                                                47.44113415492957,  # LL_y
                                                11.073066098484848,  # UR_x
                                                47.54576584507042]  # UR_y

        self.expected_dem_area_extent_utm = [626938.928052,  # LL_x
                                             5256253.56579,  # LL_y
                                             656188.928052,  # UR_x
                                             5267203.56579]  # UR_y

    def test_to_map_geometry_lonlat_3D_geolayer(self):
        for rsp_alg in rsp_algs:
            SMGT = SensorMapGeometryTransformer3D(lons=self.lons_3D,
                                                  lats=self.lats_3D,
                                                  # resamp_alg='nearest',
                                                  resamp_alg=rsp_alg,
                                                  )

            # to Lon/Lat
            data_mapgeo_3D, dem_gt, dem_prj = SMGT.to_map_geometry(self.data_sensor_geo_3D, tgt_prj=4326)
            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
Daniel Scheffler's avatar
Daniel Scheffler committed
178
179
            self.assertGreater(data_mapgeo_3D.shape[0], self.data_sensor_geo_3D.shape[0])
            self.assertGreater(data_mapgeo_3D.shape[1], self.data_sensor_geo_3D.shape[1])
180
            self.assertEqual(data_mapgeo_3D.shape[2], 2)
Daniel Scheffler's avatar
Daniel Scheffler committed
181
182
183
184
185
186
            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)))

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

Daniel Scheffler's avatar
Daniel Scheffler committed
192
193
194
195
196
197
198
199
200
    def test_to_sensor_geometry(self):
        for rsp_alg in rsp_algs:
            SMGT = SensorMapGeometryTransformer3D(lons=self.lons_3D,
                                                  lats=self.lats_3D,
                                                  resamp_alg=rsp_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)
201
202
            self.assertEqual(dem_sensors_geo.shape, (150, 1000,  2))
            self.assertEqual(self.data_map_geo_3D.dtype, dem_sensors_geo.dtype)