using modifier nir_reflectance gives error

[hproe]

A clear and concise description of what the bug is.

Loading SLSTR band S7 with modifier nir_reflectance fails with error: ValueError: operands could not be broadcast together with shapes (2400, 3000) (1200, 1500) To Reproduce

import satpy
satpy.config.set(config_path=['/home/hp/Documents/mySatpyConfigs/'])

from satpy import Scene
from satpy.composites import GenericCompositor
import pyresample
from IPython.display import Image
from pathlib import Path
from glob import glob

dataPath=str(Path().home())+'/DATA/s3/x'
files=glob(dataPath+'/*/*')

scn=Scene(reader='slstr_l1b',filenames=files)
area='pyrenees'
cmp='s7refl'
code=' SLSTR S7refl'
scn.load([cmp])

Expected behavior

code to load/display the reflective part of the band

Actual results

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[2], line 20
     18 cmp='s7refl'
     19 code=' SLSTR S7refl'
---> 20 scn.load([cmp])
     21 scn=scn.resample(area)
     23 prodPath=str(Path().home())+'/y/'

File ~/anaconda3/envs/satpy/lib/python3.10/site-packages/satpy/scene.py:1360, in Scene.load(self, wishlist, calibration, resolution, polarization, level, modifiers, generate, unload, **kwargs)
   1358 self._read_datasets_from_storage(**kwargs)
   1359 if generate:
-> 1360     self.generate_possible_composites(unload)

File ~/anaconda3/envs/satpy/lib/python3.10/site-packages/satpy/scene.py:1423, in Scene.generate_possible_composites(self, unload)
   1416 def generate_possible_composites(self, unload):
   1417     """See which composites can be generated and generate them.
   1418 
   1419     Args:
   1420         unload (bool): if the dependencies of the composites
   1421                        should be unloaded after successful generation.
   1422     """
-> 1423     keepables = self._generate_composites_from_loaded_datasets()
   1425     if self.missing_datasets:
   1426         self._remove_failed_datasets(keepables)

File ~/anaconda3/envs/satpy/lib/python3.10/site-packages/satpy/scene.py:1442, in Scene._generate_composites_from_loaded_datasets(self)
   1439 trunk_nodes = self._dependency_tree.trunk(limit_nodes_to=self.missing_datasets,
   1440                                           limit_children_to=self._datasets.keys())
   1441 needed_comp_nodes = set(self._filter_loaded_datasets_from_trunk_nodes(trunk_nodes))
-> 1442 return self._generate_composites_nodes_from_loaded_datasets(needed_comp_nodes)

File ~/anaconda3/envs/satpy/lib/python3.10/site-packages/satpy/scene.py:1448, in Scene._generate_composites_nodes_from_loaded_datasets(self, compositor_nodes)
   1446 keepables = set()
   1447 for node in compositor_nodes:
-> 1448     self._generate_composite(node, keepables)
   1449 return keepables

File ~/anaconda3/envs/satpy/lib/python3.10/site-packages/satpy/scene.py:1506, in Scene._generate_composite(self, comp_node, keepables)
   1503     return
   1505 try:
-> 1506     composite = compositor(prereq_datasets,
   1507                            optional_datasets=optional_datasets,
   1508                            **comp_node.name.to_dict())
   1509     cid = DataID.new_id_from_dataarray(composite)
   1510     self._datasets[cid] = composite

File ~/anaconda3/envs/satpy/lib/python3.10/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 ~/anaconda3/envs/satpy/lib/python3.10/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 ~/anaconda3/envs/satpy/lib/python3.10/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 ~/anaconda3/envs/satpy/lib/python3.10/site-packages/pyspectral/near_infrared_reflectance.py:275, in Calculator.reflectance_from_tbs(self, sun_zenith, tb_near_ir, tb_thermal, **kwargs)
    273 corrected_thermal_emiss_one = thermal_emiss_one * self._rad3x_correction
    274 nomin = l_nir - corrected_thermal_emiss_one
--> 275 denom = self._solar_radiance - corrected_thermal_emiss_one
    276 data = nomin / denom
    277 mask = denom < EPSILON

File ~/anaconda3/envs/satpy/lib/python3.10/site-packages/dask/array/core.py:223, in check_if_handled_given_other.<locals>.wrapper(self, other)
    221     return NotImplemented
    222 else:
--> 223     return f(self, other)

File ~/anaconda3/envs/satpy/lib/python3.10/site-packages/dask/array/core.py:2387, in Array.__sub__(self, other)
   2385 @check_if_handled_given_other
   2386 def __sub__(self, other):
-> 2387     return elemwise(operator.sub, self, other)

File ~/anaconda3/envs/satpy/lib/python3.10/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 ~/anaconda3/envs/satpy/lib/python3.10/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 (2400, 3000) (1200, 1500)

Screenshots not applicable

Environment Info:

• Linux
• Satpy Version: 0.42.0]
• PyResample Version:

Additional context

composite.yaml:

  s7refl:
    compositor: !!python/name:satpy.composites.SingleBandCompositor
    prerequisites:
      - name: S7
        modifiers: [nir_reflectance]
    standard_name: s7

[djhoese]

As discussed on slack (and this is mostly as context for @mraspaud if he has any other ideas), the primary issue is a but in the NIRReflectance modifier here:

https://github.com/pytroll/satpy/blob/2dba2258b616b6298369581dfeda94f059d62319/satpy/modifiers/spectral.py#L64-L70

where it "matches" DataArrays only for the hard requirements (prerequisites) and does not include the optional ones. As we can see here:

https://github.com/pytroll/satpy/blob/2dba2258b616b6298369581dfeda94f059d62319/satpy/etc/composites/slstr.yaml#L19-L26

this configurationg for slstr has optional inputs. So what is happening is that the resolutions are different and the modifier is supposed to raise a IncompatibleAreas exception to let the Scene know that resampling is required to perform this modification, but this isn't happening.

One suggested solution was for @hproe to specify resolution=1000 in the .load call which should force solar_zenith_angle (available in 500 and 1000 meter resolutions) to be loaded. @hproe said on slack that this does not work. @hproe do you get a different error message or the exact same one? Make sure this resolution is provided to the Scene.load call.

As another test, @hproe you could do:

scn.load(['S8', 'solar_zenith_angle'], resolution=1000)
print(scn['S8'].attrs['_satpy_id'])
print(scn['solar_zenith_angle'].attrs['_satpy_id'])

And let us know the output.

[hproe]

With and without adding resolution=1000 I get the same error message: ValueError: operands could not be broadcast together with shapes (2400, 3000) (1200, 1500) I.e. the proposed code above for printouts does not run, it quits when loading the scene.

[djhoese]

With my above requested test code, did you comment out your existing scn.load call? Meaning, don't try to load anything else, only the load call I have in my code.

[hproe]

Hi Ryan -

Sorry for my not correctly executing your snippet. Here is what I get when using your load command:

DataID(name='S8', wavelength=WavelengthRange(min=10.4, central=10.85, max=11.3, unit='µm'), resolution=1000, calibration=, view=, stripe=, modifiers=()) DataID(name='solar_zenith_angle', resolution=1000, view=, modifiers=())

cheers, HP

On 03.05.23 15:55, David Hoese wrote:

|scn.load(['S8', 'solar_zenith_angle'], resolution=1000) print(scn['S8'].attrs['_satpy_id']) print(scn['solar_zenith_angle'].attrs['_satpy_id'])|

[djhoese]

Thanks. So both of them say they are the same resolution. Could you change the print statements to:

print(scn['S8'].shape, scn['S8'].attrs['area'].shape)
print(scn['solar_zenith_angle'].shape, scn['solar_zenith_angle'].attrs['area'].shape)

And tell me what that outputs?

[hproe]

(1200, 1500) (1200, 1500) (1200, 1500) (2400, 3000)

On 10.05.23 13:58, David Hoese wrote:

|print(scn['S8'].shape, scn['S8'].attrs['area'].shape) print(scn['solar_zenith_angle'].shape, scn['solar_zenith_angle'].attrs['area'].shape) |

[djhoese]

There it is! The area for the solar zenith angle does not match the size of the solar zenith angle data! Now to find out where/why that happens...

[djhoese]

Oh it's a swath. This could be a bigger bug (there is still the bug in the modifier) than we thought. @hproe could you do one more print for me:

print(scn['S8'].attrs['area'].lons.shape)
print(scn['solar_zenith_angle'].attrs['area'].lons.shape)

[hproe]

Looks like there is a discrepncy:

(1200, 1500) (2400, 3000)

On 11.05.23 16:49, David Hoese wrote:

Oh it's a swath. This could be a bigger bug (there is still the bug in the modifier) than we thought. @hproe https://github.com/hproe could you do one more print for me:

|print(scn['S8'].attrs['area'].lons.shape) print(scn['solar_zenith_angle'].attrs['area'].lons.shape) |

— Reply to this email directly, view it on GitHub https://github.com/pytroll/satpy/issues/2460#issuecomment-1544126881, or unsubscribe https://github.com/notifications/unsubscribe-auth/AITAXJL6O2RME3SMHYSC2MDXFT355ANCNFSM6AAAAAAXUAGBIU. You are receiving this because you were mentioned.Message ID: @.***>

[djhoese]

Oops, too much going on this morning. I actually wanted you to check if there was a _satpy_id on the longitude:

print(scn['S8'].attrs['area'].lons.attrs["_satpy_id"])
print(scn['solar_zenith_angle'].attrs['area'].lons.attrs["_satpy_id"])

I can't remember if the lons/lats in the SwathDefinition (the area) are xarray DataArrays or dask arrays so this might fail if they aren't DataArrays.

[hproe]

DataID(name='longitude', resolution=1000, view=, stripe=, modifiers=()) DataID(name='longitude', resolution=1000, view=, stripe=, modifiers=())

On 11.05.23 17:29, David Hoese wrote:

Oops, too much going on this morning. I actually wanted you to check if there was a |_satpy_id| on the longitude:

|print(scn['S8'].attrs['area'].lons.attrs["_satpy_id"]) print(scn['solar_zenith_angle'].attrs['area'].lons.attrs["_satpy_id"]) |

I can't remember if the lons/lats in the SwathDefinition (the area) are xarray DataArrays or dask arrays so this might fail if they aren't DataArrays.

— Reply to this email directly, view it on GitHub https://github.com/pytroll/satpy/issues/2460#issuecomment-1544203778, or unsubscribe https://github.com/notifications/unsubscribe-auth/AITAXJKXO5MRH5GPXFEQCOLXFUAXPANCNFSM6AAAAAAXUAGBIU. You are receiving this because you were mentioned.Message ID: @.***>

[djhoese]

Do you have any idea if the geolocation files contain multiple resolutions of the longitude and latitude arrays? Also, is it a concern that the stripe is different between two datasets?

[hproe]

I'll have to think about it (and tend biomass around the house). Give me some time. cheers, HP

On 11.05.23 17:40, David Hoese wrote:

Do you have any idea if the geolocation files contain multiple resolutions of the longitude and latitude arrays? Also, is it a concern that the |stripe| is different between two datasets?

— Reply to this email directly, view it on GitHub https://github.com/pytroll/satpy/issues/2460#issuecomment-1544221714, or unsubscribe https://github.com/notifications/unsubscribe-auth/AITAXJJ7R6RTZJ5NYP5GQ4DXFUB6FANCNFSM6AAAAAAXUAGBIU. You are receiving this because you were mentioned.Message ID: @.***>

[hproe]

Unfortunately, I am not very much versed in SLSTR. For my purpose up to now I took the data as is.

From what I see in the netCDF files coming with a dataset there is just one resolution - each set for bands S7 and higher (NIR-IR) is 1200 rows by 1500 columns (bands S1 to S6 (VIS) have resolution halved -> 2400 by 3000).

From the COPERNICUS archive (THE source) the datasets are all of the same size. I can send you a sample file (320MB). If yes, pls tell me where to upload it.

cheers,HP

On 11.05.23 17:40, David Hoese wrote:

Do you have any idea if the geolocation files contain multiple resolutions of the longitude and latitude arrays? Also, is it a concern that the |stripe| is different between two datasets?

— Reply to this email directly, view it on GitHub https://github.com/pytroll/satpy/issues/2460#issuecomment-1544221714, or unsubscribe https://github.com/notifications/unsubscribe-auth/AITAXJJ7R6RTZJ5NYP5GQ4DXFUB6FANCNFSM6AAAAAAXUAGBIU. You are receiving this because you were mentioned.Message ID: @.***>

[mraspaud]

Different stripes should not be an issue for channels with same resolution.

[djhoese]

@hproe What do your filenames look like? I've never used sentinel data. I'm trying to download some from copernicus but I have no idea what options to choose.

[hproe]

S3A_SL_1_RBT____20211203T100925_20211203T101225_20211203T124052_0179_079_179_2160_LN2_O_NR_004.SEN3

which actually is a folder containing the (mostly nc) files listed in the attached list. How to order see attached screenshot.

HP

On 12.05.23 17:19, David Hoese wrote:

@hproe https://github.com/hproe What do your filenames look like? I've never used sentinel data. I'm trying to download some from copernicus but I have no idea what options to choose.

— Reply to this email directly, view it on GitHub https://github.com/pytroll/satpy/issues/2460#issuecomment-1545912895, or unsubscribe https://github.com/notifications/unsubscribe-auth/AITAXJLBTE2LEE7GOJQD6GLXFZIJTANCNFSM6AAAAAAXUAGBIU. You are receiving this because you were mentioned.Message ID: @.***>

cartesian_an.nc cartesian_ao.nc cartesian_bn.nc cartesian_bo.nc cartesian_fn.nc cartesian_fo.nc cartesian_in.nc cartesian_io.nc cartesian_tx.nc F1_BT_fn.nc F1_BT_fo.nc F1_quality_fn.nc F1_quality_fo.nc F2_BT_in.nc F2_BT_io.nc F2_quality_in.nc F2_quality_io.nc flags_an.nc flags_ao.nc flags_bn.nc flags_bo.nc flags_fn.nc flags_fo.nc flags_in.nc flags_io.nc geodetic_an.nc geodetic_ao.nc geodetic_bn.nc geodetic_bo.nc geodetic_fn.nc geodetic_fo.nc geodetic_in.nc geodetic_io.nc geodetic_tx.nc geometry_tn.nc geometry_to.nc indices_an.nc indices_ao.nc indices_bn.nc indices_bo.nc indices_fn.nc indices_fo.nc indices_in.nc indices_io.nc list.txt met_tx.nc S1_quality_an.nc S1_quality_ao.nc S1_radiance_an.nc S1_radiance_ao.nc S2_quality_an.nc S2_quality_ao.nc S2_radiance_an.nc S2_radiance_ao.nc S3_quality_an.nc S3_quality_ao.nc S3_radiance_an.nc S3_radiance_ao.nc S4_quality_an.nc S4_quality_ao.nc S4_quality_bn.nc S4_quality_bo.nc S4_radiance_an.nc S4_radiance_ao.nc S4_radiance_bn.nc S4_radiance_bo.nc S5_quality_an.nc S5_quality_ao.nc S5_quality_bn.nc S5_quality_bo.nc S5_radiance_an.nc S5_radiance_ao.nc S5_radiance_bn.nc S5_radiance_bo.nc S6_quality_an.nc S6_quality_ao.nc S6_quality_bn.nc S6_quality_bo.nc S6_radiance_an.nc S6_radiance_ao.nc S6_radiance_bn.nc S6_radiance_bo.nc S7_BT_in.nc S7_BT_io.nc S7_quality_in.nc S7_quality_io.nc S8_BT_in.nc S8_BT_io.nc S8_quality_in.nc S8_quality_io.nc S9_BT_in.nc S9_BT_io.nc S9_quality_in.nc S9_quality_io.nc time_an.nc time_bn.nc time_in.nc viscal.nc xfdumanifest.xml

[hproe]

or ... I can stack some somewhere.

On 12.05.23 17:19, David Hoese wrote:

@hproe https://github.com/hproe What do your filenames look like? I've never used sentinel data. I'm trying to download some from copernicus but I have no idea what options to choose.

— Reply to this email directly, view it on GitHub https://github.com/pytroll/satpy/issues/2460#issuecomment-1545912895, or unsubscribe https://github.com/notifications/unsubscribe-auth/AITAXJLBTE2LEE7GOJQD6GLXFZIJTANCNFSM6AAAAAAXUAGBIU. You are receiving this because you were mentioned.Message ID: @.***>

[djhoese]

@mraspaud It seems the i stripe doesn't have a 500m resolution:

netcdf geodetic_an {
dimensions:
        columns = 3000 ;
        orphan_pixels = 374 ;
        rows = 2187 ;

        int longitude_an(rows, columns) ;
                longitude_an:_FillValue = -2147483648 ;
                longitude_an:add_offset = 0. ;
                longitude_an:long_name = "Longitude of detector FOV centre on the earth\'s surface" ;
                longitude_an:scale_factor = 1.e-06 ;
                longitude_an:standard_name = "longitude" ;
                longitude_an:units = "degrees_east" ;
                longitude_an:valid_max = 180000000 ;
                longitude_an:valid_min = -180000000 ;

netcdf geodetic_in {
dimensions:
        columns = 1500 ;
        orphan_pixels = 187 ;
        rows = 1094 ;
        int longitude_in(rows, columns) ;
                longitude_in:_FillValue = -2147483648 ;
                longitude_in:add_offset = 0. ;
                longitude_in:long_name = "Longitude of detector FOV centre on the earth\'s surface" ;
                longitude_in:scale_factor = 1.e-06 ;
                longitude_in:standard_name = "longitude" ;
                longitude_in:units = "degrees_east" ;
                longitude_in:valid_max = 180000000 ;
                longitude_in:valid_min = -180000000 ;

So it does matter which stripe is used when requesting the lon/lats. If you specify resolution=1000 you get:

In [12]: scn.load(['S8', 'solar_zenith_angle'], resolution=1000)

In [13]: print(scn['S8'].attrs['area'].lons.attrs["_satpy_id"])
    ...: print(scn['solar_zenith_angle'].attrs['area'].lons.attrs["_satpy_id"])
DataID(name='longitude', resolution=1000, view=<view.nadir>, stripe=<stripe.i>, modifiers=())
DataID(name='longitude', resolution=1000, view=<view.nadir>, stripe=<stripe.a>, modifiers=())

In [14]:

In [14]: print(scn['S8'].shape, scn['S8'].attrs['area'].shape)
    ...: print(scn['solar_zenith_angle'].shape, scn['solar_zenith_angle'].attrs['area'].shape)
(1094, 1500) (1094, 1500)
(1094, 1500) (2187, 3000)

If you let it go with the default then you get solar_zenith_angle at 500m:

In [16]: scn.load(['S8', 'solar_zenith_angle'])

In [17]: print(scn['S8'].attrs['area'].lons.attrs["_satpy_id"])
    ...: print(scn['solar_zenith_angle'].attrs['area'].lons.attrs["_satpy_id"])
DataID(name='longitude', resolution=1000, view=<view.nadir>, stripe=<stripe.i>, modifiers=())
DataID(name='longitude', resolution=500, view=<view.nadir>, stripe=<stripe.a>, modifiers=())

In [18]: print(scn['S8'].shape, scn['S8'].attrs['area'].shape)
    ...: print(scn['solar_zenith_angle'].shape, scn['solar_zenith_angle'].attrs['area'].shape)
(1094, 1500) (1094, 1500)
(2187, 3000) (2187, 3000)

Since solar zenith angle doesn't have any inherent stripe to it, I think the default (first) stripe enum value is chosen a.

I think this could be fixed if the combination of a and 500 didn't exist (and possibly other stripes if they don't have that resolution) in the YAML.

[djhoese]

a and b have 500m only, f and i are 1000m only

...unless I'm missing something.

[gerritholl]

The same (?) problem occurs when loading FCI FDHSI + HRFI data:

import hdf5plugin
from satpy import Scene
from glob import glob
from satpy.utils import debug_on; debug_on()
import tempfile
import os
import pathlib

config = """sensor/name: visir

composites:
  dwd_ir38:
    compositor: !!python/name:satpy.composites.SingleBandCompositor
    prerequisites:
      - name: ir_38
        modifiers: [atm_correction]
    standard_name: dwd_ir38

  dwd_ir38_reflectance:
    compositor: !!python/name:satpy.composites.SingleBandCompositor
    prerequisites:
      - name: ir_38
        modifiers: [nir_reflectance]
    standard_name: dwd_ir38_reflectance
"""
with tempfile.TemporaryDirectory() as td:
    p = pathlib.Path(td)
    fn = p / "composites" / "seviri.yaml"
    fn.parent.mkdir(exist_ok=True, parents=True)
    with fn.open(mode="wt", encoding="ascii") as fp:
        fp.write(config)
os.environ["SATPY_CONFIG_PATH"] = td

fci_files = glob(f"/media/nas/x23352/MTG/FCI/2023/12/03/W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-*-FD--CHK-BODY--DIS-NC4E_C_EUMT_20231203*_IDPFI_OPE_20231203*_20231203*_N_JLS_C_0050_*.nc")
sc = Scene(filenames={"fci_l1c_nc": fci_files})
sc.load(["dwd_ir38_reflectance"])

Fails with:

[DEBUG: 2023-12-08 14:54:14 : pyspectral.near_infrared_reflectance] Derive the 3.9 micron radiance CO2 correction coefficent
Traceback (most recent call last):
  File "/data/gholl/checkouts/protocode/mwe/fci-dust-and-ir38-failure.py", line 36, in <module>
    sc.load(["dust", "dwd_ir38_reflectance"])
  File "/data/gholl/checkouts/satpy/satpy/scene.py", line 1445, in load
    self.generate_possible_composites(unload)
  File "/data/gholl/checkouts/satpy/satpy/scene.py", line 1508, in generate_possible_composites
    keepables = self._generate_composites_from_loaded_datasets()
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/data/gholl/checkouts/satpy/satpy/scene.py", line 1527, in _generate_composites_from_loaded_datasets
    return self._generate_composites_nodes_from_loaded_datasets(needed_comp_nodes)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/data/gholl/checkouts/satpy/satpy/scene.py", line 1533, in _generate_composites_nodes_from_loaded_datasets
    self._generate_composite(node, keepables)
  File "/data/gholl/checkouts/satpy/satpy/scene.py", line 1557, in _generate_composite
    prereq_datasets = self._get_prereq_datasets(
                      ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/data/gholl/checkouts/satpy/satpy/scene.py", line 1639, in _get_prereq_datasets
    self._generate_composite(prereq_node, keepables)
  File "/data/gholl/checkouts/satpy/satpy/scene.py", line 1591, in _generate_composite
    composite = compositor(prereq_datasets,
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/data/gholl/checkouts/satpy/satpy/modifiers/spectral.py", line 71, in __call__
    return self._get_reflectance_as_dataarray(*inputs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/data/gholl/checkouts/satpy/satpy/modifiers/spectral.py", line 76, in _get_reflectance_as_dataarray
    reflectance = self._get_reflectance_as_dask(nir.data, da_tb11, da_tb13_4, da_sun_zenith, nir.attrs)
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/data/gholl/checkouts/satpy/satpy/modifiers/spectral.py", line 125, in _get_reflectance_as_dask
    return reflectance_3x_calculator.reflectance_from_tbs(da_sun_zenith, da_nir, da_tb11, tb_ir_co2=da_tb13_4) * 100
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/data/gholl/checkouts/pyspectral/pyspectral/near_infrared_reflectance.py", line 267, in reflectance_from_tbs
    self.derive_rad39_corr(tb_therm, tbco2)
  File "/data/gholl/checkouts/pyspectral/pyspectral/near_infrared_reflectance.py", line 159, in derive_rad39_corr
    self._rad3x_correction = (bt11 - 0.25 * (bt11 - bt13)) ** 4 / bt11 ** 4
                                             ~~~~~^~~~~~
  File "/data/gholl/mambaforge/envs/py312/lib/python3.12/site-packages/dask/array/core.py", line 222, in wrapper
    return f(self, other)
           ^^^^^^^^^^^^^^
  File "/data/gholl/mambaforge/envs/py312/lib/python3.12/site-packages/dask/array/core.py", line 2404, in __sub__
    return elemwise(operator.sub, self, other)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/data/gholl/mambaforge/envs/py312/lib/python3.12/site-packages/dask/array/core.py", line 4787, in elemwise
    broadcast_shapes(*shapes)
  File "/data/gholl/mambaforge/envs/py312/lib/python3.12/site-packages/dask/array/core.py", line 4715, in broadcast_shapes
    raise ValueError(
ValueError: operands could not be broadcast together with shapes (11136, 11136) (5568, 5568)

Loading works fine when using only FDHSI and no HRFI.

I added this as a comment because I think it's the same problem. If it turns out not to be, I can open a new issue.

[pnuu]

Because I haven't seen this with FCI, using scn.load(..., generate=False) should be a valid workaround (it is generally also much faster if the channels are used in multiple composites). I already hunted down the place where a check should be added, but can't now found where I commented it. But with a very fast browse of the code, this line should be wrapped around try/except ValueError: raise IncompatibleAreas.

[gerritholl]

The generate=False workaround works for FCI. I have not tested SLSTR.

[djhoese]

@pnuu why is that try/except needed? What is the ValueError that is being caught and why isn't it caught by match_data_arrays above that line?

[djhoese]

Ooohhh optional_datasets isn't included in the projectables but it should be. That's the "proper" fix for this I think.

[pnuu]

@djhoese my recollection from the previous inspection was that the match_data_arrays() didn't raise, but it seems that part is further down in check_geolocation(). And the ValueError would've been the one in the trace Gerrit posted about the shapes.

[djhoese]

Right. The match_data_arrays is supposed to check all shape and area compatibility. The problem here is that the optionals are not included in that check and should be. In that case match_data_arrays would have raised IncompatibleAreas and the "derive" function never would have been called.

[djhoese]

Kind of related: #2694