- ENH: Allowing stack to take single band in input instead of a list

- FIX: Fixing a regression loading optical bands which have been previously cleaned (Landsat, Theia, possibly PlanetScope)
- FIX: `load` and `stack` always returns `float32` arrays

CHANGES.md

@@ -2,9 +2,11 @@
 
 ## X.Y.Z (YYYY-MM-DD)
 
-## 0.4.7.post1 (2021-07-23)
-- Allowing `stack` to take single band in input instead of a list
-- Fixing a regression loading Landsat bands which have been previously cleaned (again)
+## 0.4.8 (2021-07-23)
+- ENH: Allowing `stack` to take single band in input instead of a list
+- FIX: Fixing a regression loading optical bands which have been previously cleaned (Landsat, Theia, possibly PlanetScope)
+- FIX: `load` and `stack` always returns `float32` arrays
+- CI: Testing if loading 2 times a band gives the same result
 
 ## 0.4.7.post0 (2021-07-23)
 

CI/SCRIPTS/test_satellites.py

@@ -3,6 +3,7 @@
 import os
 import tempfile
 
+import numpy as np
 import xarray as xr
 from cloudpathlib import AnyPath
 from geopandas import gpd
@@ -213,13 +214,22 @@ def _test_core(pattern: str, prod_dir: str, possible_bands: list, debug=False):
                     # Remove DEM tifs if existing
                     remove_dem(prod)
 
-                    # Get stack bands
+                    # BAND TESTS
                     LOGGER.info("Checking load and stack")
                     # DO NOT RECOMPUTE BANDS WITH SNAP --> WAY TOO SLOW
                     stack_bands = [
                         band for band in possible_bands if prod.has_band(band)
                     ]
+                    first_band = stack_bands[0]
+
+                    # Check that band loaded 2 times gives the same results (disregarding float uncertainties)
+                    band_arr1 = prod.load(first_band, resolution=res)[first_band]
+                    band_arr2 = prod.load(first_band, resolution=res)[first_band]
+                    np.testing.assert_array_almost_equal(band_arr1, band_arr2)
+                    assert band_arr1.dtype == np.float32
+                    assert band_arr2.dtype == np.float32
 
+                    # Get stack bands
                     # Stack data
                     ci_stack = get_ci_data_dir().joinpath(
                         prod.condensed_name, f"{prod.condensed_name}_stack.tif"
@@ -231,6 +241,7 @@ def _test_core(pattern: str, prod_dir: str, possible_bands: list, debug=False):
                     stack = prod.stack(
                         stack_bands, resolution=res, stack_path=curr_path
                     )
+                    assert stack.dtype == np.float32
 
                     # Write to path if needed
                     if not ci_stack.exists():
@@ -241,20 +252,19 @@ def _test_core(pattern: str, prod_dir: str, possible_bands: list, debug=False):
 
                     # Load a band with the size option
                     LOGGER.info("Checking load with size keyword")
-                    band = stack_bands[0]
                     ci_band = get_ci_data_dir().joinpath(
                         prod.condensed_name,
-                        f"{prod.condensed_name}_{band.name}_test.tif",
+                        f"{prod.condensed_name}_{first_band.name}_test.tif",
                     )
                     curr_path_band = os.path.join(
-                        tmp_dir, f"{prod.condensed_name}_{band.name}_test.tif"
+                        tmp_dir, f"{prod.condensed_name}_{first_band.name}_test.tif"
                     )
                     if not ci_band.exists():
                         ci_band = curr_path_band
 
                     band_arr = prod.load(
-                        band, size=(stack.rio.width, stack.rio.height)
-                    )[band]
+                        first_band, size=(stack.rio.width, stack.rio.height)
+                    )[first_band]
                     rasters.write(band_arr, curr_path_band)
                     assert_raster_almost_equal(curr_path_band, ci_band, decimal=4)
 

eoreader/__init__.py

@@ -17,7 +17,7 @@
 """
 **EOReader** library
 """
-__version__ = "0.4.7-1"
+__version__ = "0.4.8"
 __title__ = "eoreader"
 __description__ = (
     "Remote-sensing opensource python library reading optical and SAR sensors, "

eoreader/products/optical/pla_product.py

@@ -414,7 +414,13 @@ def _read_band(
         )
 
         # Compute the correct radiometry of the band
-        band_xda = band_xda / 10000.0
+        original_dtype = band_xda.encoding.get("dtype", band_xda.dtype)
+        if original_dtype == "uint16":
+            band_xda /= 10000.0
+
+        # Convert type if needed
+        if band_xda.dtype != np.float32:
+            band_xda = band_xda.astype(np.float32)
 
         return band_xda
 

eoreader/products/optical/s2_theia_product.py

@@ -239,10 +239,16 @@ def _read_band(
         # Read band
         band_xda = rasters.read(
             path, resolution=resolution, size=size, resampling=Resampling.bilinear
-        ).astype(np.float32)
+        )
+
+        # Compute the correct radiometry of the band (Theia product are stored into int16 bits)
+        original_dtype = band_xda.encoding.get("dtype", band_xda.dtype)
+        if original_dtype == "int16":
+            band_xda /= 10000.0
 
-        # Compute the correct radiometry of the band
-        band_xda = band_xda / 10000.0
+        # Convert type if needed
+        if band_xda.dtype != np.float32:
+            band_xda = band_xda.astype(np.float32)
 
         return band_xda
 

eoreader/products/optical/s3_product.py

@@ -433,7 +433,7 @@ def _read_band(
         # Read band
         return rasters.read(
             path, resolution=resolution, size=size, resampling=Resampling.bilinear
-        )
+        ).astype(np.float32)
 
     # pylint: disable=R0913
     # R0913: Too many arguments (6/5) (too-many-arguments)