VIIRS day_microphysics array broadcast error

[martin-rdz]

Describe the bug
Generating the day_microphysics RGB with NOAA-20 VIIRS data fails due to an array broadcast error.

To Reproduce

fls = glob(f'data_laads_netcdf/NOAA-20/20221029_1456/*1454*')
scene = Scene(filenames=fls, reader='viirs_l1b')
scene.load(['day_microphysics'], generate=True)
scene.save_dataset('day_microphysics', 
                   filename='day_microphysics.png')

Expected behavior
Successful generation of the day_microphysics RGB

Actual results
A value error occurs when calculating the reflectance_from_tbs.

[DEBUG: 2022-11-03 10:44:39 : satpy.readers.yaml_reader] Reading ('/dev/2022-10-05_polar_orbiting_satellites/satproc-venv/lib/python3.8/site-packages/satpy/etc/readers/viirs_l1b.yaml',)
[DEBUG: 2022-11-03 10:44:40 : satpy.readers.yaml_reader] Assigning to viirs_l1b: ['data_laads_netcdf/NOAA-20/20221029_1456/VJ103IMG_NRT.A2022302.1454.021.2022302162534.nc', 'data_laads_netcdf/NOAA-20/20221029_1456/VJ103MOD_NRT.A2022302.1454.021.2022302162534.nc', 'data_laads_netcdf/NOAA-20/20221029_1456/VJ103DNB_NRT.A2022302.1454.021.2022302162534.nc', 'data_laads_netcdf/NOAA-20/20221029_1456/VJ102IMG_NRT.A2022302.1454.021.2022302164936.nc', 'data_laads_netcdf/NOAA-20/20221029_1456/VJ102MOD_NRT.A2022302.1454.021.2022302164936.nc', 'data_laads_netcdf/NOAA-20/20221029_1456/VJ102DNB_NRT.A2022302.1454.021.2022302164936.nc']
[DEBUG: 2022-11-03 10:44:40 : satpy.node] Skipping optional DataQuery(wavelength=13.4, calibration='brightness_temperature'): Unknown dataset DataQuery(wavelength=13.4, calibration='brightness_temperature')
[DEBUG: 2022-11-03 10:44:40 : satpy.readers.yaml_reader] No coordinates found for DataID(name='m_lon', resolution=742, modifiers=())
[DEBUG: 2022-11-03 10:44:40 : satpy.readers.yaml_reader] No coordinates found for DataID(name='i_lat', resolution=371, modifiers=())
[DEBUG: 2022-11-03 10:44:40 : satpy.readers.yaml_reader] No coordinates found for DataID(name='i_lon', resolution=371, modifiers=())
[DEBUG: 2022-11-03 10:44:40 : satpy.readers.yaml_reader] No coordinates found for DataID(name='m_lat', resolution=742, modifiers=())
[DEBUG: 2022-11-03 10:44:40 : satpy.readers.viirs_l1b] Adjusting scaling factors to convert '1' to '%'
[DEBUG: 2022-11-03 10:44:40 : satpy.readers.viirs_l1b] File units and output units are the same (degrees)
[INFO: 2022-11-03 10:44:40 : satpy.modifiers.spectral] Getting reflective part of I04
[DEBUG: 2022-11-03 10:44:40 : pyspectral.rsr_reader] Filename: .local/share/pyspectral/rsr_viirs_NOAA-20.h5
[DEBUG: 2022-11-03 10:44:40 : pyspectral.rsr_reader] Filename: .local/share/pyspectral/rsr_viirs_NOAA-20.h5
[DEBUG: 2022-11-03 10:44:40 : pyspectral.solar] Begin and end wavelength/wavenumber: 1.000004 7.133945 
[INFO: 2022-11-03 10:44:40 : pyspectral.near_infrared_reflectance] No lut filename available in config file. Will generate filename automatically
[INFO: 2022-11-03 10:44:40 : pyspectral.near_infrared_reflectance] lut filename: /tmp/tb2rad_lut_noaa-20_viirs_i4.npz
[INFO: 2022-11-03 10:44:40 : pyspectral.near_infrared_reflectance] File was there and has been read!
[DEBUG: 2022-11-03 10:44:40 : pyspectral.near_infrared_reflectance] Apply sun-zenith angle clipping between 0 and 88.00

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In [28], line 14
     12 fls = glob(f'data_laads_netcdf/{sat}/20221029_1456/*1454*')
     13 scene = Scene(filenames=fls, reader='viirs_l1b')
---> 14 scene.load(['day_microphysics'], generate=True)
     15 scene.save_dataset('day_microphysics', 
     16                    filename='day_microphysics.png')

File /dev/2022-10-05_polar_orbiting_satellites/satproc-venv/lib/python3.8/site-packages/satpy/scene.py:1352, in Scene.load(self, wishlist, calibration, resolution, polarization, level, generate, unload, **kwargs)
   1350 self._read_datasets_from_storage(**kwargs)
   1351 if generate:
-> 1352     self.generate_possible_composites(unload)

File /dev/2022-10-05_polar_orbiting_satellites/satproc-venv/lib/python3.8/site-packages/satpy/scene.py:1415, in Scene.generate_possible_composites(self, unload)
   1408 def generate_possible_composites(self, unload):
   1409     """See which composites can be generated and generate them.
   1410 
   1411     Args:
   1412         unload (bool): if the dependencies of the composites
   1413                        should be unloaded after successful generation.
   1414     """
-> 1415     keepables = self._generate_composites_from_loaded_datasets()
   1417     if self.missing_datasets:
   1418         self._remove_failed_datasets(keepables)

File /dev/2022-10-05_polar_orbiting_satellites/satproc-venv/lib/python3.8/site-packages/satpy/scene.py:1434, in Scene._generate_composites_from_loaded_datasets(self)
   1431 trunk_nodes = self._dependency_tree.trunk(limit_nodes_to=self.missing_datasets,
   1432                                           limit_children_to=self._datasets.keys())
   1433 needed_comp_nodes = set(self._filter_loaded_datasets_from_trunk_nodes(trunk_nodes))
-> 1434 return self._generate_composites_nodes_from_loaded_datasets(needed_comp_nodes)

File /dev/2022-10-05_polar_orbiting_satellites/satproc-venv/lib/python3.8/site-packages/satpy/scene.py:1440, in Scene._generate_composites_nodes_from_loaded_datasets(self, compositor_nodes)
   1438 keepables = set()
   1439 for node in compositor_nodes:
-> 1440     self._generate_composite(node, keepables)
   1441 return keepables

File /dev/2022-10-05_polar_orbiting_satellites/satproc-venv/lib/python3.8/site-packages/satpy/scene.py:1498, in Scene._generate_composite(self, comp_node, keepables)
   1495     return
   1497 try:
-> 1498     composite = compositor(prereq_datasets,
   1499                            optional_datasets=optional_datasets,
   1500                            **comp_node.name.to_dict())
   1501     cid = DataID.new_id_from_dataarray(composite)
   1502     self._datasets[cid] = composite

File /dev/2022-10-05_polar_orbiting_satellites/satproc-venv/lib/python3.8/site-packages/satpy/modifiers/spectral.py:70, in NIRReflectance.__call__(self, projectables, optional_datasets, **info)
     65 """Get the reflectance part of an NIR channel.
     66 
     67 Not supposed to be used for wavelength outside [3, 4] µm.
     68 """
     69 projectables = self.match_data_arrays(projectables)
---> 70 return self._get_reflectance_as_dataarray(projectables, optional_datasets)

File /dev/2022-10-05_polar_orbiting_satellites/satproc-venv/lib/python3.8/site-packages/satpy/modifiers/spectral.py:81, in NIRReflectance._get_reflectance_as_dataarray(self, projectables, optional_datasets)
     78 da_sun_zenith = self._get_sun_zenith_from_provided_data(projectables, optional_datasets)
     80 logger.info('Getting reflective part of %s', _nir.attrs['name'])
---> 81 reflectance = self._get_reflectance_as_dask(da_nir, da_tb11, da_tb13_4, da_sun_zenith, _nir.attrs)
     83 proj = self._create_modified_dataarray(reflectance, base_dataarray=_nir)
     84 proj.attrs['units'] = '%'

File /dev/2022-10-05_polar_orbiting_satellites/satproc-venv/lib/python3.8/site-packages/satpy/modifiers/spectral.py:125, in NIRReflectance._get_reflectance_as_dask(self, da_nir, da_tb11, da_tb13_4, da_sun_zenith, metadata)
    123 """Calculate 3.x reflectance in % with pyspectral from dask arrays."""
    124 reflectance_3x_calculator = self._init_reflectance_calculator(metadata)
--> 125 return reflectance_3x_calculator.reflectance_from_tbs(da_sun_zenith, da_nir, da_tb11, tb_ir_co2=da_tb13_4) * 100

File /dev/2022-10-05_polar_orbiting_satellites/satproc-venv/lib/python3.8/site-packages/pyspectral/near_infrared_reflectance.py:279, in Calculator.reflectance_from_tbs(self, sun_zenith, tb_near_ir, tb_thermal, **kwargs)
    277 corrected_thermal_emiss_one = thermal_emiss_one * self._rad3x_correction
    278 nomin = l_nir - corrected_thermal_emiss_one
--> 279 denom = self._solar_radiance - corrected_thermal_emiss_one
    280 data = nomin / denom
    281 mask = denom < EPSILON

File /dev/2022-10-05_polar_orbiting_satellites/satproc-venv/lib/python3.8/site-packages/dask/array/core.py:224, in check_if_handled_given_other.<locals>.wrapper(self, other)
    222     return NotImplemented
    223 else:
--> 224     return f(self, other)

File /dev/2022-10-05_polar_orbiting_satellites/satproc-venv/lib/python3.8/site-packages/dask/array/core.py:2384, in Array.__sub__(self, other)
   2382 @check_if_handled_given_other
   2383 def __sub__(self, other):
-> 2384     return elemwise(operator.sub, self, other)

File /dev/2022-10-05_polar_orbiting_satellites/satproc-venv/lib/python3.8/site-packages/dask/array/core.py:4762, in elemwise(op, out, where, dtype, name, *args, **kwargs)
   4758     shapes.append(out.shape)
   4760 shapes = [s if isinstance(s, Iterable) else () for s in shapes]
   4761 out_ndim = len(
-> 4762     broadcast_shapes(*shapes)
   4763 )  # Raises ValueError if dimensions mismatch
   4764 expr_inds = tuple(range(out_ndim))[::-1]
   4766 if dtype is not None:

File /dev/2022-10-05_polar_orbiting_satellites/satproc-venv/lib/python3.8/site-packages/dask/array/core.py:4690, in broadcast_shapes(*shapes)
   4688         dim = 0 if 0 in sizes else np.max(sizes)
   4689     if any(i not in [-1, 0, 1, dim] and not np.isnan(i) for i in sizes):
-> 4690         raise ValueError(
   4691             "operands could not be broadcast together with "
   4692             "shapes {}".format(" ".join(map(str, shapes)))
   4693         )
   4694     out.append(dim)
   4695 return tuple(reversed(out))

ValueError: operands could not be broadcast together with shapes (3232, 3200) (6464, 6400)

Environment Info:

• OS: WSL ubuntu 20.04
• Satpy Version: 0.37.1
• PyResample Version: 1.25.1

Additional context
The data files were obtained via the LAADS DAAC NRT server.

[pnuu]

The day_microphysics composite is defined in visir.yaml using wavelengths. With multi-resolution data, this often causes the data to be from a mixture of different resolution channels. The options you have:

• do not set compute=True and resample the data
• supply Scene only either IMG or MOD data
• define the composite using channel names (see example for M channels below) in your $SATPY_CONFIG_PATH/composites/viirs.yaml

  day_microphysics:
    compositor: !!python/name:satpy.composites.GenericCompositor
    prerequisites:
    - name: M07
      modifiers: [sunz_corrected]
    - name: M12
      modifiers: [nir_reflectance]
    - M15
    standard_name: day_microphysics

[pnuu]

And references on adding your own configuration and composites:

• https://satpy.readthedocs.io/en/stable/config.html
• https://satpy.readthedocs.io/en/stable/composites.html#creating-composite-configuration-files

[martin-rdz]

The definition using the explicit channel names solved the issue.

For the sake of completeness, here my solution for the high resolution composite

  day_microphysics_hr:
    compositor: !!python/name:satpy.composites.GenericCompositor
    prerequisites:
    - name: I02
      modifiers: [sunz_corrected_iband]
    - name: I04
      modifiers: [nir_reflectance_hires]
    - I05
    standard_name: day_microphysics

One caveat with the LAADS DAAC data (reader='viirs_l1b') is, that the solar zenith angle is named differently and causes
an error:

KeyError: "Unknown datasets: {DataQuery(name='I02', modifiers=('sunz_corrected_iband',)): 
{DataQuery(name='solar_zenith_angle', resolution=371, calibration='reflectance'): 
{DataQuery(name='solar_zenith_angle', resolution=371, calibration='reflectance')}}}"

Following issue #1190, this works (also in viirs.yaml):

  sunz_corrected_iband:
    modifier: !!python/name:satpy.modifiers.SunZenithCorrector
    prerequisites:
    - name: i_solar_zenith_angle
      resolution: 371

Thanks a lot for the quick and helpful reply.

[djhoese]

I'm going to reopen this. In my opinion there should be a better error message for this. @pnuu is right that changing what resolutions you provide should get around this issue, but generate=True is the default so you didn't do anything wrong in your code. There should have been a more obvious warning telling you to resample the data.

[djhoese]

I'm able to reproduce this error. Interestingly I don't get this error with VIIRS SDR data, only VIIRS L1B. This suggests it may be a reader issue.

[djhoese]

I've created #2260 to fix some of the issues involved here. I had missed in the #1190 issue that changes to the reader were used. My PR adds those changes which now fix the error you were getting originally @martin-rdz. I think this means it is still doing more processing than necessary (trying to do sunz correction with M band angles on an I band data), but it at least doesn't error out. I think your high resolution version and even a low resolution version of this composite that are specific to VIIRS would be a good addition to Satpy. How would you feel about make a pull request with these additions?