metadata.py 10.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
# -*- coding: utf-8 -*-

from datetime import datetime
from xml.etree import ElementTree
import logging
import numpy as np
from scipy.interpolate import interp2d

L1B_product_props = dict(
    xml_detector_label=dict(
        VNIR='detector1',
Daniel Scheffler's avatar
Daniel Scheffler committed
12
        SWIR='detector2'),
13
14
    fn_detector_suffix=dict(
        VNIR='D1',
Daniel Scheffler's avatar
Daniel Scheffler committed
15
        SWIR='D2')
16
17
18
19
20
21
22
23
24
)


##############
# BASE CLASSES
##############


class _EnMAP_Metadata_Detector_ImGeo(object):
25
26
27
28
29
30
31
32
    """Base class for all EnMAP metadata associated with a single EnMAP detector in image geometry.

    NOTE:
        - All metadata that have VNIR and SWIR detector in image geometry in common should be included here.
        - All classes representing VNIR or SWIR specific metadata should inherit from _EnMAP_Metadata_Detector_ImGeo.

    """
    def __init__(self, detector_name: str, logger: logging.Logger=None):
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
        """

        :param detector_name: Name of the detector (VNIR or SWIR)
        :param logger:
        """

        self.detector_name = detector_name  # type: str
        self.logger = logger or logging.getLogger()

        self.fwhm = None  # type: np.ndarray  # Full width half maximmum for each EnMAP band
        self.wvl_center = None  # type: np.ndarray  # Center wavelengths for each EnMAP band
        self.nwvl = None  # type: int  # Number of wave bands
        self.nrows = None  # type: int  # number of rows
        self.ncols = None  # type: int  # number of columns
        self.smile_coef = None  # type: np.ndarray  # smile coefficients needed for smile computation
        self.nsmile_coef = None  # type: int  # number of smile coefficients
        self.smile = None  # type: np.ndarray  # smile for each EnMAP image column
        self.l_min = None  # type: np.ndarray
        self.l_max = None  # type: np.ndarray
        self.geom_view_zenith = None  # type: float  # viewing zenith angle
        self.geom_view_azimuth = None  # type: float  # viewing azimuth angle
        self.geom_sun_zenith = None  # type: float  # sun zenith angle
        self.geom_sun_azimuth = None  # type: float  # sun azimuth angle
        self.mu_sun = None  # type: float  # earth-sun distance # TODO doc correct?
Daniel Scheffler's avatar
Daniel Scheffler committed
57
58
        self.lat_UL_UR_LL_LR = None  # type:  list  # latitude coordinates for UL, UR, LL, LR
        self.lon_UL_UR_LL_LR = None  # type:  list  # longitude coordinates for UL, UR, LL, LR
59
60
61
62
63
64
        self.lats = None  # type: np.ndarray  # 2D array of latitude coordinates according to given lon/lat sampling
        self.lons = None  # type: np.ndarray  # 2D array of longitude coordinates according to given lon/lat sampling
        self.unit = ''  # type: str  # radiometric unit of pixel values
        self.snr = None  # type: np.ndarray  # Signal to noise ratio as computed from radiance data

    def _read_metadata(self, path_xml, detector_label_xml, lon_lat_smpl=(15, 15), nsmile_coef=4):
65
66
67
68
69
70
71
        """Read the metadata of a specific EnMAP detector in image geometry.

        :param path_xml:  file path of the metadata XML file
        :param detector_label_xml:
        :param lon_lat_smpl:  number if sampling points in lon, lat fields
        :param nsmile_coef:  number of polynomial coefficients for smile
        """
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
        xml = ElementTree.parse(path_xml).getroot()
        lbl = detector_label_xml
        self.logger.info("Load data for: %s" % lbl)

        self.fwhm = np.array(xml.findall("%s/fwhm" % lbl)[0].text.replace("\n", "").split(), dtype=np.float)
        self.wvl_center = np.array(
            xml.findall("%s/centre_wavelength" % lbl)[0].text.replace("\n", "").split(), dtype=np.float)
        self.nwvl = len(self.wvl_center)
        self.nrows = np.int(xml.findall("%s/rows" % lbl)[0].text)
        self.ncols = np.int(xml.findall("%s/columns" % lbl)[0].text)
        self.smile_coef = np.array(xml.findall("%s/smile" % lbl)[0].text.replace("\n", "").split(), dtype=np.float) \
                            .reshape((-1, nsmile_coef + 1))[:, 1:]
        self.nsmile_coef = nsmile_coef
        self.smile = self.calc_smile()
        self.l_min = np.array(xml.findall("%s/L_min" % lbl)[0].text.split(), dtype=np.float)
        self.l_max = np.array(xml.findall("%s/L_max" % lbl)[0].text.split(), dtype=np.float)
        self.geom_view_zenith = np.float(
            xml.findall("%s/observation_geometry/zenith_angle" % lbl)[0].text.split()[0])
        self.geom_view_azimuth = np.float(
Daniel Scheffler's avatar
Daniel Scheffler committed
91
            xml.findall("%s/observation_geometry/azimuth_angle" % lbl)[0].text.split()[0])
92
93
        self.geom_sun_zenith = np.float(
            xml.findall("%s/illumination_geometry/zenith_angle" % lbl)[0].text.split()[0])
Daniel Scheffler's avatar
Daniel Scheffler committed
94
95
        self.geom_sun_azimuth = np.float(
            xml.findall("%s/illumination_geometry/azimuth_angle" % lbl)[0].text.split()[0])
96
97
        self.mu_sun = np.cos(np.deg2rad(self.geom_sun_zenith))
        self.lat_UL_UR_LL_LR = \
Daniel Scheffler's avatar
Daniel Scheffler committed
98
99
            [float(xml.findall("%s/geometry/bounding_box/%s_northing" % (lbl, corner))[0].text.split()[0])
             for corner in ("UL", "UR", "LL", "LR")]
100
        self.lon_UL_UR_LL_LR = \
Daniel Scheffler's avatar
Daniel Scheffler committed
101
102
            [float(xml.findall("%s/geometry/bounding_box/%s_easting" % (lbl, corner))[0].text.split()[0])
             for corner in ("UL", "UR", "LL", "LR")]
103
104
        self.lats = self.interpolate_corners(*self.lat_UL_UR_LL_LR, *lon_lat_smpl)
        self.lons = self.interpolate_corners(*self.lon_UL_UR_LL_LR, *lon_lat_smpl)
105
106
        self.unit = 'none'  # '" ".join(xml.findall("%s/radiance_unit" % lbl)[0].text.split())
        self.unitcode = 'DN'
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149

    def calc_smile(self):
        """Compute smile for each EnMAP column.
        The sum in (1) is expressed as inner product of two arrays with inner dimension as the polynomial smile
        coefficients shape = (ncols, nsmile_coef) of polynomial x

        :return:
        """
        # smile(icol, iwvl) = sum_(p=0)^(nsmile_coef-1) smile_coef[iwvl, p] * icol**p (1)
        return np.inner(
            np.array([[icol ** p for p in np.arange(self.nsmile_coef)] for icol in np.arange(self.ncols)]),
            self.smile_coef  # shape = (nwvl, nsmile_coef)
        )  # shape = (ncols, nwvl)

    def calc_snr(self, data: np.ndarray):
        """Compute EnMAP snr from radiance data.

        :param data: Numpy array with radiance for scene
        """
        self.logger.info("Compute snr for: %s" % self.detector_name)
        self.logger.warning("SNR model missing -> const. value of 500 is returned")
        return 500 * np.ones(data.shape, dtype=np.float)

    @staticmethod
    def interpolate_corners(ul: float, ur: float, ll: float, lr: float, nx: int, ny: int):
        """Compute interpolated field from corner values of a scalar field given at: ul, ur, ll, lr.

        :param nx, ny: final shape
        """
        ff = interp2d(x=[0, 1], y=[0, 1], z=[[ul, ur], [ll, lr]], kind='linear')
        rr = np.zeros((nx, ny), dtype=np.float)
        for i, x in enumerate(np.linspace(0, 1, nx)):
            for j, y in enumerate(np.linspace(0, 1, ny)):
                rr[i, j] = ff(x, y)
        return rr


########################################################
# EnPT metadata objects for EnMAP data in image geometry
########################################################


class EnMAP_Metadata_ImGeo(object):
150
151
152
153
154
155
156
157
158
    """EnMAP Metadata class holding the metadata of the complete EnMAP product in image geometry
    including VNIR and SWIR detector.

    Attributes:
        - logger(logging.Logger):  None or logging instance
        - observation_datetime(datetime.datetime):  datetime of observation time (currently missing in metadata)
        - vnir(EnMAP_Metadata_VNIR_ImGeo)
        - swir(EnMAP_Metadata_SWIR_ImGeo)
    """
159
160
    def __init__(self, path_metaxml, logger=None):
        self.logger = logger or logging.getLogger()
Daniel Scheffler's avatar
Daniel Scheffler committed
161
        self._path_xml = path_metaxml
162

163
        # defaults
164
165
166
167
        self.observation_datetime = None  # type: datetime  # Date and Time of image observation
        self.vnir = None  # type: EnMAP_Metadata_VNIR_ImGeo # metadata of VNIR only
        self.swir = None  # type: EnMAP_Metadata_SWIR_ImGeo # metadata of SWIR only

Daniel Scheffler's avatar
Daniel Scheffler committed
168
    def read_common_meta(self, observation_time: datetime=None):
169
170
171
        # FIXME observation time is currently missing in the XML
        self.observation_datetime = observation_time

Daniel Scheffler's avatar
Daniel Scheffler committed
172
173
174
    def read_metadata(self, observation_time: datetime=None, lon_lat_smpl=(15, 15), nsmile_coef=4):
        self.read_common_meta(observation_time)
        self.vnir = EnMAP_Metadata_VNIR_ImGeo(self._path_xml)
175
        self.vnir.read_metadata(lon_lat_smpl=lon_lat_smpl, nsmile_coef=nsmile_coef)
Daniel Scheffler's avatar
Daniel Scheffler committed
176
        self.swir = EnMAP_Metadata_SWIR_ImGeo(self._path_xml)
177
178
179
180
        self.swir.read_metadata(lon_lat_smpl=lon_lat_smpl, nsmile_coef=nsmile_coef)


class EnMAP_Metadata_VNIR_ImGeo(_EnMAP_Metadata_Detector_ImGeo):
181
182
183
184
185
    """EnMAP Metadata class holding the metadata of the VNIR detector in image geometry.

    NOTE:
        - inherits all attributes from base class _EnMAP_Metadata_Detector_ImGeo.
    """
186
187
188
    def __init__(self, path_metaxml, logger=None):
        # get all attributes from base class '_EnMAP_Metadata_Detector_ImGeo'
        super(EnMAP_Metadata_VNIR_ImGeo, self).__init__('VNIR', logger=logger)
Daniel Scheffler's avatar
Daniel Scheffler committed
189
        self._path_xml = path_metaxml
190
191
192
193
        self.detector_label = L1B_product_props['xml_detector_label']['VNIR']

    def read_metadata(self, lon_lat_smpl=(15, 15), nsmile_coef=4):
        super(EnMAP_Metadata_VNIR_ImGeo, self)\
Daniel Scheffler's avatar
Daniel Scheffler committed
194
            ._read_metadata(self._path_xml, self.detector_label, lon_lat_smpl=lon_lat_smpl, nsmile_coef=nsmile_coef)
195
196
197


class EnMAP_Metadata_SWIR_ImGeo(_EnMAP_Metadata_Detector_ImGeo):
198
199
200
201
202
    """EnMAP Metadata class holding the metadata of the SWIR detector in image geometry.

    NOTE:
        - inherits all attributes from base class _EnMAP_Metadata_Detector_ImGeo.
    """
203
204
    def __init__(self, path_metaxml, logger=None):
        # get all attributes from base class '_EnMAP_Metadata_Detector_ImGeo'
Daniel Scheffler's avatar
Daniel Scheffler committed
205
206
        super(EnMAP_Metadata_SWIR_ImGeo, self).__init__('SWIR', logger=logger)
        self._path_xml = path_metaxml
207
208
209
        self.detector_label = L1B_product_props['xml_detector_label']['SWIR']

    def read_metadata(self, lon_lat_smpl=(15, 15), nsmile_coef=4):
210
211
212
213
214
215
        """Read the metadata of the SWIR detector in image geometry.

        :param lon_lat_smpl: number if sampling points in lon, lat fields
        :param nsmile_coef: number of polynomial coefficients for smile
        :return:
        """
216
        super(EnMAP_Metadata_SWIR_ImGeo, self)\
Daniel Scheffler's avatar
Daniel Scheffler committed
217
            ._read_metadata(self._path_xml, self.detector_label, lon_lat_smpl=lon_lat_smpl, nsmile_coef=nsmile_coef)