planck() - Code Metrics - Inspection of "Clean up for easier maintanance" - pytroll/pyspectral - Measure and Improve Code Quality continuously with Scrutinizer

Test Failed

Push — develop ( 43e3bd...6a143d )

by Adam

created 2017-07-14 15:07 UTC

planck() F

↳ Parent: Project

Complexity

Conditions

Size

Total Lines

Duplication

Lines	0
Ratio	0 %

Importance

Changes	1
Bugs	0	Features	0

Metric	Value
cc	11
c	1
b	0
f	0
dl	0
loc	82
rs	3.1764

How to fix Long Method Complexity

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

# Copyright (c) 2013, 2014, 2015, 2016, 2017 Adam.Dybbroe

# Author(s):

#   Adam.Dybbroe <[email protected]>
#   Panu Lahtinen <[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/>.

"""Planck radiation equation"""
import numpy as np

import logging
LOG = logging.getLogger(__name__)

H_PLANCK = 6.62606957 * 1e-34  # SI-unit = [J*s]
K_BOLTZMANN = 1.3806488 * 1e-23  # SI-unit = [J/K]
C_SPEED = 2.99792458 * 1e8  # SI-unit = [m/s]

EPSILON = 0.000001


def blackbody_rad2temp(wavelength, radiance):
    """Derive brightness temperatures from radiance using the Planck
    function. Wavelength space. Assumes SI units as input and returns
    temperature in Kelvin

    """

    mask = False
    if np.isscalar(radiance):
        rad = np.array([radiance, ], dtype='float64')
    else:
        rad = np.array(radiance, dtype='float64')
        if np.ma.is_masked(radiance):
            mask = radiance.mask

    rad = np.ma.masked_array(rad, mask=mask)
    rad = np.ma.masked_less_equal(rad, 0)

    if np.isscalar(wavelength):
        wvl = np.array([wavelength, ], dtype='float64')
    else:
        wvl = np.array(wavelength, dtype='float64')

    const1 = H_PLANCK * C_SPEED / K_BOLTZMANN
    const2 = 2 * H_PLANCK * C_SPEED**2
    res = const1 / (wvl * np.log(np.divide(const2, rad * wvl**5) + 1.0))

    shape = rad.shape
    resshape = res.shape

    if wvl.shape[0] == 1:
        if rad.shape[0] == 1:
            return res[0]
        else:
            return res[::].reshape(shape)
    else:
        if rad.shape[0] == 1:
            return res[0, :]
        else:
            if len(shape) == 1:
                return np.reshape(res, (shape[0], resshape[1]))
            else:
                return np.reshape(res, (shape[0], shape[1], resshape[1]))


def blackbody_wn_rad2temp(wavenumber, radiance):
    """Derive brightness temperatures from radiance using the Planck 
    function. Wavenumber space"""

    if np.isscalar(radiance):
        rad = np.array([radiance, ], dtype='float64')
    else:
        rad = np.array(radiance, dtype='float64')
    if np.isscalar(wavenumber):
        wavnum = np.array([wavenumber, ], dtype='float64')
    else:
        wavnum = np.array(wavenumber, dtype='float64')

    const1 = H_PLANCK * C_SPEED / K_BOLTZMANN
    const2 = 2 * H_PLANCK * C_SPEED**2
    res = const1 * wavnum / np.log(np.divide(const2 * wavnum**3, rad) + 1.0)

    shape = rad.shape
    resshape = res.shape

    if wavnum.shape[0] == 1:

        if rad.shape[0] == 1:
            return res[0]
        else:
            return res[::].reshape(shape)
    else:
        if rad.shape[0] == 1:
            return res[0, :]
        else:
            if len(shape) == 1:
                return np.reshape(res, (shape[0], resshape[1]))
            else:
                return np.reshape(res, (shape[0], shape[1], resshape[1]))


def planck(wave, temp, wavelength=True):
    """The Planck radiation or Blackbody radiation as a function of wavelength 
    or wavenumber. SI units.
    _planck(wave, temperature, wavelength=True)
    wave = Wavelength/wavenumber or a sequence of wavelengths/wavenumbers (m or m^-1)
    temp = Temperature (scalar) or a sequence of temperatures (K)


    Output: Wavelength space: The spectral radiance per meter (not micron!)
            Unit = W/m^2 sr^-1 m^-1

            Wavenumber space: The spectral radiance in Watts per square meter
            per steradian per m-1:
            Unit = W/m^2 sr^-1 (m^-1)^-1 = W/m sr^-1

            Converting from SI units to mW/m^2 sr^-1 (cm^-1)^-1:
            1.0 W/m^2 sr^-1 (m^-1)^-1 = 0.1 mW/m^2 sr^-1 (cm^-1)^-1

    """
    units = ['wavelengths', 'wavenumbers']
    if wavelength:
        LOG.debug("Using {0} when calculating the Blackbody radiance".format(
            units[(wavelength == True) - 1]))

    if np.isscalar(temp):
        temperature = np.array([temp, ], dtype='float64')
    else:
        temperature = np.array(temp, dtype='float64')

    shape = temperature.shape
    if np.isscalar(wave):
        wln = np.array([wave, ], dtype='float64')
    else:
        wln = np.array(wave, dtype='float64')

    if wavelength:
        const = 2 * H_PLANCK * C_SPEED ** 2
        nom = const / wln ** 5
        arg1 = H_PLANCK * C_SPEED / (K_BOLTZMANN * wln)
    else:
        nom = 2 * H_PLANCK * (C_SPEED ** 2) * (wln ** 3)
        arg1 = H_PLANCK * C_SPEED * wln / K_BOLTZMANN

    arg2 = np.where(np.greater(np.abs(temperature), EPSILON),
                    np.array(1. / temperature), -9).reshape(-1, 1)
    arg2 = np.ma.masked_array(arg2, mask=arg2 == -9)
    LOG.debug("Max and min - arg1: %s  %s", str(arg1.max()), str(arg1.min()))
    LOG.debug("Max and min - arg2: %s  %s", str(arg2.max()), str(arg2.min()))
    try:
        exp_arg = np.multiply(arg1.astype('float32'), arg2.astype('float32'))
    except MemoryError:
        LOG.warning(("Dimensions used in numpy.multiply probably reached "
                     "limit!\n"
                     "Make sure the Radiance<->Tb table has been created "
                     "and try running again"))
        raise

    LOG.debug("Max and min before exp: %s  %s", str(exp_arg.max()),
              str(exp_arg.min()))
    if exp_arg.min() < 0:
        LOG.warning("Something is fishy: \n" +
                    "\tDenominator might be zero or negative in radiance derivation:")
        dubious = np.where(exp_arg < 0)[0]
        LOG.warning(
            "Number of items having dubious values: " + str(dubious.shape[0]))

    denom = np.exp(exp_arg) - 1
    rad = nom / denom
    radshape = rad.shape
    if wln.shape[0] == 1:
        if temperature.shape[0] == 1:
            return rad[0, 0]
        else:
            return rad[:, 0].reshape(shape)
    else:
        if temperature.shape[0] == 1:
            return rad[0, :]
        else:
            if len(shape) == 1:
                return np.reshape(rad, (shape[0], radshape[1]))
            else:
                return np.reshape(rad, (shape[0], shape[1], radshape[1]))


def blackbody_wn(wavenumber, temp):
    """The Planck radiation or Blackbody radiation as a function of wavenumber
    SI units!
    blackbody_wn(wavnum, temperature)
    wavenumber = A wavenumber (scalar) or a sequence of wave numbers (m-1)
    temp = A temperatfure (scalar) or a sequence of temperatures (K)


    Output: The spectral radiance in Watts per square meter per steradian
            per m-1:
            Unit = W/m^2 sr^-1 (m^-1)^-1 = W/m sr^-1

            Converting from SI units to mW/m^2 sr^-1 (cm^-1)^-1:
            1.0 W/m^2 sr^-1 (m^-1)^-1 = 0.1 mW/m^2 sr^-1 (cm^-1)^-1
    """

    return planck(wavenumber, temp, wavelength=False)


def blackbody(wavel, temp):
    """The Planck radiation or Blackbody radiation as a function of wavelength
    SI units.
    blackbody(wavelength, temperature)
    wavel = Wavelength or a sequence of wavelengths (m)
    temp = Temperature (scalar) or a sequence of temperatures (K)


    Output: The spectral radiance per meter (not micron!)
            Unit = W/m^2 sr^-1 m^-1

    """

    return planck(wavel, temp, wavelength=True)


1		#!/usr/bin/env python
2		# -- coding: utf-8 --
3
4		# Copyright (c) 2013, 2014, 2015, 2016, 2017 Adam.Dybbroe
5
6		# Author(s):
7
8		# Adam.Dybbroe <[email protected]>
9		# Panu Lahtinen <[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		"""Planck radiation equation"""
25		import numpy as np
26
27		import logging
28		LOG = logging.getLogger(__name__)
29
30		H_PLANCK = 6.62606957 * 1e-34 # SI-unit = [J*s]
31		K_BOLTZMANN = 1.3806488 * 1e-23 # SI-unit = [J/K]
32		C_SPEED = 2.99792458 * 1e8 # SI-unit = [m/s]
33
34		EPSILON = 0.000001
35
36
37		def blackbody_rad2temp(wavelength, radiance):
38		"""Derive brightness temperatures from radiance using the Planck
39		function. Wavelength space. Assumes SI units as input and returns
40		temperature in Kelvin
41
42		"""
43
44		mask = False
45		if np.isscalar(radiance):
46		rad = np.array([radiance, ], dtype='float64')
47		else:
48		rad = np.array(radiance, dtype='float64')
49		if np.ma.is_masked(radiance):
50		mask = radiance.mask
51
52		rad = np.ma.masked_array(rad, mask=mask)
53		rad = np.ma.masked_less_equal(rad, 0)
54
55		if np.isscalar(wavelength):
56		wvl = np.array([wavelength, ], dtype='float64')
57		else:
58		wvl = np.array(wavelength, dtype='float64')
59
60		const1 = H_PLANCK * C_SPEED / K_BOLTZMANN
61		const2 = 2 * H_PLANCK * C_SPEED**2
62		res = const1 / (wvl * np.log(np.divide(const2, rad * wvl**5) + 1.0))
63
64		shape = rad.shape
65		resshape = res.shape
66
67		if wvl.shape[0] == 1:
68		if rad.shape[0] == 1:
69		return res[0]
70		else:
71		return res[::].reshape(shape)
72		else:
73		if rad.shape[0] == 1:
74		return res[0, :]
75		else:
76		if len(shape) == 1:
77		return np.reshape(res, (shape[0], resshape[1]))
78		else:
79		return np.reshape(res, (shape[0], shape[1], resshape[1]))
80
81
82		def blackbody_wn_rad2temp(wavenumber, radiance):
83		"""Derive brightness temperatures from radiance using the Planck
84		function. Wavenumber space"""
85
86		if np.isscalar(radiance):
87		rad = np.array([radiance, ], dtype='float64')
88		else:
89		rad = np.array(radiance, dtype='float64')
90		if np.isscalar(wavenumber):
91		wavnum = np.array([wavenumber, ], dtype='float64')
92		else:
93		wavnum = np.array(wavenumber, dtype='float64')
94
95		const1 = H_PLANCK * C_SPEED / K_BOLTZMANN
96		const2 = 2 * H_PLANCK * C_SPEED**2
97		res = const1 * wavnum / np.log(np.divide(const2 * wavnum**3, rad) + 1.0)
98
99		shape = rad.shape
100		resshape = res.shape
101
102	View Code Duplication	if wavnum.shape[0] == 1:
		0 ignored issues – show Duplication introduced 2017-07-10 09:33 UTC by Report Bug Copy Issue Report This code seems to be duplicated in your project. Loading history...
103		if rad.shape[0] == 1:
104		return res[0]
105		else:
106		return res[::].reshape(shape)
107		else:
108		if rad.shape[0] == 1:
109		return res[0, :]
110		else:
111		if len(shape) == 1:
112		return np.reshape(res, (shape[0], resshape[1]))
113		else:
114		return np.reshape(res, (shape[0], shape[1], resshape[1]))
115
116
117		def planck(wave, temp, wavelength=True):
118		"""The Planck radiation or Blackbody radiation as a function of wavelength
119		or wavenumber. SI units.
120		_planck(wave, temperature, wavelength=True)
121		wave = Wavelength/wavenumber or a sequence of wavelengths/wavenumbers (m or m^-1)
122		temp = Temperature (scalar) or a sequence of temperatures (K)
123
124
125		Output: Wavelength space: The spectral radiance per meter (not micron!)
126		Unit = W/m^2 sr^-1 m^-1
127
128		Wavenumber space: The spectral radiance in Watts per square meter
129		per steradian per m-1:
130		Unit = W/m^2 sr^-1 (m^-1)^-1 = W/m sr^-1
131
132		Converting from SI units to mW/m^2 sr^-1 (cm^-1)^-1:
133		1.0 W/m^2 sr^-1 (m^-1)^-1 = 0.1 mW/m^2 sr^-1 (cm^-1)^-1
134
135		"""
136		units = ['wavelengths', 'wavenumbers']
137		if wavelength:
138		LOG.debug("Using {0} when calculating the Blackbody radiance".format(
139		units[(wavelength == True) - 1]))
140
141		if np.isscalar(temp):
142		temperature = np.array([temp, ], dtype='float64')
143		else:
144		temperature = np.array(temp, dtype='float64')
145
146		shape = temperature.shape
147		if np.isscalar(wave):
148		wln = np.array([wave, ], dtype='float64')
149		else:
150		wln = np.array(wave, dtype='float64')
151
152		if wavelength:
153		const = 2 * H_PLANCK * C_SPEED ** 2
154		nom = const / wln ** 5
155		arg1 = H_PLANCK * C_SPEED / (K_BOLTZMANN * wln)
156		else:
157		nom = 2 * H_PLANCK * (C_SPEED ** 2) * (wln ** 3)
158		arg1 = H_PLANCK * C_SPEED * wln / K_BOLTZMANN
159
160		arg2 = np.where(np.greater(np.abs(temperature), EPSILON),
161		np.array(1. / temperature), -9).reshape(-1, 1)
162		arg2 = np.ma.masked_array(arg2, mask=arg2 == -9)
163		LOG.debug("Max and min - arg1: %s %s", str(arg1.max()), str(arg1.min()))
164		LOG.debug("Max and min - arg2: %s %s", str(arg2.max()), str(arg2.min()))
165		try:
166		exp_arg = np.multiply(arg1.astype('float32'), arg2.astype('float32'))
167		except MemoryError:
168		LOG.warning(("Dimensions used in numpy.multiply probably reached "
169		"limit!\n"
170		"Make sure the Radiance<->Tb table has been created "
171		"and try running again"))
172		raise
173
174		LOG.debug("Max and min before exp: %s %s", str(exp_arg.max()),
175		str(exp_arg.min()))
176		if exp_arg.min() < 0:
177		LOG.warning("Something is fishy: \n" +
178		"\tDenominator might be zero or negative in radiance derivation:")
179		dubious = np.where(exp_arg < 0)[0]
180		LOG.warning(
181		"Number of items having dubious values: " + str(dubious.shape[0]))
182
183		denom = np.exp(exp_arg) - 1
184		rad = nom / denom
185		radshape = rad.shape
186		if wln.shape[0] == 1:
187		if temperature.shape[0] == 1:
188		return rad[0, 0]
189		else:
190		return rad[:, 0].reshape(shape)
191		else:
192		if temperature.shape[0] == 1:
193		return rad[0, :]
194		else:
195		if len(shape) == 1:
196		return np.reshape(rad, (shape[0], radshape[1]))
197		else:
198		return np.reshape(rad, (shape[0], shape[1], radshape[1]))
199
200
201		def blackbody_wn(wavenumber, temp):
202		"""The Planck radiation or Blackbody radiation as a function of wavenumber
203		SI units!
204		blackbody_wn(wavnum, temperature)
205		wavenumber = A wavenumber (scalar) or a sequence of wave numbers (m-1)
206		temp = A temperatfure (scalar) or a sequence of temperatures (K)
207
208
209		Output: The spectral radiance in Watts per square meter per steradian
210		per m-1:
211		Unit = W/m^2 sr^-1 (m^-1)^-1 = W/m sr^-1
212
213		Converting from SI units to mW/m^2 sr^-1 (cm^-1)^-1:
214		1.0 W/m^2 sr^-1 (m^-1)^-1 = 0.1 mW/m^2 sr^-1 (cm^-1)^-1
215		"""
216
217		return planck(wavenumber, temp, wavelength=False)
218
219
220		def blackbody(wavel, temp):
221		"""The Planck radiation or Blackbody radiation as a function of wavelength
222		SI units.
223		blackbody(wavelength, temperature)
224		wavel = Wavelength or a sequence of wavelengths (m)
225		temp = Temperature (scalar) or a sequence of temperatures (K)
226
227
228		Output: The spectral radiance per meter (not micron!)
229	View Code Duplication	Unit = W/m^2 sr^-1 m^-1
		0 ignored issues – show Duplication introduced 2017-07-14 15:11 UTC by Report Bug Copy Issue Report This code seems to be duplicated in your project. Loading history...
230		"""
231
232		return planck(wavel, temp, wavelength=True)
233

pytroll / pyspectral

Push — develop ( 43e3bd...6a143d )

planck() F

Complexity

Size

Duplication

Importance

How to fix Long Method Complexity

Long Method

Complexity

Duplication Side-by-Side

Filter issues like