Skip to content

Commit

Permalink
Merge pull request #145 from simonrp84/arctica-m-srf
Browse files Browse the repository at this point in the history
Add converter for Arctica-M N1 RSR data.
  • Loading branch information
adybbroe committed Oct 10, 2022
2 parents 1cd41c9 + c45a914 commit 3a0f406
Show file tree
Hide file tree
Showing 5 changed files with 184 additions and 10 deletions.
2 changes: 2 additions & 0 deletions pyspectral/etc/pyspectral.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -146,6 +146,8 @@ download_from_internet: True
# Sentinel-3A-olci:
# path: /path/to/original/sentinel-3a/olci/data

# Arctica-M-N1-msu-gsa:
# path: /path/to/original/Arctica_M_N1_SRF.xlsx

# Sentinel-2A-msi:
# path: /path/to/original/sentinel-2a/msi/data/S2-SRF_COPE-GSEG-EOPG-TN-15-0007_3.0.xlsx
Expand Down
39 changes: 39 additions & 0 deletions pyspectral/tests/test_utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@
import pytest
import unittest
import numpy as np
import tempfile
import responses


Expand All @@ -40,6 +41,7 @@


TEST_RSR = {'20': {}, }
TEST_RSR2 = {'20': {}, }
TEST_RSR['20']['det-1'] = {}
TEST_RSR['20']['det-1']['central_wavelength'] = 3.75
TEST_RSR['20']['det-1']['wavelength'] = np.array([
Expand All @@ -60,6 +62,8 @@
1., 0.92676, 0.67429, 0.44715, 0.27762, 0.14852,
0.07141, 0.04151, 0.02925, 0.02085, 0.01414, 0.01], dtype='float32')

TEST_RSR2['20'] = TEST_RSR['20']['det-1']

RESULT_RSR = {'20': {}}
RESULT_RSR['20']['det-1'] = {}
RESULT_RSR['20']['det-1']['wavenumber'] = np.array([
Expand All @@ -84,6 +88,20 @@
dtype='float32')


class _MockedRSR:
def __init__(self, usedet=True):
self.instrument = 'test_sensor'
if usedet:
self.tmp_rsr = TEST_RSR
else:
self.tmp_rsr = TEST_RSR2
self.output_dir = tempfile.mkdtemp()

def __call__(self, chan, plat):
self.rsr = self.tmp_rsr[chan]
return self


class RsrTestData(object):
"""Container for the RSR test datasets."""

Expand Down Expand Up @@ -136,6 +154,27 @@ def test_convert2wavenumber(self):
wvn = newrsr['20']['det-1']['wavenumber']
self.assertTrue(np.allclose(wvn_res, wvn))

def test_convert2hdf5(self):
"""Test the conversion utility from original RSR to HDF5."""
mocked_rsr_nodet = _MockedRSR(usedet=False)
utils.convert2hdf5(mocked_rsr_nodet, 'Test_SAT', ['20'])
fname = f'{mocked_rsr_nodet.output_dir}/rsr_test_sensor_Test_SAT.h5'
self.assertTrue(os.path.exists(fname))
try:
os.remove(fname)
except OSError:
pass

mocked_rsr_det = _MockedRSR(usedet=True)
utils.convert2hdf5(mocked_rsr_det, 'Test_SAT', ['20'], detectors=['det-1'])
fname = f'{mocked_rsr_det.output_dir}/rsr_test_sensor_Test_SAT.h5'
self.assertTrue(os.path.exists(fname))
try:
os.remove(fname)
except OSError:
pass


def test_get_bandname_from_wavelength(self):
"""Test the right bandname is found provided the wavelength in micro meters."""
bname = utils.get_bandname_from_wavelength('abi', 0.4, self.rsr.rsr)
Expand Down
40 changes: 31 additions & 9 deletions pyspectral/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -72,7 +72,8 @@
'Meteosat-8': 'seviri',
'FY-4A': 'agri',
'GEO-KOMPSAT-2A': 'ami',
'MTG-I1': 'fci'
'MTG-I1': 'fci',
'Arctica-M-N1': 'msu-gsa'
}

HTTP_PYSPECTRAL_RSR = "https://zenodo.org/record/7116653/files/pyspectral_rsr_data.tgz"
Expand Down Expand Up @@ -240,7 +241,7 @@ def sort_data(x_vals, y_vals):
return x_vals, y_vals


def convert2hdf5(ClassIn, platform_name, bandnames, scale=1e-06):
def convert2hdf5(ClassIn, platform_name, bandnames, scale=1e-06, detectors=None):
"""Retrieve original RSR data and convert to internal hdf5 format.
*scale* is the number which has to be multiplied to the wavelength data in
Expand All @@ -264,17 +265,38 @@ def convert2hdf5(ClassIn, platform_name, bandnames, scale=1e-06):
for chname in bandnames:
sensor = ClassIn(chname, platform_name)
grp = h5f.create_group(chname)
wvl = sensor.rsr['wavelength'][~np.isnan(sensor.rsr['wavelength'])]
rsp = sensor.rsr['response'][~np.isnan(sensor.rsr['wavelength'])]
grp.attrs['central_wavelength'] = get_central_wave(wvl, rsp)
arr = sensor.rsr['wavelength']

# If multiple detectors, assume all have same wavelength range in SRF.
if detectors is not None:
wvl = sensor.rsr[detectors[0]]['wavelength'][~np.isnan(sensor.rsr[detectors[0]]['wavelength'])]
arr = sensor.rsr[detectors[0]]['wavelength']
grp.attrs['number_of_detectors'] = len(detectors)
else:
wvl = sensor.rsr['wavelength'][~np.isnan(sensor.rsr['wavelength'])]
arr = sensor.rsr['wavelength']

# Save wavelengths to file
dset = grp.create_dataset('wavelength', arr.shape, dtype='f')
dset.attrs['unit'] = 'm'
dset.attrs['scale'] = scale
dset[...] = arr
arr = sensor.rsr['response']
dset = grp.create_dataset('response', arr.shape, dtype='f')
dset[...] = arr

# Now to do the responses
if detectors is None:
rsp = sensor.rsr['response'][~np.isnan(sensor.rsr['wavelength'])]
grp.attrs['central_wavelength'] = get_central_wave(wvl, rsp)
arr = sensor.rsr['response']
dset = grp.create_dataset('response', arr.shape, dtype='f')
dset[...] = arr
else:
for cur_det in detectors:
det_grp = grp.create_group(cur_det)
rsp = sensor.rsr[cur_det]['response'][~np.isnan(sensor.rsr[cur_det]['wavelength'])]
det_grp.attrs['central_wavelength'] = get_central_wave(wvl, rsp)
arr = sensor.rsr[cur_det]['response']
dset = det_grp.create_dataset('response', arr.shape, dtype='f')
dset[...] = arr



def download_rsr(dest_dir=None, dry_run=False):
Expand Down
7 changes: 6 additions & 1 deletion rsr_convert_scripts/README.rst
Original file line number Diff line number Diff line change
Expand Up @@ -174,7 +174,12 @@ the pyspectral.yaml file:
.. code::
%> python msu_gs_reader.py
%> python msu_gsa_reader.py
Converts RSRs for the MSU-GS/A sensors aboard Arctica-M N1 satellite.
RSRs were retrieved from Roshydromet. Filenames look like:

``rtcoef_electro-l_2_msugs_srf_ch01.txt``

Converts RSRs for the MSU-GS sensor aboard Electro-L N2. RSRs were retrieved from the NWP-SAF.
Filenames look like:
Expand Down
106 changes: 106 additions & 0 deletions rsr_convert_scripts/msu_gsa_reader.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,106 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# Copyright (c) 2022 Pytroll developers
#
# 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 Arctica-M N1 MSU-GS/A spectral response functions.
Data from Roshydromet via personal communication.
"""

from pyspectral.utils import convert2hdf5 as tohdf5
from pyspectral.raw_reader import InstrumentRSR
from pyspectral.config import get_config
import pandas as pd
import logging

LOG = logging.getLogger(__name__)

# Names of the individual bands
MSUGSA_BAND_NAMES = ['ch1', 'ch2', 'ch3', 'ch4', 'ch5',
'ch6', 'ch7', 'ch8', 'ch9', 'ch10']

# Set up VIS and IR bands, needed for selecting sheet in RSR file
VISBANDS = {'ch1', 'ch2', 'ch3'}
IRBANDS = {'ch4', 'ch5', 'ch6', 'ch7', 'ch8', 'ch9', 'ch10'}

#: 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 MsugsaRSR(InstrumentRSR):
"""Container for the Arctica-M-N1 MSU-GS/A relative spectral response data"""

def __init__(self, bandname, platform_name):

super(MsugsaRSR, self).__init__(
bandname, platform_name, MSUGSA_BAND_NAMES)

self.instrument = 'msu-gsa'
self.platform_name = platform_name

options = get_config()

self.msugsa_path = options[
self.platform_name + '-' + self.instrument].get('path')

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

self._load()

def _load(self, scale=10000.0):
"""Load the MSU-GS/A RSR data for the band requested"""
detectors = {}
# The Arctica satellites have two instruments on them. Pyspectral isn't set up
# to handle this, so instead we call them separate detectors.
for detnum in (1, 2):
if self.bandname in VISBANDS:
data = pd.read_excel(self.msugsa_path, sheet_name=f'MSU-{detnum} vis', skiprows=2,
names=['wvl', 'ch1', 'ch2', 'ch3'])
elif self.bandname in IRBANDS:
data = pd.read_excel(self.msugsa_path, sheet_name=f'MSU-{detnum} IR', skiprows=2,
names=['wvl', 'ch4', 'ch5', 'ch6', 'ch7', 'ch8', 'ch9', 'ch10'])
else:
raise KeyError(f"{self.bandname} is not a valid VIS or IR band!")
wavelength = data['wvl']
response = data[self.bandname]

detectors[f'det-{detnum}'] = {'wavelength': wavelength, 'response': response}

self.rsr = detectors


def main():
"""Main"""
for platform_name in ['Arctica-M-N1', ]:
tohdf5(MsugsaRSR, platform_name, MSUGSA_BAND_NAMES, detectors=['det-1', 'det-2'])


if __name__ == "__main__":
import sys
LOG = logging.getLogger('msu_gsa_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)

main()

0 comments on commit 3a0f406

Please sign in to comment.