get_reflectance_lut() - Code Metrics - Inspection of "Improve unit test coverage" - pytroll/pyspectral - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Push — develop ( f92b69...710c13 )

by Adam

created 2018-01-08 10:48 UTC

get_reflectance_lut() A

↳ Parent: Project

Complexity

Conditions

Size

Total Lines

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
cc	2
dl	0
loc	14
rs	9.4285
c	0
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/>.

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

"""

import logging
import os
from six import integer_types

import h5py
import numpy as np
from scipy.interpolate import interpn

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__)

WITH_CYTHON = True
try:
    from geotiepoints.multilinear import MultilinearInterpolator
except ImportError:
    LOG.warning(
        "Couldn't import fast multilinear interpolation with Cython.")
    LOG.warning("Check your geotiepoints installation!")
    WITH_CYTHON = False


class BandFrequencyOutOfRange(Exception):

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


class Rayleigh(object):

    """Container for the Rayleigh scattering correction of satellite imager short
    wave bands

    """

    def __init__(self, platform_name, sensor, **kwargs):
        self.platform_name = platform_name
        self.sensor = sensor
        self.coeff_filename = None

        atm_type = kwargs.get('atmosphere', 'subarctic summer')
        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', 'rayleigh_only')

        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))

    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(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(bandname, rsr.rsr)

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

        return get_central_wave(wvl, resp, weight=1. / wvl**4)

    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

        """

        return get_reflectance_lut(self.reflectance_lut_filename)

    def get_reflectance(self, sun_zenith, sat_zenith, azidiff, bandname,
                        blueband=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 = np.clip(np.asarray(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 = np.clip(np.asarray(sat_zenith), 0, clip_angle)
        satzsec = 1. / np.cos(np.deg2rad(np.asarray(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!")
            return np.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, :, :, :]

        import time
        tic = time.time()

        if WITH_CYTHON:
            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)]
            f_3d_grid = raylwvl

            minterp = MultilinearInterpolator(smin, smax, orders)
            minterp.set_values(np.atleast_2d(f_3d_grid.ravel()))

            interp_points2 = np.vstack((np.asarray(sunzsec).ravel(),
                                        np.asarray(180 - azidiff).ravel(),
                                        np.asarray(satzsec).ravel()))

            ipn = minterp(interp_points2).reshape(shape)
        else:
            interp_points = np.dstack((np.asarray(sunzsec),
                                       np.asarray(180. - azidiff),
                                       np.asarray(satzsec)))

            ipn = interpn((sunz_sec_coord, azid_coord, satz_sec_coord),
                          raylwvl[:, :, :], interp_points)

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

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

        return np.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

    """

    with h5py.File(filename, 'r') as h5f:
        tab = h5f['reflectance'][:]
        wvl = h5f['wavelengths'][:]
        azidiff = h5f['azimuth_difference'][:]
        satellite_zenith_secant = h5f['satellite_zenith_secant'][:]
        sun_zenith_secant = h5f['sun_zenith_secant'][:]

    return tab, wvl, azidiff, satellite_zenith_secant, sun_zenith_secant

if __name__ == "__main__":

    this = Rayleigh('Suomi-NPP', 'viirs')
    # SUNZ = np.arange(200000).reshape(400, 500) * 0.0004
    # SATZ = np.arange(200000).reshape(400, 500) * 0.00025
    # AZIDIFF = np.arange(200000).reshape(400, 500) * 0.0009
    # rfl = this.get_reflectance(SUNZ, SATZ, AZIDIFF, 'M4')

    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			"""Rayleigh correction of shortwave imager bands in the wavelength range 400 to
25			800 nm
26
27			"""
28
29			import logging
30			import os
31			from six import integer_types
32
33			import h5py
34			import numpy as np
35			from scipy.interpolate import interpn
36
37			from pyspectral.rsr_reader import RelativeSpectralResponse
38			from pyspectral.utils import (INSTRUMENTS, RAYLEIGH_LUT_DIRS,
39			AEROSOL_TYPES, ATMOSPHERES,
40			BANDNAMES,
41			download_luts, get_central_wave,
42			get_bandname_from_wavelength)
43
44			LOG = logging.getLogger(__name__)
45
46			WITH_CYTHON = True
47			try:
48			from geotiepoints.multilinear import MultilinearInterpolator
49			except ImportError:
50			LOG.warning(
51			"Couldn't import fast multilinear interpolation with Cython.")
52			LOG.warning("Check your geotiepoints installation!")
53			WITH_CYTHON = False
54
55
56			class BandFrequencyOutOfRange(Exception):
57
58			"""Exception when the band frequency is out of the visible range"""
59			pass
60
61
62			class Rayleigh(object):
63
64			"""Container for the Rayleigh scattering correction of satellite imager short
65			wave bands
66
67			"""
68
69			def __init__(self, platform_name, sensor, **kwargs):
70			self.platform_name = platform_name
71			self.sensor = sensor
72			self.coeff_filename = None
73
74			atm_type = kwargs.get('atmosphere', 'subarctic summer')
75			if atm_type not in ATMOSPHERES:
76			raise AttributeError('Atmosphere type not supported! ' +
77			'Need to be one of {}'.format(str(ATMOSPHERES)))
78
79			aerosol_type = kwargs.get('aerosol_type', 'rayleigh_only')
80
81			if aerosol_type not in AEROSOL_TYPES:
82			raise AttributeError('Aerosol type not supported! ' +
83			'Need to be one of {0}'.format(str(AEROSOL_TYPES)))
84
85			rayleigh_dir = RAYLEIGH_LUT_DIRS[aerosol_type]
86
87			if atm_type not in ATMOSPHERES.keys():
88			LOG.error("Atmosphere type %s not supported", atm_type)
89
90			LOG.info("Atmosphere chosen: %s", atm_type)
91
92			# Try fix instrument naming
93			instr = INSTRUMENTS.get(platform_name, sensor)
94			if instr != sensor:
95			sensor = instr
96			LOG.warning("Inconsistent sensor/satellite input - " +
97			"sensor set to %s", sensor)
98
99			self.sensor = sensor.replace('/', '')
100
101			ext = atm_type.replace(' ', '_')
102			lutname = "rayleigh_lut_{0}.h5".format(ext)
103			self.reflectance_lut_filename = os.path.join(rayleigh_dir, lutname)
104			if not os.path.exists(self.reflectance_lut_filename):
105			LOG.warning(
106			"No lut file %s on disk", self.reflectance_lut_filename)
107			LOG.info("Will download from internet...")
108			download_luts(aerosol_type=aerosol_type)
109
110			if (not os.path.exists(self.reflectance_lut_filename) or
111			not os.path.isfile(self.reflectance_lut_filename)):
112			raise IOError('pyspectral file for Rayleigh scattering correction ' +
113			'does not exist! Filename = ' +
114			str(self.reflectance_lut_filename))
115
116			LOG.debug('LUT filename: %s', str(self.reflectance_lut_filename))
117
118			def get_effective_wavelength(self, bandname):
119			"""Get the effective wavelength with Rayleigh scattering in mind"""
120
121			try:
122			rsr = RelativeSpectralResponse(self.platform_name, self.sensor)
123			except IOError:
124			LOG.exception(
125			"No spectral responses for this platform and sensor: %s %s", self.platform_name, self.sensor)
126			if isinstance(bandname, (float, integer_types)):
127			LOG.warning(
128			"Effective wavelength is set to the requested band wavelength = %f", bandname)
129			return bandname
130			raise
131
132			if isinstance(bandname, str):
133			bandname = BANDNAMES.get(bandname, bandname)
134			elif isinstance(bandname, (float, integer_types)):
135			if not(0.4 < bandname < 0.8):
136			raise BandFrequencyOutOfRange(
137			'Requested band frequency should be between 0.4 and 0.8 microns!')
138			bandname = get_bandname_from_wavelength(bandname, rsr.rsr)
139
140			wvl, resp = rsr.rsr[bandname][
141			'det-1']['wavelength'], rsr.rsr[bandname]['det-1']['response']
142
143			return get_central_wave(wvl, resp, weight=1. / wvl**4)
144
145			def get_reflectance_lut(self):
146			"""Read the LUT with reflectances as a function of wavelength, satellite zenith
147			secant, azimuth difference angle, and sun zenith secant
148
149			"""
150
151			return get_reflectance_lut(self.reflectance_lut_filename)
152
153			def get_reflectance(self, sun_zenith, sat_zenith, azidiff, bandname,
154			blueband=None):
155			"""Get the reflectance from the three sun-sat angles."""
156			# Get wavelength in nm for band:
157			wvl = self.get_effective_wavelength(bandname) * 1000.0
158			rayl, wvl_coord, azid_coord, satz_sec_coord, sunz_sec_coord = self.get_reflectance_lut()
159
160			clip_angle = np.rad2deg(np.arccos(1. / sunz_sec_coord.max()))
161			sun_zenith = np.clip(np.asarray(sun_zenith), 0, clip_angle)
162			sunzsec = 1. / np.cos(np.deg2rad(sun_zenith))
163			clip_angle = np.rad2deg(np.arccos(1. / satz_sec_coord.max()))
164			sat_zenith = np.clip(np.asarray(sat_zenith), 0, clip_angle)
165			satzsec = 1. / np.cos(np.deg2rad(np.asarray(sat_zenith)))
166
167			shape = sun_zenith.shape
168
169			if not(wvl_coord.min() < wvl < wvl_coord.max()):
170			LOG.warning(
171			"Effective wavelength for band %s outside 400-800 nm range!", str(bandname))
172			LOG.info(
173			"Set the rayleigh/aerosol reflectance contribution to zero!")
174			return np.zeros(shape)
175
176			idx = np.searchsorted(wvl_coord, wvl)
177			wvl1 = wvl_coord[idx - 1]
178			wvl2 = wvl_coord[idx]
179
180			fac = (wvl2 - wvl) / (wvl2 - wvl1)
181
182			raylwvl = fac * rayl[idx - 1, :, :, :] + (1 - fac) * rayl[idx, :, :, :]
183
184			import time
185			tic = time.time()
186
187			if WITH_CYTHON:
188			smin = [sunz_sec_coord[0], azid_coord[0], satz_sec_coord[0]]
189			smax = [sunz_sec_coord[-1], azid_coord[-1], satz_sec_coord[-1]]
190			orders = [
191			len(sunz_sec_coord), len(azid_coord), len(satz_sec_coord)]
192			f_3d_grid = raylwvl
193
194			minterp = MultilinearInterpolator(smin, smax, orders)
195			minterp.set_values(np.atleast_2d(f_3d_grid.ravel()))
196
197			interp_points2 = np.vstack((np.asarray(sunzsec).ravel(),
198			np.asarray(180 - azidiff).ravel(),
199			np.asarray(satzsec).ravel()))
200
201			ipn = minterp(interp_points2).reshape(shape)
202			else:
203			interp_points = np.dstack((np.asarray(sunzsec),
204			np.asarray(180. - azidiff),
205			np.asarray(satzsec)))
206
207			ipn = interpn((sunz_sec_coord, azid_coord, satz_sec_coord),
208			raylwvl[:, :, :], interp_points)
209
210			LOG.debug("Time - Interpolation: {0:f}".format(time.time() - tic))
211
212			ipn *= 100
213			res = ipn
214			if blueband is not None:
215			res = np.where(np.less(blueband, 20.), res,
216			(1 - (blueband - 20) / 80) * res)
217
218			return np.clip(res, 0, 100)
219
220
221			def get_reflectance_lut(filename):
222			"""Read the LUT with reflectances as a function of wavelength, satellite
223			zenith secant, azimuth difference angle, and sun zenith secant
224
225			"""
226
227			with h5py.File(filename, 'r') as h5f:
228			tab = h5f['reflectance'][:]
229			wvl = h5f['wavelengths'][:]
230			azidiff = h5f['azimuth_difference'][:]
231			satellite_zenith_secant = h5f['satellite_zenith_secant'][:]
232			sun_zenith_secant = h5f['sun_zenith_secant'][:]
233
234			return tab, wvl, azidiff, satellite_zenith_secant, sun_zenith_secant
235
236			if __name__ == "__main__":
237
238			this = Rayleigh('Suomi-NPP', 'viirs')
239			# SUNZ = np.arange(200000).reshape(400, 500) * 0.0004
240			# SATZ = np.arange(200000).reshape(400, 500) * 0.00025
241			# AZIDIFF = np.arange(200000).reshape(400, 500) * 0.0009
242			# rfl = this.get_reflectance(SUNZ, SATZ, AZIDIFF, 'M4')
243
244			SHAPE = (1000, 3000)
245			NDIM = SHAPE[0] * SHAPE[1]
246			SUNZ = np.ma.arange(
247			NDIM / 2, NDIM + NDIM / 2).reshape(SHAPE) * 60. / float(NDIM)
248			SATZ = np.ma.arange(NDIM).reshape(SHAPE) * 60. / float(NDIM)
249			AZIDIFF = np.ma.arange(NDIM).reshape(SHAPE) * 179.9 / float(NDIM)
250			rfl = this.get_reflectance(SUNZ, SATZ, AZIDIFF, 'M4')
251

pytroll / pyspectral

Push — develop ( f92b69...710c13 )

get_reflectance_lut() A

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like