pytroll
using modifier nir_reflectance gives error
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
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.
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.
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.
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=<calibration.brightness_temperature>, view=<view.nadir>, stripe=<stripe.i>, modifiers=()) DataID(name='solar_zenith_angle', resolution=1000, view=<view.nadir>, 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'])|
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?
(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) |
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...
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)
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: @.***>
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.
DataID(name='longitude', resolution=1000, view=<view.nadir>, stripe=<stripe.i>, modifiers=()) DataID(name='longitude', resolution=1000, view=<view.nadir>, stripe=<stripe.a>, 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: @.***>
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?
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: @.***>
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: @.***>
Different stripes should not be an issue for channels with same resolution.
@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.
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
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: @.***>
@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.
a and b have 500m only, f and i are 1000m only
...unless I'm missing something.
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.
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.
@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?
Ooohhh optional_datasets isn't included in the projectables but it should be. That's the "proper" fix for this I think.
@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.
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.