Merge pull request #2642 from pnuu/preserve-nir-reflectance-dtype

Set dtype for get_lonlats() in NIR reflectance calculation

.github/workflows/ci.yaml

@@ -64,27 +64,36 @@ jobs:
       - name: Install unstable dependencies
         if: matrix.experimental == true
         shell: bash -l {0}
+        # Install pykdtree with --no-build-isolation so it builds with numpy 2.0
         # We must get LD_PRELOAD for stdlibc++ or else the manylinux wheels
         # may break the conda-forge libraries trying to use newer glibc versions
         run: |
+          python -m pip install versioneer extension-helpers setuptools-scm configobj pkgconfig
           python -m pip install \
           --index-url https://pypi.anaconda.org/scientific-python-nightly-wheels/simple/ \
           --trusted-host pypi.anaconda.org \
           --no-deps --pre --upgrade \
           matplotlib \
           numpy \
           pandas \
-          scipy; \
-          python -m pip install \
-          --no-deps --upgrade \
+          scipy
+          mamba remove --force-remove -y pykdtree pyresample trollimage pyhdf netcdf4 h5py
+          python -m pip install --upgrade --no-deps --pre --no-build-isolation \
+          git+https://github.com/storpipfugl/pykdtree \
+          git+https://github.com/pytroll/pyresample \
+          git+https://github.com/pytroll/trollimage \
+          git+https://github.com/fhs/pyhdf \
+          git+https://github.com/takluyver/h5py@cython-3 \
+          git+https://github.com/Unidata/netcdf4-python \
           git+https://github.com/dask/dask \
           git+https://github.com/dask/distributed \
           git+https://github.com/zarr-developers/zarr \
           git+https://github.com/Unidata/cftime \
           git+https://github.com/rasterio/rasterio \
           git+https://github.com/pydata/bottleneck \
           git+https://github.com/pydata/xarray \
-          git+https://github.com/astropy/astropy;
+          git+https://github.com/shapely/shapely \
+          git+https://github.com/astropy/astropy
           LD_PRELOAD=$(python -c "import sys; print(sys.prefix)")/lib/libstdc++.so
           echo "LD_PRELOAD=${LD_PRELOAD}" >> $GITHUB_ENV
 

satpy/modifiers/spectral.py

@@ -67,23 +67,24 @@ def __call__(self, projectables, optional_datasets=None, **info):
         Not supposed to be used for wavelength outside [3, 4] µm.
         """
         projectables = self.match_data_arrays(projectables)
-        return self._get_reflectance_as_dataarray(projectables, optional_datasets)
+        inputs = self._get_nir_inputs(projectables, optional_datasets)
+        return self._get_reflectance_as_dataarray(*inputs)
 
-    def _get_reflectance_as_dataarray(self, projectables, optional_datasets):
+    def _get_reflectance_as_dataarray(self, nir, da_tb11, da_tb13_4, da_sun_zenith):
         """Get the reflectance as a dataarray."""
-        _nir, _tb11 = projectables
-        da_nir = _nir.data
-        da_tb11 = _tb11.data
-        da_tb13_4 = self._get_tb13_4_from_optionals(optional_datasets)
-        da_sun_zenith = self._get_sun_zenith_from_provided_data(projectables, optional_datasets)
-
-        logger.info("Getting reflective part of %s", _nir.attrs["name"])
-        reflectance = self._get_reflectance_as_dask(da_nir, da_tb11, da_tb13_4, da_sun_zenith, _nir.attrs)
-
-        proj = self._create_modified_dataarray(reflectance, base_dataarray=_nir)
+        logger.info("Getting reflective part of %s", nir.attrs["name"])
+        reflectance = self._get_reflectance_as_dask(nir.data, da_tb11, da_tb13_4, da_sun_zenith, nir.attrs)
+        proj = self._create_modified_dataarray(reflectance, base_dataarray=nir)
         proj.attrs["units"] = "%"
         return proj
 
+    def _get_nir_inputs(self, projectables, optional_datasets):
+        nir, tb11 = projectables
+        da_tb11 = tb11.data
+        da_tb13_4 = self._get_tb13_4_from_optionals(optional_datasets)
+        da_sun_zenith = self._get_sun_zenith_from_provided_data(nir, optional_datasets, nir.dtype)
+        return (nir, da_tb11, da_tb13_4, da_sun_zenith)
+
     @staticmethod
     def _get_tb13_4_from_optionals(optional_datasets):
         tb13_4 = None
@@ -95,7 +96,7 @@ def _get_tb13_4_from_optionals(optional_datasets):
         return tb13_4
 
     @staticmethod
-    def _get_sun_zenith_from_provided_data(projectables, optional_datasets):
+    def _get_sun_zenith_from_provided_data(nir, optional_datasets, dtype):
         """Get the sunz from available data or compute it if unavailable."""
         sun_zenith = None
 
@@ -106,9 +107,8 @@ def _get_sun_zenith_from_provided_data(projectables, optional_datasets):
         if sun_zenith is None:
             if sun_zenith_angle is None:
                 raise ImportError("Module pyorbital.astronomy needed to compute sun zenith angles.")
-            _nir = projectables[0]
-            lons, lats = _nir.attrs["area"].get_lonlats(chunks=_nir.data.chunks)
-            sun_zenith = sun_zenith_angle(_nir.attrs["start_time"], lons, lats)
+            lons, lats = nir.attrs["area"].get_lonlats(chunks=nir.data.chunks, dtype=dtype)
+            sun_zenith = sun_zenith_angle(nir.attrs["start_time"], lons, lats)
         return sun_zenith
 
     def _create_modified_dataarray(self, reflectance, base_dataarray):
@@ -159,20 +159,16 @@ def __call__(self, projectables, optional_datasets=None, **info):
 
         """
         projectables = self.match_data_arrays(projectables)
-        return self._get_emissivity_as_dataarray(projectables, optional_datasets)
+        inputs = self._get_nir_inputs(projectables, optional_datasets)
+        return self._get_emissivity_as_dataarray(*inputs)
 
-    def _get_emissivity_as_dataarray(self, projectables, optional_datasets):
+    def _get_emissivity_as_dataarray(self, nir, da_tb11, da_tb13_4, da_sun_zenith):
         """Get the emissivity as a dataarray."""
-        _nir, _tb11 = projectables
-        da_nir = _nir.data
-        da_tb11 = _tb11.data
-        da_tb13_4 = self._get_tb13_4_from_optionals(optional_datasets)
-        da_sun_zenith = self._get_sun_zenith_from_provided_data(projectables, optional_datasets)
 
-        logger.info("Getting emissive part of %s", _nir.attrs["name"])
-        emissivity = self._get_emissivity_as_dask(da_nir, da_tb11, da_tb13_4, da_sun_zenith, _nir.attrs)
+        logger.info("Getting emissive part of %s", nir.attrs["name"])
+        emissivity = self._get_emissivity_as_dask(nir.data, da_tb11, da_tb13_4, da_sun_zenith, nir.attrs)
 
-        proj = self._create_modified_dataarray(emissivity, base_dataarray=_nir)
+        proj = self._create_modified_dataarray(emissivity, base_dataarray=nir)
         proj.attrs["units"] = "K"
         return proj
 

satpy/tests/test_modifiers.py

@@ -273,7 +273,7 @@ def test_no_sunz_no_co2(self, calculator, apply_modifier_info, sza):
 
         # due to copying of DataArrays, self.get_lonlats is not the same as the one that was called
         # we must used the area from the final result DataArray
-        res.attrs["area"].get_lonlats.assert_called()
+        res.attrs["area"].get_lonlats.assert_called_with(chunks=((2,), (2,)), dtype=self.nir.dtype)
         sza.assert_called_with(self.start_time, self.lons, self.lats)
         self.refl_from_tbs.assert_called_with(self.da_sunz, self.nir.data, self.ir_.data, tb_ir_co2=None)
         assert np.allclose(res.data, self.refl * 100).compute()