Skip to content

Commit

Permalink
Merge branch 'develop' into pre-master
Browse files Browse the repository at this point in the history
  • Loading branch information
Adam.Dybbroe committed Jun 26, 2015
2 parents a002405 + 7b45589 commit 43ff100
Show file tree
Hide file tree
Showing 8 changed files with 310 additions and 17 deletions.
26 changes: 26 additions & 0 deletions etc/pyspectral.cfg_template
Expand Up @@ -49,6 +49,12 @@ tb2rad_lut_filename = /path/to/radiance/tb/lut/data/tb2rad_lut_noaa18_avhrr_ir3.
[NOAA-19-avhrr]
path = /path/to/original/avhrr/rsr/data/
tb2rad_lut_filename = /path/to/radiance/tb/lut/data/tb2rad_lut_noaa19_avhrr_ir3.7.npz
ch1 = NOAA_19_A308C001.txt
ch2 = NOAA_19_A308C002.txt
ch3a = NOAA_19_A308C03A.txt
ch3b = NOAA_19_A308C03B.txt
ch4 = NOAA_19_A308C004.txt
ch5 = NOAA_19_A308C005.txt


[Metop-B-avhrr]
Expand All @@ -62,3 +68,23 @@ path = /path/to/original/modis/rst/data/L1B_RSR_LUT
[EOS-Aqua-modis]
path = /path/to/original/modis/rst/data/aqua


[Himawari-8-ahi]
path = /path/to/original/ahi/data
ch1 = AHI_Sep2013_CH01.SRF
ch2 = AHI_Sep2013_CH02.SRF
ch3 = AHI_Sep2013_CH03.SRF
ch4 = AHI_Sep2013_CH04.SRF
ch5 = AHI_Sep2013_CH05.SRF
ch6 = AHI_Sep2013_CH06.SRF
ch7 = AHI_Sep2013_CH07.SRF
ch8 = AHI_Sep2013_CH08.SRF
ch9 = AHI_Sep2013_CH09.SRF
ch10 = AHI_Sep2013_CH10.SRF
ch11 = AHI_Sep2013_CH11.SRF
ch12 = AHI_Sep2013_CH12.SRF
ch13 = AHI_Sep2013_CH13.SRF
ch14 = AHI_Sep2013_CH14.SRF
ch15 = AHI_Sep2013_CH15.SRF
ch16 = AHI_Sep2013_CH16.SRF

180 changes: 180 additions & 0 deletions pyspectral/ahi_reader.py
@@ -0,0 +1,180 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-

# Copyright (c) 2015 Adam.Dybbroe

# Author(s):

# Adam.Dybbroe <a000680@c14526.ad.smhi.se>

# 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/>.

"""Read the Himawari AHI spectral response functions. Data from
http://cimss.ssec.wisc.edu/goes/calibration/SRF/ahi/
"""

import logging
LOG = logging.getLogger(__name__)

import ConfigParser
import os
import numpy as np

from pyspectral.utils import get_central_wave

try:
CONFIG_FILE = os.environ['PSP_CONFIG_FILE']
except KeyError:
LOG.exception('Environment variable PSP_CONFIG_FILE not set!')
raise

if not os.path.exists(CONFIG_FILE) or not os.path.isfile(CONFIG_FILE):
raise IOError(str(CONFIG_FILE) + " pointed to by the environment " +
"variable PSP_CONFIG_FILE is not a file or does not exist!")

AHI_BAND_NAMES = ['ch1', 'ch2', 'ch3', 'ch4', 'ch5',
'ch6', 'ch7', 'ch8', 'ch9', 'ch10',
'ch11', 'ch12', 'ch13', 'ch14', 'ch15', 'ch16']

#: Default time format
_DEFAULT_TIME_FORMAT = '%Y-%m-%d %H:%M:%S'

#: Default log format
_DEFAULT_LOG_FORMAT = '[%(levelname)s: %(asctime)s : %(name)s] %(message)s'


class AhiRSR(object):

"""Container for the Himawari AHI relative spectral response data"""

def __init__(self, bandname, satname='himawari8'):
"""
"""
self.satname = satname
self.bandname = bandname
self.filenames = {}
self.rsr = None

conf = ConfigParser.ConfigParser()
try:
conf.read(CONFIG_FILE)
except ConfigParser.NoSectionError:
LOG.exception('Failed reading configuration file: ' +
str(CONFIG_FILE))
raise

options = {}
for option, value in conf.items(self.satname + '-ahi', raw=True):
options[option] = value

for option, value in conf.items('general', raw=True):
options[option] = value

self.output_dir = options.get('rsr_dir', './')

self._get_bandfilenames(**options)
LOG.debug("Filenames: " + str(self.filenames))
if os.path.exists(self.filenames[bandname]):
self.requested_band_filename = self.filenames[bandname]
self._load()
else:
raise IOError("Couldn't find an existing file for this band: " +
str(self.bandname))

# To be compatible with VIIRS....
self.filename = self.requested_band_filename

def _get_bandfilenames(self, **options):
"""Get the AHI rsr filenames"""

path = options["path"]
for band in AHI_BAND_NAMES:
LOG.debug("Band = " + str(band))
self.filenames[band] = os.path.join(path, options[band])
LOG.debug(self.filenames[band])
if not os.path.exists(self.filenames[band]):
LOG.warning("Couldn't find an existing file for this band: " +
str(self.filenames[band]))

def _load(self):
"""Load the Himawari AHI RSR data for the band requested
"""

data = np.genfromtxt(self.requested_band_filename,
unpack=True,
names=['wavenumber',
'response'])
# Data seems to be wavenumbers (in some unit...)
# Now, provide wavelengths in micro meters
wavelength = 1.e4 / data['wavenumber']
response = data['response']

idx = np.argsort(wavelength)
wavelength = np.take(wavelength, idx)
response = np.take(response, idx)
self.rsr = {}
self.rsr['det-1'] = {'wavelength': wavelength,
'response': response}


def convert2hdf5(platform_id, sat_number):
"""Retrieve original RSR data and convert to internal hdf5 format"""

import h5py

satellite_id = platform_id + str(sat_number)

ahi = AhiRSR('ch1', satellite_id)
filename = os.path.join(ahi.output_dir,
"rsr_ahi_%s%d.h5" % (platform_id, sat_number))

with h5py.File(filename, "w") as h5f:
h5f.attrs['description'] = 'Relative Spectral Responses for AHI'
h5f.attrs['platform'] = platform_id
h5f.attrs['sat_number'] = sat_number
h5f.attrs['band_names'] = AHI_BAND_NAMES

for chname in AHI_BAND_NAMES:
ahi = AhiRSR(chname, satellite_id)
grp = h5f.create_group(chname)
wvl = ahi.rsr[
'det-1']['wavelength'][~np.isnan(ahi.rsr['det-1']['wavelength'])]
rsp = ahi.rsr[
'det-1']['response'][~np.isnan(ahi.rsr['det-1']['wavelength'])]
grp.attrs['central_wavelength'] = get_central_wave(wvl, rsp)
arr = ahi.rsr['det-1']['wavelength']
dset = grp.create_dataset('wavelength', arr.shape, dtype='f')
dset.attrs['unit'] = 'm'
dset.attrs['scale'] = 1e-06
dset[...] = arr
arr = ahi.rsr['det-1']['response']
dset = grp.create_dataset('response', arr.shape, dtype='f')
dset[...] = arr


if __name__ == "__main__":

import sys
LOG = logging.getLogger('ahi_rsr')
handler = logging.StreamHandler(sys.stderr)

formatter = logging.Formatter(fmt=_DEFAULT_LOG_FORMAT,
datefmt=_DEFAULT_TIME_FORMAT)
handler.setFormatter(formatter)
handler.setLevel(logging.DEBUG)
LOG.setLevel(logging.DEBUG)
LOG.addHandler(handler)

#ahi = AhiRSR('ch1', 'himawari8')
convert2hdf5('himawari', 8)
2 changes: 1 addition & 1 deletion pyspectral/avhrr_rsr.py
Expand Up @@ -128,7 +128,7 @@ def convert2hdf5(platform_name):
h5f.attrs['band_names'] = AVHRR_BAND_NAMES

for chname in AVHRR_BAND_NAMES:
avhrr = AvhrrRSR(chname, 'noaa18')
avhrr = AvhrrRSR(chname, satellite_id)
grp = h5f.create_group(chname)
wvl = avhrr.rsr['wavelength'][~np.isnan(avhrr.rsr['wavelength'])]
rsp = avhrr.rsr['response'][~np.isnan(avhrr.rsr['wavelength'])]
Expand Down
25 changes: 16 additions & 9 deletions pyspectral/near_infrared_reflectance.py
Expand Up @@ -67,13 +67,14 @@ class Calculator(RadTbConverter):
"""

def __init__(self, platform_name, instrument, bandname,
solar_flux=None, **options):
solar_flux=None, **kwargs):
"""Init"""
super(Calculator, self).__init__(platform_name, instrument,
bandname, method=1, **options)
bandname, method=1, **kwargs)

self.bandname = BANDNAMES.get(bandname, bandname)

options = {}
if CONFIG_FILE:
conf = ConfigParser.ConfigParser()
try:
Expand All @@ -82,7 +83,6 @@ def __init__(self, platform_name, instrument, bandname,
LOG.warning('Failed reading configuration file: %s',
str(CONFIG_FILE))

options = {}
for option, value in conf.items(platform_name + '-' + instrument,
raw=True):
options[option] = value
Expand All @@ -96,8 +96,8 @@ def __init__(self, platform_name, instrument, bandname,
self._solar_radiance = None
self._rad39_correction = 1.0

if 'detector' in options:
self.detector = options['detector']
if 'detector' in kwargs:
self.detector = kwargs['detector']
else:
self.detector = 'det-1'

Expand All @@ -107,7 +107,10 @@ def __init__(self, platform_name, instrument, bandname,

if 'tb2rad_lut_filename' in options:
self.lutfile = options['tb2rad_lut_filename']
LOG.info("lut filename: %s", str(self.lutfile))
if not self.lutfile.endswith('.npz'):
self.lutfile = self.lutfile + '.npz'

LOG.info("lut filename: " + str(self.lutfile))
if not os.path.exists(self.lutfile):
self.make_tb2rad_lut(self.bandname,
self.lutfile)
Expand Down Expand Up @@ -197,6 +200,8 @@ def reflectance_from_tbs(self, sun_zenith, tb_near_ir, tb_thermal,
ch37name = '20'
elif self.instrument in ["avhrr", "avhrr3", "avhrr/3"]:
ch37name = 'ch3b'
elif self.instrument == 'ahi':
ch37name = 'ch7'
else:
raise NotImplementedError('Not yet support for this '
'instrument %s' % str(self.instrument))
Expand All @@ -215,10 +220,10 @@ def reflectance_from_tbs(self, sun_zenith, tb_near_ir, tb_thermal,
# FIXME!
scale = self.rsr_integral * 1e-6
retv = self.tb2radiance(tb_therm, ch37name, self.lut)
#print("tb2radiance conversion: " + str(retv))
# print("tb2radiance conversion: " + str(retv))
thermal_emiss_one = retv['radiance'] * scale
retv = self.tb2radiance(tb_nir, ch37name, self.lut)
#print("tb2radiance conversion: " + str(retv))
# print("tb2radiance conversion: " + str(retv))
l_nir = retv['radiance'] * scale

if thermal_emiss_one.ravel().shape[0] < 10:
Expand All @@ -231,7 +236,7 @@ def reflectance_from_tbs(self, sun_zenith, tb_near_ir, tb_thermal,
sunz = sunz.filled(88.)

mu0 = np.cos(np.deg2rad(sunz))
#mu0 = np.where(np.less(mu0, 0.1), 0.1, mu0)
# mu0 = np.where(np.less(mu0, 0.1), 0.1, mu0)
self._rad37 = l_nir
self._rad37_t11 = thermal_emiss_one
self._solar_radiance = 1. / np.pi * self.solar_flux * mu0
Expand All @@ -240,6 +245,8 @@ def reflectance_from_tbs(self, sun_zenith, tb_near_ir, tb_thermal,
# 13.4 micron is provided:
if co2corr:
self.derive_rad39_corr(tb_therm, tbco2)
else:
self._rad39_correction = 1.0

# mask = thermal_emiss_one > l_nir

Expand Down
4 changes: 3 additions & 1 deletion pyspectral/radiance_tb_conversion.py
Expand Up @@ -176,7 +176,9 @@ def tb2radiance(self, tb_, bandname, lut=None):
ntb = (tb_ * self.tb_scale).astype('int16')
start = int(lut['tb'][0] * self.tb_scale)
retv = {}
retv['radiance'] = lut['radiance'][ntb - start]
bounds = 0, lut['radiance'].shape[0] - 1
index = np.clip(ntb - start, bounds[0], bounds[1])
retv['radiance'] = lut['radiance'][index]
retv['unit'] = unit
retv['scale'] = scale
return retv
Expand Down
2 changes: 1 addition & 1 deletion pyspectral/tests/test_reflectance.py
Expand Up @@ -119,7 +119,7 @@ def test_rsr_integral(self):

refl37 = ProductionClassWN(
'EOS-Aqua', 'modis', '20', wavespace='wavenumber')
expected = 130.0039
expected = 130.00391
self.assertAlmostEqual(refl37.rsr_integral, expected, 4)

def test_reflectance(self):
Expand Down
12 changes: 7 additions & 5 deletions pyspectral/tests/test_solarflux.py
@@ -1,7 +1,7 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-

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

# Author(s):

Expand All @@ -22,7 +22,7 @@

"""Unit testing the solar_flux calculations"""

from pyspectral.solar import (SolarIrradianceSpectrum,
from pyspectral.solar import (SolarIrradianceSpectrum,
TOTAL_IRRADIANCE_SPECTRUM_2000ASTM)

import os
Expand All @@ -49,6 +49,7 @@
1., 0.92676, 0.67429, 0.44715, 0.27762, 0.14852,
0.07141, 0.04151, 0.02925, 0.02085, 0.01414, 0.01], dtype='double')


class TestSolarflux(unittest.TestCase):

"""Unit testing the solar flux calculations"""
Expand All @@ -60,9 +61,9 @@ def setUp(self):
return

def test_read(self):
"""Test that solar irradiance spctrum"""
"""Test that solar irradiance spctrum"""
self.solar_irr = \
SolarIrradianceSpectrum(TOTAL_IRRADIANCE_SPECTRUM_2000ASTM,
SolarIrradianceSpectrum(TOTAL_IRRADIANCE_SPECTRUM_2000ASTM,
dlambda=0.005)
self.assertTrue(os.path.exists(self.solar_irr.filename))

Expand All @@ -72,8 +73,9 @@ def test_read(self):
def test_solar_flux(self):
"""Calculate the solar-flux"""
self.solar_irr = \
SolarIrradianceSpectrum(TOTAL_IRRADIANCE_SPECTRUM_2000ASTM,
SolarIrradianceSpectrum(TOTAL_IRRADIANCE_SPECTRUM_2000ASTM,
dlambda=0.005)

# rsr function (se above) is given in micronsm therefore the scale
# factor is 1.0 and not 1e+6 (default)!
sflux = self.solar_irr.inband_solarflux(self.rsr, scale=1.0)
Expand Down

0 comments on commit 43ff100

Please sign in to comment.