get_reflectance_lut() - Code Metrics - Inspection of "Fix Rayleigh corrector to work with dask" - pytroll/pyspectral - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Pull Request — develop (#33)

unknown

created 2018-04-08 21:27 UTC

get_reflectance_lut() B

↳ Parent: Project

Complexity

Conditions

Size

Total Lines

Duplication

Lines	0
Ratio	0 %

Importance

Changes	4
Bugs	0	Features	0

Metric	Value
cc	2
dl	0
loc	32
rs	8.8571
c	4
b	0
f	0

#!/usr/bin/env python
# -*- coding: utf-8 -*-

# Copyright (c) 2016-2018 Pytroll

# Author(s):

#   Adam.Dybbroe <[email protected]>
#   Martin Raspaud <[email protected]>

# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU 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 General Public License for more details.

# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

"""Atmospheric correction of shortwave imager bands in the wavelength range 400
to 800 nm

"""

import os
import time
import logging
from six import integer_types

import h5py
import numpy as np

try:
    from dask.array import where, zeros, map_blocks, from_array, clip
    HAVE_DASK = True
except ImportError:
    from numpy import where, zeros, clip
    map_blocks = None
    from_array = None
    HAVE_DASK = False

from geotiepoints.multilinear import MultilinearInterpolator
from pyspectral.rsr_reader import RelativeSpectralResponse
from pyspectral.utils import (INSTRUMENTS, RAYLEIGH_LUT_DIRS,
                              AEROSOL_TYPES, ATMOSPHERES,
                              BANDNAMES,
                              download_luts, get_central_wave,
                              get_bandname_from_wavelength)

LOG = logging.getLogger(__name__)


class BandFrequencyOutOfRange(ValueError):

    """Exception when the band frequency is out of the visible range"""

    pass


class Rayleigh(object):

    """Container for the atmospheric correction of satellite imager bands.

    This class removes background contributions of Rayleigh scattering of
    molecules and Mie scattering and absorption by aerosols.

    """

    def __init__(self, platform_name, sensor, **kwargs):
        """Initialize class and determine LUT to use."""
        self.platform_name = platform_name
        self.sensor = sensor
        self.coeff_filename = None

        atm_type = kwargs.get('atmosphere', 'us-standard')
        if atm_type not in ATMOSPHERES:
            raise AttributeError('Atmosphere type not supported! ' +
                                 'Need to be one of {}'.format(str(ATMOSPHERES)))

        aerosol_type = kwargs.get('aerosol_type', 'marine_clean_aerosol')

        if aerosol_type not in AEROSOL_TYPES:
            raise AttributeError('Aerosol type not supported! ' +
                                 'Need to be one of {0}'.format(str(AEROSOL_TYPES)))

        rayleigh_dir = RAYLEIGH_LUT_DIRS[aerosol_type]

        if atm_type not in ATMOSPHERES.keys():
            LOG.error("Atmosphere type %s not supported", atm_type)

        LOG.info("Atmosphere chosen: %s", atm_type)

        # Try fix instrument naming
        instr = INSTRUMENTS.get(platform_name, sensor)
        if instr != sensor:
            sensor = instr
            LOG.warning("Inconsistent sensor/satellite input - " +
                        "sensor set to %s", sensor)

        self.sensor = sensor.replace('/', '')

        ext = atm_type.replace(' ', '_')
        lutname = "rayleigh_lut_{0}.h5".format(ext)
        self.reflectance_lut_filename = os.path.join(rayleigh_dir, lutname)
        if not os.path.exists(self.reflectance_lut_filename):
            LOG.warning(
                "No lut file %s on disk", self.reflectance_lut_filename)
            LOG.info("Will download from internet...")
            download_luts(aerosol_type=aerosol_type)

        if (not os.path.exists(self.reflectance_lut_filename) or
                not os.path.isfile(self.reflectance_lut_filename)):
            raise IOError('pyspectral file for Rayleigh scattering correction ' +
                          'does not exist! Filename = ' +
                          str(self.reflectance_lut_filename))

        LOG.debug('LUT filename: %s', str(self.reflectance_lut_filename))
        self._rayl = None
        self._wvl_coord = None
        self._azid_coord = None
        self._satz_sec_coord = None
        self._sunz_sec_coord = None

    def get_effective_wavelength(self, bandname):
        """Get the effective wavelength with Rayleigh scattering in mind"""
        try:
            rsr = RelativeSpectralResponse(self.platform_name, self.sensor)
        except IOError:
            LOG.exception(
                "No spectral responses for this platform and sensor: %s %s", self.platform_name, self.sensor)
            if isinstance(bandname, (float, integer_types)):
                LOG.warning(
                    "Effective wavelength is set to the requested band wavelength = %f", bandname)
                return bandname
            raise

        if isinstance(bandname, str):
            bandname = BANDNAMES.get(self.sensor, BANDNAMES['generic']).get(bandname, bandname)
        elif isinstance(bandname, (float, integer_types)):
            if not(0.4 < bandname < 0.8):
                raise BandFrequencyOutOfRange(
                    'Requested band frequency should be between 0.4 and 0.8 microns!')
            bandname = get_bandname_from_wavelength(self.sensor, bandname, rsr.rsr)

        wvl, resp = rsr.rsr[bandname][
            'det-1']['wavelength'], rsr.rsr[bandname]['det-1']['response']

        cwvl = get_central_wave(wvl, resp, weight=1. / wvl**4)
        LOG.debug("Band name: %s  Effective wavelength: %f", bandname, cwvl)

        return cwvl

    def get_reflectance_lut(self):
        """Read the LUT with reflectances as a function of wavelength, satellite zenith
        secant, azimuth difference angle, and sun zenith secant

        """
        if self._rayl is None:
            lut_vars = get_reflectance_lut(self.reflectance_lut_filename)
            self._rayl = lut_vars[0]
            self._wvl_coord = lut_vars[1]
            self._azid_coord = lut_vars[2]
            self._satz_sec_coord = lut_vars[3]
            self._sunz_sec_coord = lut_vars[4]
        return self._rayl, self._wvl_coord, self._azid_coord,\
            self._satz_sec_coord, self._sunz_sec_coord

    def get_reflectance(self, sun_zenith, sat_zenith, azidiff, bandname,
                        redband=None):
        """Get the reflectance from the three sun-sat angles."""
        # Get wavelength in nm for band:
        wvl = self.get_effective_wavelength(bandname) * 1000.0
        rayl, wvl_coord, azid_coord, satz_sec_coord, sunz_sec_coord = \
            self.get_reflectance_lut()

        clip_angle = np.rad2deg(np.arccos(1. / sunz_sec_coord.max()))
        sun_zenith = clip(sun_zenith, 0, clip_angle)
        sunzsec = 1. / np.cos(np.deg2rad(sun_zenith))
        clip_angle = np.rad2deg(np.arccos(1. / satz_sec_coord.max()))
        sat_zenith = clip(np.asarray(sat_zenith), 0, clip_angle)
        satzsec = 1. / np.cos(np.deg2rad(sat_zenith))

        shape = sun_zenith.shape

        if not(wvl_coord.min() < wvl < wvl_coord.max()):
            LOG.warning(
                "Effective wavelength for band %s outside 400-800 nm range!",
                str(bandname))
            LOG.info(
                "Set the rayleigh/aerosol reflectance contribution to zero!")
            chunks = sun_zenith.chunks if redband is None else redband.chunks
            return zeros(shape, chunks=chunks)

        idx = np.searchsorted(wvl_coord, wvl)
        wvl1 = wvl_coord[idx - 1]
        wvl2 = wvl_coord[idx]

        fac = (wvl2 - wvl) / (wvl2 - wvl1)
        raylwvl = fac * rayl[idx - 1, :, :, :] + (1 - fac) * rayl[idx, :, :, :]
        tic = time.time()

        smin = [sunz_sec_coord[0], azid_coord[0], satz_sec_coord[0]]
        smax = [sunz_sec_coord[-1], azid_coord[-1], satz_sec_coord[-1]]
        orders = [
            len(sunz_sec_coord), len(azid_coord), len(satz_sec_coord)]
        minterp = MultilinearInterpolator(smin, smax, orders)

        f_3d_grid = raylwvl
        minterp.set_values(np.atleast_2d(f_3d_grid.ravel()))

        def _do_interp(minterp, sunzsec, azidiff, satzsec):
            interp_points2 = np.vstack((sunzsec.ravel(),
                                        180 - azidiff.ravel(),
                                        satzsec.ravel()))
            res = minterp(interp_points2)
            return res.reshape(sunzsec.shape)

        if HAVE_DASK:
            ipn = map_blocks(_do_interp, minterp, sunzsec, azidiff,
                             satzsec, dtype=raylwvl.dtype,
                             chunks=azidiff.chunks)
        else:
            ipn = _do_interp(minterp, sunzsec, azidiff, satzsec)

        LOG.debug("Time - Interpolation: {0:f}".format(time.time() - tic))

        ipn *= 100
        res = ipn
        if redband is not None:
            res = where(redband < 20., res,
                        (1 - (redband - 20) / 80) * res)

        return clip(res, 0, 100)


def get_reflectance_lut(filename):
    """Read the LUT with reflectances as a function of wavelength, satellite
    zenith secant, azimuth difference angle, and sun zenith secant

    """

    h5f = h5py.File(filename, 'r')

    tab = h5f['reflectance']
    wvl = h5f['wavelengths']
    azidiff = h5f['azimuth_difference']
    satellite_zenith_secant = h5f['satellite_zenith_secant']
    sun_zenith_secant = h5f['sun_zenith_secant']

    if HAVE_DASK:
        tab = from_array(tab, chunks=(10, 10, 10, 10))
        # wvl_coord is used in a lot of non-dask functions, keep in memory
        wvl = from_array(wvl, chunks=(100,)).persist()
        azidiff = from_array(azidiff, chunks=(1000,))
        satellite_zenith_secant = from_array(satellite_zenith_secant,
                                             chunks=(1000,))
        sun_zenith_secant = from_array(sun_zenith_secant,
                                       chunks=(1000,))
    else:
        # load all of the data we are going to use in to memory
        tab = tab[:]
        wvl = wvl[:]
        azidiff = azidiff[:]
        satellite_zenith_secant = satellite_zenith_secant[:]
        sun_zenith_secant = sun_zenith_secant[:]

    return tab, wvl, azidiff, satellite_zenith_secant, sun_zenith_secant


# if __name__ == "__main__":
#     SHAPE = (1000, 3000)
#     NDIM = SHAPE[0] * SHAPE[1]
#     SUNZ = np.ma.arange(
#         NDIM / 2, NDIM + NDIM / 2).reshape(SHAPE) * 60. / float(NDIM)
#     SATZ = np.ma.arange(NDIM).reshape(SHAPE) * 60. / float(NDIM)
#     AZIDIFF = np.ma.arange(NDIM).reshape(SHAPE) * 179.9 / float(NDIM)
#     rfl = this.get_reflectance(SUNZ, SATZ, AZIDIFF, 'M4')


1			#!/usr/bin/env python
2			# -- coding: utf-8 --
3
4			# Copyright (c) 2016-2018 Pytroll
5
6			# Author(s):
7
8			# Adam.Dybbroe <[email protected]>
9			# Martin Raspaud <[email protected]>
10
11			# This program is free software: you can redistribute it and/or modify
12			# it under the terms of the GNU General Public License as published by
13			# the Free Software Foundation, either version 3 of the License, or
14			# (at your option) any later version.
15
16			# This program is distributed in the hope that it will be useful,
17			# but WITHOUT ANY WARRANTY; without even the implied warranty of
18			# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19			# GNU General Public License for more details.
20
21			# You should have received a copy of the GNU General Public License
22			# along with this program. If not, see <http://www.gnu.org/licenses/>.
23
24			"""Atmospheric correction of shortwave imager bands in the wavelength range 400
25			to 800 nm
26
27			"""
28
29			import os
30			import time
31			import logging
32			from six import integer_types
33
34			import h5py
35			import numpy as np
36
37			try:
38			from dask.array import where, zeros, map_blocks, from_array, clip
39			HAVE_DASK = True
40			except ImportError:
41			from numpy import where, zeros, clip
42			map_blocks = None
43			from_array = None
44			HAVE_DASK = False
45
46			from geotiepoints.multilinear import MultilinearInterpolator
47			from pyspectral.rsr_reader import RelativeSpectralResponse
48			from pyspectral.utils import (INSTRUMENTS, RAYLEIGH_LUT_DIRS,
49			AEROSOL_TYPES, ATMOSPHERES,
50			BANDNAMES,
51			download_luts, get_central_wave,
52			get_bandname_from_wavelength)
53
54			LOG = logging.getLogger(__name__)
55
56
57			class BandFrequencyOutOfRange(ValueError):
58
59			"""Exception when the band frequency is out of the visible range"""
60
61			pass
62
63
64			class Rayleigh(object):
65
66			"""Container for the atmospheric correction of satellite imager bands.
67
68			This class removes background contributions of Rayleigh scattering of
69			molecules and Mie scattering and absorption by aerosols.
70
71			"""
72
73			def __init__(self, platform_name, sensor, **kwargs):
74			"""Initialize class and determine LUT to use."""
75			self.platform_name = platform_name
76			self.sensor = sensor
77			self.coeff_filename = None
78
79			atm_type = kwargs.get('atmosphere', 'us-standard')
80			if atm_type not in ATMOSPHERES:
81			raise AttributeError('Atmosphere type not supported! ' +
82			'Need to be one of {}'.format(str(ATMOSPHERES)))
83
84			aerosol_type = kwargs.get('aerosol_type', 'marine_clean_aerosol')
85
86			if aerosol_type not in AEROSOL_TYPES:
87			raise AttributeError('Aerosol type not supported! ' +
88			'Need to be one of {0}'.format(str(AEROSOL_TYPES)))
89
90			rayleigh_dir = RAYLEIGH_LUT_DIRS[aerosol_type]
91
92			if atm_type not in ATMOSPHERES.keys():
93			LOG.error("Atmosphere type %s not supported", atm_type)
94
95			LOG.info("Atmosphere chosen: %s", atm_type)
96
97			# Try fix instrument naming
98			instr = INSTRUMENTS.get(platform_name, sensor)
99			if instr != sensor:
100			sensor = instr
101			LOG.warning("Inconsistent sensor/satellite input - " +
102			"sensor set to %s", sensor)
103
104			self.sensor = sensor.replace('/', '')
105
106			ext = atm_type.replace(' ', '_')
107			lutname = "rayleigh_lut_{0}.h5".format(ext)
108			self.reflectance_lut_filename = os.path.join(rayleigh_dir, lutname)
109			if not os.path.exists(self.reflectance_lut_filename):
110			LOG.warning(
111			"No lut file %s on disk", self.reflectance_lut_filename)
112			LOG.info("Will download from internet...")
113			download_luts(aerosol_type=aerosol_type)
114
115			if (not os.path.exists(self.reflectance_lut_filename) or
116			not os.path.isfile(self.reflectance_lut_filename)):
117			raise IOError('pyspectral file for Rayleigh scattering correction ' +
118			'does not exist! Filename = ' +
119			str(self.reflectance_lut_filename))
120
121			LOG.debug('LUT filename: %s', str(self.reflectance_lut_filename))
122			self._rayl = None
123			self._wvl_coord = None
124			self._azid_coord = None
125			self._satz_sec_coord = None
126			self._sunz_sec_coord = None
127
128			def get_effective_wavelength(self, bandname):
129			"""Get the effective wavelength with Rayleigh scattering in mind"""
130			try:
131			rsr = RelativeSpectralResponse(self.platform_name, self.sensor)
132			except IOError:
133			LOG.exception(
134			"No spectral responses for this platform and sensor: %s %s", self.platform_name, self.sensor)
135			if isinstance(bandname, (float, integer_types)):
136			LOG.warning(
137			"Effective wavelength is set to the requested band wavelength = %f", bandname)
138			return bandname
139			raise
140
141			if isinstance(bandname, str):
142			bandname = BANDNAMES.get(self.sensor, BANDNAMES['generic']).get(bandname, bandname)
143			elif isinstance(bandname, (float, integer_types)):
144			if not(0.4 < bandname < 0.8):
145			raise BandFrequencyOutOfRange(
146			'Requested band frequency should be between 0.4 and 0.8 microns!')
147			bandname = get_bandname_from_wavelength(self.sensor, bandname, rsr.rsr)
148
149			wvl, resp = rsr.rsr[bandname][
150			'det-1']['wavelength'], rsr.rsr[bandname]['det-1']['response']
151
152			cwvl = get_central_wave(wvl, resp, weight=1. / wvl**4)
153			LOG.debug("Band name: %s Effective wavelength: %f", bandname, cwvl)
154
155			return cwvl
156
157			def get_reflectance_lut(self):
158			"""Read the LUT with reflectances as a function of wavelength, satellite zenith
159			secant, azimuth difference angle, and sun zenith secant
160
161			"""
162			if self._rayl is None:
163			lut_vars = get_reflectance_lut(self.reflectance_lut_filename)
164			self._rayl = lut_vars[0]
165			self._wvl_coord = lut_vars[1]
166			self._azid_coord = lut_vars[2]
167			self._satz_sec_coord = lut_vars[3]
168			self._sunz_sec_coord = lut_vars[4]
169			return self._rayl, self._wvl_coord, self._azid_coord,\
170			self._satz_sec_coord, self._sunz_sec_coord
171
172			def get_reflectance(self, sun_zenith, sat_zenith, azidiff, bandname,
173			redband=None):
174			"""Get the reflectance from the three sun-sat angles."""
175			# Get wavelength in nm for band:
176			wvl = self.get_effective_wavelength(bandname) * 1000.0
177			rayl, wvl_coord, azid_coord, satz_sec_coord, sunz_sec_coord = \
178			self.get_reflectance_lut()
179
180			clip_angle = np.rad2deg(np.arccos(1. / sunz_sec_coord.max()))
181			sun_zenith = clip(sun_zenith, 0, clip_angle)
182			sunzsec = 1. / np.cos(np.deg2rad(sun_zenith))
183			clip_angle = np.rad2deg(np.arccos(1. / satz_sec_coord.max()))
184			sat_zenith = clip(np.asarray(sat_zenith), 0, clip_angle)
185			satzsec = 1. / np.cos(np.deg2rad(sat_zenith))
186
187			shape = sun_zenith.shape
188
189			if not(wvl_coord.min() < wvl < wvl_coord.max()):
190			LOG.warning(
191			"Effective wavelength for band %s outside 400-800 nm range!",
192			str(bandname))
193			LOG.info(
194			"Set the rayleigh/aerosol reflectance contribution to zero!")
195			chunks = sun_zenith.chunks if redband is None else redband.chunks
196			return zeros(shape, chunks=chunks)
197
198			idx = np.searchsorted(wvl_coord, wvl)
199			wvl1 = wvl_coord[idx - 1]
200			wvl2 = wvl_coord[idx]
201
202			fac = (wvl2 - wvl) / (wvl2 - wvl1)
203			raylwvl = fac * rayl[idx - 1, :, :, :] + (1 - fac) * rayl[idx, :, :, :]
204			tic = time.time()
205
206			smin = [sunz_sec_coord[0], azid_coord[0], satz_sec_coord[0]]
207			smax = [sunz_sec_coord[-1], azid_coord[-1], satz_sec_coord[-1]]
208			orders = [
209			len(sunz_sec_coord), len(azid_coord), len(satz_sec_coord)]
210			minterp = MultilinearInterpolator(smin, smax, orders)
211
212			f_3d_grid = raylwvl
213			minterp.set_values(np.atleast_2d(f_3d_grid.ravel()))
214
215			def _do_interp(minterp, sunzsec, azidiff, satzsec):
216			interp_points2 = np.vstack((sunzsec.ravel(),
217			180 - azidiff.ravel(),
218			satzsec.ravel()))
219			res = minterp(interp_points2)
220			return res.reshape(sunzsec.shape)
221
222			if HAVE_DASK:
223			ipn = map_blocks(_do_interp, minterp, sunzsec, azidiff,
224			satzsec, dtype=raylwvl.dtype,
225			chunks=azidiff.chunks)
226			else:
227			ipn = _do_interp(minterp, sunzsec, azidiff, satzsec)
228
229			LOG.debug("Time - Interpolation: {0:f}".format(time.time() - tic))
230
231			ipn *= 100
232			res = ipn
233			if redband is not None:
234			res = where(redband < 20., res,
235			(1 - (redband - 20) / 80) * res)
236
237			return clip(res, 0, 100)
238
239
240			def get_reflectance_lut(filename):
241			"""Read the LUT with reflectances as a function of wavelength, satellite
242			zenith secant, azimuth difference angle, and sun zenith secant
243
244			"""
245
246			h5f = h5py.File(filename, 'r')
247
248			tab = h5f['reflectance']
249			wvl = h5f['wavelengths']
250			azidiff = h5f['azimuth_difference']
251			satellite_zenith_secant = h5f['satellite_zenith_secant']
252			sun_zenith_secant = h5f['sun_zenith_secant']
253
254			if HAVE_DASK:
255			tab = from_array(tab, chunks=(10, 10, 10, 10))
256			# wvl_coord is used in a lot of non-dask functions, keep in memory
257			wvl = from_array(wvl, chunks=(100,)).persist()
258			azidiff = from_array(azidiff, chunks=(1000,))
259			satellite_zenith_secant = from_array(satellite_zenith_secant,
260			chunks=(1000,))
261			sun_zenith_secant = from_array(sun_zenith_secant,
262			chunks=(1000,))
263			else:
264			# load all of the data we are going to use in to memory
265			tab = tab[:]
266			wvl = wvl[:]
267			azidiff = azidiff[:]
268			satellite_zenith_secant = satellite_zenith_secant[:]
269			sun_zenith_secant = sun_zenith_secant[:]
270
271			return tab, wvl, azidiff, satellite_zenith_secant, sun_zenith_secant
272
273
274			# if __name__ == "__main__":
275			# SHAPE = (1000, 3000)
276			# NDIM = SHAPE[0] * SHAPE[1]
277			# SUNZ = np.ma.arange(
278			# NDIM / 2, NDIM + NDIM / 2).reshape(SHAPE) * 60. / float(NDIM)
279			# SATZ = np.ma.arange(NDIM).reshape(SHAPE) * 60. / float(NDIM)
280			# AZIDIFF = np.ma.arange(NDIM).reshape(SHAPE) * 179.9 / float(NDIM)
281			# rfl = this.get_reflectance(SUNZ, SATZ, AZIDIFF, 'M4')
282

pytroll / pyspectral

Pull Request — develop (#33)

get_reflectance_lut() B

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like