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:07 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	53
rs	5.8625

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 short
    wave bands. Removing background contributions of Rayleigh scattering of
    molecules and Mie scattering and absorption by aerosols.

    """

    def __init__(self, platform_name, sensor, **kwargs):
        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 short
67			wave bands. Removing background contributions of Rayleigh scattering of
68			molecules and Mie scattering and absorption by aerosols.
69
70			"""
71
72			def __init__(self, platform_name, sensor, **kwargs):
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
129			try:
130			rsr = RelativeSpectralResponse(self.platform_name, self.sensor)
131			except IOError:
132			LOG.exception(
133			"No spectral responses for this platform and sensor: %s %s", self.platform_name, self.sensor)
134			if isinstance(bandname, (float, integer_types)):
135			LOG.warning(
136			"Effective wavelength is set to the requested band wavelength = %f", bandname)
137			return bandname
138			raise
139
140			if isinstance(bandname, str):
141			bandname = BANDNAMES.get(self.sensor, BANDNAMES['generic']).get(bandname, bandname)
142			elif isinstance(bandname, (float, integer_types)):
143			if not(0.4 < bandname < 0.8):
144			raise BandFrequencyOutOfRange(
145			'Requested band frequency should be between 0.4 and 0.8 microns!')
146			bandname = get_bandname_from_wavelength(self.sensor, bandname, rsr.rsr)
147
148			wvl, resp = rsr.rsr[bandname][
149			'det-1']['wavelength'], rsr.rsr[bandname]['det-1']['response']
150
151			cwvl = get_central_wave(wvl, resp, weight=1. / wvl**4)
152			LOG.debug("Band name: %s Effective wavelength: %f", bandname, cwvl)
153
154			return cwvl
155
156			def get_reflectance_lut(self):
157			"""Read the LUT with reflectances as a function of wavelength, satellite zenith
158			secant, azimuth difference angle, and sun zenith secant
159
160			"""
161			if self._rayl is None:
162			lut_vars = get_reflectance_lut(self.reflectance_lut_filename)
163			self._rayl = lut_vars[0]
164			self._wvl_coord = lut_vars[1]
165			self._azid_coord = lut_vars[2]
166			self._satz_sec_coord = lut_vars[3]
167			self._sunz_sec_coord = lut_vars[4]
168			return self._rayl, self._wvl_coord, self._azid_coord,\
169			self._satz_sec_coord, self._sunz_sec_coord
170
171			def get_reflectance(self, sun_zenith, sat_zenith, azidiff, bandname,
172			redband=None):
173			"""Get the reflectance from the three sun-sat angles."""
174			# Get wavelength in nm for band:
175			wvl = self.get_effective_wavelength(bandname) * 1000.0
176			rayl, wvl_coord, azid_coord, satz_sec_coord, sunz_sec_coord = \
177			self.get_reflectance_lut()
178
179			clip_angle = np.rad2deg(np.arccos(1. / sunz_sec_coord.max()))
180			sun_zenith = clip(sun_zenith, 0, clip_angle)
181			sunzsec = 1. / np.cos(np.deg2rad(sun_zenith))
182			clip_angle = np.rad2deg(np.arccos(1. / satz_sec_coord.max()))
183			sat_zenith = clip(np.asarray(sat_zenith), 0, clip_angle)
184			satzsec = 1. / np.cos(np.deg2rad(sat_zenith))
185
186			shape = sun_zenith.shape
187
188			if not(wvl_coord.min() < wvl < wvl_coord.max()):
189			LOG.warning(
190			"Effective wavelength for band %s outside 400-800 nm range!",
191			str(bandname))
192			LOG.info(
193			"Set the rayleigh/aerosol reflectance contribution to zero!")
194			chunks = sun_zenith.chunks if redband is None else redband.chunks
195			return zeros(shape, chunks=chunks)
196
197			idx = np.searchsorted(wvl_coord, wvl)
198			wvl1 = wvl_coord[idx - 1]
199			wvl2 = wvl_coord[idx]
200
201			fac = (wvl2 - wvl) / (wvl2 - wvl1)
202			raylwvl = fac * rayl[idx - 1, :, :, :] + (1 - fac) * rayl[idx, :, :, :]
203			tic = time.time()
204
205			smin = [sunz_sec_coord[0], azid_coord[0], satz_sec_coord[0]]
206			smax = [sunz_sec_coord[-1], azid_coord[-1], satz_sec_coord[-1]]
207			orders = [
208			len(sunz_sec_coord), len(azid_coord), len(satz_sec_coord)]
209			minterp = MultilinearInterpolator(smin, smax, orders)
210
211			f_3d_grid = raylwvl
212			minterp.set_values(np.atleast_2d(f_3d_grid.ravel()))
213
214			def _do_interp(minterp, sunzsec, azidiff, satzsec):
215			interp_points2 = np.vstack((sunzsec.ravel(),
216			180 - azidiff.ravel(),
217			satzsec.ravel()))
218			res = minterp(interp_points2)
219			return res.reshape(sunzsec.shape)
220
221			if HAVE_DASK:
222			ipn = map_blocks(_do_interp, minterp, sunzsec, azidiff,
223			satzsec, dtype=raylwvl.dtype,
224			chunks=azidiff.chunks)
225			else:
226			ipn = _do_interp(minterp, sunzsec, azidiff, satzsec)
227
228			LOG.debug("Time - Interpolation: {0:f}".format(time.time() - tic))
229
230			ipn *= 100
231			res = ipn
232			if redband is not None:
233			res = where(redband < 20., res,
234			(1 - (redband - 20) / 80) * res)
235
236			return clip(res, 0, 100)
237
238
239			def get_reflectance_lut(filename):
240			"""Read the LUT with reflectances as a function of wavelength, satellite
241			zenith secant, azimuth difference angle, and sun zenith secant
242
243			"""
244
245			h5f = h5py.File(filename, 'r')
246
247			tab = h5f['reflectance']
248			wvl = h5f['wavelengths']
249			azidiff = h5f['azimuth_difference']
250			satellite_zenith_secant = h5f['satellite_zenith_secant']
251			sun_zenith_secant = h5f['sun_zenith_secant']
252
253			if HAVE_DASK:
254			tab = from_array(tab, chunks=(10, 10, 10, 10))
255			# wvl_coord is used in a lot of non-dask functions, keep in memory
256			wvl = from_array(wvl, chunks=(100,)).persist()
257			azidiff = from_array(azidiff, chunks=(1000,))
258			satellite_zenith_secant = from_array(satellite_zenith_secant,
259			chunks=(1000,))
260			sun_zenith_secant = from_array(sun_zenith_secant,
261			chunks=(1000,))
262			else:
263			# load all of the data we are going to use in to memory
264			tab = tab[:]
265			wvl = wvl[:]
266			azidiff = azidiff[:]
267			satellite_zenith_secant = satellite_zenith_secant[:]
268			sun_zenith_secant = sun_zenith_secant[:]
269
270			return tab, wvl, azidiff, satellite_zenith_secant, sun_zenith_secant
271
272
273			# if __name__ == "__main__":
274			# SHAPE = (1000, 3000)
275			# NDIM = SHAPE[0] * SHAPE[1]
276			# SUNZ = np.ma.arange(
277			# NDIM / 2, NDIM + NDIM / 2).reshape(SHAPE) * 60. / float(NDIM)
278			# SATZ = np.ma.arange(NDIM).reshape(SHAPE) * 60. / float(NDIM)
279			# AZIDIFF = np.ma.arange(NDIM).reshape(SHAPE) * 179.9 / float(NDIM)
280			# rfl = this.get_reflectance(SUNZ, SATZ, AZIDIFF, 'M4')
281

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