
This notebook is developed during the pythia cook-off at NCAR Mesa-Lab Boulder Colorado, June 12-14, 2024
Participants in the workshop event have the chance to practice collaborative problem-solving and hands-on learning in the field of Python programming.
This notebook is part of the Breakout Topic: Geostationary on AWS, lead by Jorge Humberto Bravo Mendez jbravo2@stevens.edu, from Stevens Institute of Technology
Advanced Himawari Imager (AHI) data with Satpy¶
Using Satpy to read and Advanced Baseline Imager (ABI) data from GOES-R satellites. Here’s a step-by-step guide:
Imports¶
import warnings
warnings.filterwarnings('ignore')
warnings.simplefilter('ignore', SyntaxWarning)
from satpy.scene import Scene
from satpy.utils import debug_on
from datetime import datetime
from glob import globStarting to create satpy scenes¶
sat_files = glob("input/HI9_AHI_L1b_JP/*")
sat_files['input/HI9_AHI_L1b_JP/HS_H09_20230815_0430_B10_JP02_R20_S0101.DAT.bz2',
'input/HI9_AHI_L1b_JP/HS_H09_20230815_0430_B06_JP02_R20_S0101.DAT.bz2',
'input/HI9_AHI_L1b_JP/HS_H09_20230815_0430_B08_JP02_R20_S0101.DAT.bz2',
'input/HI9_AHI_L1b_JP/HS_H09_20230815_0430_B13_JP02_R20_S0101.DAT.bz2',
'input/HI9_AHI_L1b_JP/HS_H09_20230815_0430_B14_JP02_R20_S0101.DAT.bz2',
'input/HI9_AHI_L1b_JP/HS_H09_20230815_0430_B03_JP02_R05_S0101.DAT.bz2',
'input/HI9_AHI_L1b_JP/HS_H09_20230815_0430_B09_JP02_R20_S0101.DAT.bz2',
'input/HI9_AHI_L1b_JP/HS_H09_20230815_0430_B07_JP02_R20_S0101.DAT.bz2',
'input/HI9_AHI_L1b_JP/HS_H09_20230815_0430_B04_JP02_R10_S0101.DAT.bz2',
'input/HI9_AHI_L1b_JP/HS_H09_20230815_0430_B05_JP02_R20_S0101.DAT.bz2',
'input/HI9_AHI_L1b_JP/HS_H09_20230815_0430_B02_JP02_R10_S0101.DAT.bz2',
'input/HI9_AHI_L1b_JP/HS_H09_20230815_0430_B15_JP02_R20_S0101.DAT.bz2',
'input/HI9_AHI_L1b_JP/HS_H09_20230815_0430_B12_JP02_R20_S0101.DAT.bz2',
'input/HI9_AHI_L1b_JP/HS_H09_20230815_0430_B11_JP02_R20_S0101.DAT.bz2',
'input/HI9_AHI_L1b_JP/HS_H09_20230815_0430_B01_JP02_R10_S0101.DAT.bz2',
'input/HI9_AHI_L1b_JP/HS_H09_20230815_0430_B16_JP02_R20_S0101.DAT.bz2']scn = Scene(filenames = sat_files, reader='ahi_hsd')
dataset_names = scn.all_dataset_names()
print(dataset_names)['B01', 'B02', 'B03', 'B04', 'B05', 'B06', 'B07', 'B08', 'B09', 'B10', 'B11', 'B12', 'B13', 'B14', 'B15', 'B16']
scn.load(dataset_names)print(scn.available_composite_names())['airmass', 'ash', 'cloud_phase_distinction', 'cloud_phase_distinction_raw', 'cloudtop', 'colorized_ir_clouds', 'convection', 'day_microphysics_ahi', 'day_microphysics_eum', 'day_severe_storms', 'dust', 'fire_temperature', 'fire_temperature_39refl', 'fire_temperature_awips', 'fire_temperature_eumetsat', 'fog', 'geo_color', 'geo_color_background_with_low_clouds', 'geo_color_high_clouds', 'geo_color_low_clouds', 'geo_color_night', 'hybrid_green', 'hybrid_green_nocorr', 'ir_cloud_day', 'mid_vapor', 'natural_color', 'natural_color_nocorr', 'natural_color_raw', 'natural_color_raw_with_night_ir', 'ndvi_hybrid_green', 'night_ir_alpha', 'night_ir_with_background', 'night_ir_with_background_hires', 'night_microphysics', 'night_microphysics_tropical', 'overview', 'overview_raw', 'reproduced_green', 'reproduced_green_uncorr', 'rocket_plume_night', 'true_color', 'true_color_ndvi_green', 'true_color_nocorr', 'true_color_raw', 'true_color_reproduction', 'true_color_reproduction_corr', 'true_color_reproduction_night_ir', 'true_color_reproduction_uncorr', 'true_color_with_night_ir', 'true_color_with_night_ir_hires', 'water_vapors1', 'water_vapors2']
rgb_im = 'airmass'
scn.load([rgb_im])### Uncomment to show it
#scn.show(rgb_im)result = scn[rgb_im]
resultLoading...
keys = scn.keys()
keys[DataID(name='B01', wavelength=WavelengthRange(min=0.45, central=0.47, max=0.49, unit='µm'), resolution=1000, calibration=<1>, modifiers=()),
DataID(name='B02', wavelength=WavelengthRange(min=0.49, central=0.51, max=0.53, unit='µm'), resolution=1000, calibration=<1>, modifiers=()),
DataID(name='B03', wavelength=WavelengthRange(min=0.62, central=0.64, max=0.66, unit='µm'), resolution=500, calibration=<1>, modifiers=()),
DataID(name='B04', wavelength=WavelengthRange(min=0.85, central=0.86, max=0.87, unit='µm'), resolution=1000, calibration=<1>, modifiers=()),
DataID(name='B05', wavelength=WavelengthRange(min=1.5, central=1.6, max=1.7, unit='µm'), resolution=2000, calibration=<1>, modifiers=()),
DataID(name='B06', wavelength=WavelengthRange(min=2.2, central=2.3, max=2.4, unit='µm'), resolution=2000, calibration=<1>, modifiers=()),
DataID(name='B07', wavelength=WavelengthRange(min=3.7, central=3.9, max=4.1, unit='µm'), resolution=2000, calibration=<2>, modifiers=()),
DataID(name='B08', wavelength=WavelengthRange(min=6.0, central=6.2, max=6.4, unit='µm'), resolution=2000, calibration=<2>, modifiers=()),
DataID(name='B09', wavelength=WavelengthRange(min=6.7, central=6.9, max=7.1, unit='µm'), resolution=2000, calibration=<2>, modifiers=()),
DataID(name='B10', wavelength=WavelengthRange(min=7.1, central=7.3, max=7.5, unit='µm'), resolution=2000, calibration=<2>, modifiers=()),
DataID(name='B11', wavelength=WavelengthRange(min=8.4, central=8.6, max=8.8, unit='µm'), resolution=2000, calibration=<2>, modifiers=()),
DataID(name='B12', wavelength=WavelengthRange(min=9.4, central=9.6, max=9.8, unit='µm'), resolution=2000, calibration=<2>, modifiers=()),
DataID(name='B13', wavelength=WavelengthRange(min=10.2, central=10.4, max=10.6, unit='µm'), resolution=2000, calibration=<2>, modifiers=()),
DataID(name='B14', wavelength=WavelengthRange(min=11.0, central=11.2, max=11.4, unit='µm'), resolution=2000, calibration=<2>, modifiers=()),
DataID(name='B15', wavelength=WavelengthRange(min=12.2, central=12.4, max=12.6, unit='µm'), resolution=2000, calibration=<2>, modifiers=()),
DataID(name='B16', wavelength=WavelengthRange(min=13.1, central=13.3, max=13.5, unit='µm'), resolution=2000, calibration=<2>, modifiers=()),
DataID(name='airmass', resolution=2000)]area_info = scn["B13"].area
area_infoLoading...
scn.load(["natural_color"])The following datasets were not created and may require resampling to be generated: DataID(name='natural_color')
rs = scn["B13"].area
lscn = scn.resample(rs)lscn.load(["natural_color"])
### Uncomment to show it
#lscn.show("natural_color")lscn.load(['true_color'])
### Uncomment to show it
#lscn.show('true_color') 0%| | 0/52 [00:00<?, ?kB/s]100%|██████████| 52/52 [00:00<00:00, 1414.83kB/s]
No rsr file /home/runner/.local/share/pyspectral/rsr_ahi_Himawari-9.h5 on disk
0kB [00:00, ?kB/s]1kB [00:00, 22310.13kB/s]
---------------------------------------------------------------------------
ReadError Traceback (most recent call last)
Cell In[14], line 1
----> 1 lscn.load(['true_color'])
3 ### Uncomment to show it
4 #lscn.show('true_color')
File ~/micromamba/envs/geostationary-cookbook/lib/python3.14/site-packages/satpy/scene.py:1484, in Scene.load(self, wishlist, calibration, resolution, polarization, level, modifiers, generate, unload, **kwargs)
1482 self._read_datasets_from_storage(**kwargs)
1483 if generate:
-> 1484 self.generate_possible_composites(unload)
File ~/micromamba/envs/geostationary-cookbook/lib/python3.14/site-packages/satpy/scene.py:1547, in Scene.generate_possible_composites(self, unload)
1540 def generate_possible_composites(self, unload):
1541 """See which composites can be generated and generate them.
1542
1543 Args:
1544 unload (bool): if the dependencies of the composites
1545 should be unloaded after successful generation.
1546 """
-> 1547 keepables = self._generate_composites_from_loaded_datasets()
1549 if self.missing_datasets:
1550 self._remove_failed_datasets(keepables)
File ~/micromamba/envs/geostationary-cookbook/lib/python3.14/site-packages/satpy/scene.py:1566, in Scene._generate_composites_from_loaded_datasets(self)
1563 trunk_nodes = self._dependency_tree.trunk(limit_nodes_to=self.missing_datasets,
1564 limit_children_to=self._datasets.keys())
1565 needed_comp_nodes = set(self._filter_loaded_datasets_from_trunk_nodes(trunk_nodes))
-> 1566 return self._generate_composites_nodes_from_loaded_datasets(needed_comp_nodes)
File ~/micromamba/envs/geostationary-cookbook/lib/python3.14/site-packages/satpy/scene.py:1572, in Scene._generate_composites_nodes_from_loaded_datasets(self, compositor_nodes)
1570 keepables = set()
1571 for node in compositor_nodes:
-> 1572 self._generate_composite(node, keepables)
1573 return keepables
File ~/micromamba/envs/geostationary-cookbook/lib/python3.14/site-packages/satpy/scene.py:1630, in Scene._generate_composite(self, comp_node, keepables)
1627 return
1629 try:
-> 1630 composite = compositor(prereq_datasets,
1631 optional_datasets=optional_datasets,
1632 **comp_node.name.to_dict())
1633 cid = DataID.new_id_from_dataarray(composite)
1634 self._datasets[cid] = composite
File ~/micromamba/envs/geostationary-cookbook/lib/python3.14/site-packages/satpy/modifiers/atmosphere.py:112, in PSPRayleighReflectance.__call__(self, projectables, optional_datasets, **info)
107 corrector = Rayleigh(vis.attrs["platform_name"], vis.attrs["sensor"],
108 atmosphere=atmosphere,
109 aerosol_type=aerosol_type)
111 try:
--> 112 refl_cor_band = corrector.get_reflectance(sunz, satz, ssadiff,
113 vis.attrs["name"],
114 red.data)
115 except (KeyError, IOError):
116 logger.warning("Could not get the reflectance correction using band name: %s", vis.attrs["name"])
File ~/micromamba/envs/geostationary-cookbook/lib/python3.14/site-packages/pyspectral/rayleigh.py:252, in Rayleigh.get_reflectance(self, sun_zenith, sat_zenith, azidiff, band_name_or_wavelength, redband)
248 # if the user gave us non-dask arrays we'll use the dask
249 # version of the algorithm but return numpy arrays back
250 compute = da is not None and not isinstance(sun_zenith, da.Array)
--> 252 wvl, band_name = self._get_effective_wavelength_and_band_name(band_name_or_wavelength)
253 repr_arr = sun_zenith if redband is None else redband
254 try:
File ~/micromamba/envs/geostationary-cookbook/lib/python3.14/site-packages/pyspectral/rayleigh.py:208, in Rayleigh._get_effective_wavelength_and_band_name(self, band_name_or_wavelength)
206 else:
207 band_name = band_name_or_wavelength
--> 208 wvl = _get_rsr_wavelength_from_band_name(
209 self.platform_name,
210 self.sensor,
211 band_name,
212 )
213 band_name_map = BANDNAMES.get(self.sensor, BANDNAMES['generic'])
214 band_name_or_wavelength = band_name_map.get(band_name, band_name)
File ~/micromamba/envs/geostationary-cookbook/lib/python3.14/site-packages/pyspectral/rayleigh.py:306, in _get_rsr_wavelength_from_band_name(platform_name, sensor, band_name)
304 def _get_rsr_wavelength_from_band_name(platform_name, sensor, band_name):
305 try:
--> 306 rsr = RelativeSpectralResponse(platform_name, sensor)
307 except OSError:
308 LOG.exception(
309 "No spectral responses for this platform and sensor: %s %s", platform_name, sensor)
File ~/micromamba/envs/geostationary-cookbook/lib/python3.14/site-packages/pyspectral/rsr_reader.py:136, in RelativeSpectralResponse.__init__(self, platform_name, instrument, filename)
133 self.instrument = instrument
134 self.filename = filename
--> 136 self._check_consistent_input()
138 self.load()
139 self.rsr.instrument = self.instrument
File ~/micromamba/envs/geostationary-cookbook/lib/python3.14/site-packages/pyspectral/rsr_reader.py:149, in RelativeSpectralResponse._check_consistent_input(self)
145 raise AttributeError(
146 "Either platform name and sensor, or filename, must be specified")
148 self._check_instrument()
--> 149 self._get_filename()
151 self._check_filename_exist()
File ~/micromamba/envs/geostationary-cookbook/lib/python3.14/site-packages/pyspectral/rsr_reader.py:173, in RelativeSpectralResponse._get_filename(self)
171 if self.do_download:
172 LOG.info("Will download from internet...")
--> 173 download_rsr()
174 if self._get_rsr_data_version() == RSR_DATA_VERSION:
175 self._rsr_data_version_uptodate = True
File ~/micromamba/envs/geostationary-cookbook/lib/python3.14/site-packages/pyspectral/utils.py:361, in download_rsr(dest_dir, dry_run)
358 if dry_run:
359 return
--> 361 _download_tarball_and_extract(HTTP_PYSPECTRAL_RSR, filename, dest_dir)
File ~/micromamba/envs/geostationary-cookbook/lib/python3.14/site-packages/pyspectral/utils.py:420, in _download_tarball_and_extract(tarball_url, local_pathname, extract_dir)
415 for data in _tqdm_or_iter(response.iter_content(chunk_size=chunk_size),
416 total=(int(total_size / chunk_size + 0.5)),
417 unit='kB'):
418 handle.write(data)
--> 420 tar = tarfile.open(local_pathname)
421 tar_kwargs = {} if sys.version_info < (3, 12) else {"filter": "data"}
422 tar.extractall(extract_dir, **tar_kwargs)
File ~/micromamba/envs/geostationary-cookbook/lib/python3.14/tarfile.py:1904, in TarFile.open(cls, name, mode, fileobj, bufsize, **kwargs)
1902 continue
1903 error_msgs_summary = '\n'.join(error_msgs)
-> 1904 raise ReadError(f"file could not be opened successfully:\n{error_msgs_summary}")
1906 elif ":" in mode:
1907 filemode, comptype = mode.split(":", 1)
ReadError: file could not be opened successfully:
- method gz: ReadError('not a gzip file')
- method bz2: ReadError('not a bzip2 file')
- method xz: ReadError('not an lzma file')
- method zst: ReadError('not a zstd file')
- method tar: ReadError('truncated header')