Rayleigh.__init__() - 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:23 UTC

Rayleigh.init() C

↳ Parent: Rayleigh

Complexity

Conditions

Size

Total Lines

Duplication

Lines	0
Ratio	0 %

Importance

Changes	2
Bugs	0	Features	0

Metric	Value
cc	8
c	2
b	0
f	0
dl	0
loc	54
rs	6.2167

How to fix Long Method

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

pytroll / pyspectral

Pull Request — develop (#33)

Rayleigh.__init__() C

Complexity

Size

Duplication

Importance

How to fix Long Method

Long Method

Duplication Side-by-Side

Filter issues like

Rayleigh.init() C