Rayleigh.get_effective_wavelength() - 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-09 15:50 UTC

Rayleigh.get_effective_wavelength() B

↳ Parent: Rayleigh

Complexity

Conditions

Size

Total Lines

Duplication

Lines	0
Ratio	0 %

Importance

Changes	2
Bugs	1	Features	0

Metric	Value
cc	6
c	2
b	1
f	0
dl	0
loc	28
rs	7.5384

#!/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, Array
    HAVE_DASK = True
except ImportError:
    from numpy import where, zeros, clip
    map_blocks = None
    from_array = None
    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()

        # force dask arrays
        compute = False
        if HAVE_DASK and not isinstance(sun_zenith, Array):
            compute = True
            sun_zenith = from_array(sun_zenith, chunks=sun_zenith.shape)
            sat_zenith = from_array(sat_zenith, chunks=sat_zenith.shape)
            azidiff = from_array(azidiff, chunks=azidiff.shape)
            if redband is not None:
                redband = from_array(redband, chunks=redband.shape)

        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(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!")
            if HAVE_DASK:
                chunks = sun_zenith.chunks if redband is None else redband.chunks
                res = zeros(shape, chunks=chunks)
                return res.compute() if compute else res
            else:
                return zeros(shape)

        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)

        res = clip(res, 0, 100)
        if compute:
            res = res.compute()
        return res


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

pytroll / pyspectral

Pull Request — develop (#33)

Rayleigh.get_effective_wavelength() B

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like