Blockwise operation (#702)

* Use blockwise operation for round_mantissa to not overload memory

* Set chunksize for round_mantissa to 1024

* remove upperbound for gdal

* add docstring for chunk_rows

* cast los_east, los_north, and incidence_angle to float32

---------

Co-authored-by: mirzaees <smirzaee@aurora.jpl.nasa.gov>

Author: mirzaees

conda-env.yml

@@ -5,7 +5,7 @@ dependencies:
   - python>=3.9
   - pip>=21.3  # https://pip.pypa.io/en/stable/reference/build-system/pyproject-toml/#editable-installation
   - git  # for pip install, due to setuptools_scm
-  - gdal>=3.3,<3.11
+  - gdal>=3.3
   - libgdal-netcdf
   - libgdal-hdf5
   - h5py>=3.6

src/dolphin/atmosphere/ionosphere.py

@@ -88,12 +88,12 @@ def estimate_ionospheric_delay(
     # Read the incidence angle
     if "los_east" in geom_files:
         # ISCE3 geocoded products
-        los_east = io.load_gdal(geom_files["los_east"])
-        los_north = io.load_gdal(geom_files["los_north"])
+        los_east = io.load_gdal(geom_files["los_east"]).astype(np.float32)
+        los_north = io.load_gdal(geom_files["los_north"]).astype(np.float32)
         inc_angle = np.arccos(np.sqrt(1 - los_east**2 - los_north**2)) * 180 / np.pi
     else:
         # ISCE2 radar coordinate
-        inc_angle = io.load_gdal(geom_files["incidence_angle"])
+        inc_angle = io.load_gdal(geom_files["incidence_angle"]).astype(np.float32)
 
     iono_inc_angle = incidence_angle_ground_to_iono(inc_angle)
 

src/dolphin/io/_utils.py

@@ -233,7 +233,7 @@ def repack_rasters(
     )
 
 
-def round_mantissa(z: np.ndarray, keep_bits: int = 10) -> None:
+def round_mantissa(z: np.ndarray, keep_bits: int = 10, chunk_rows: int = 1024) -> None:
     """Zero out mantissa bits of elements of array in place.
 
     Drops a specified number of bits from the floating point mantissa,
@@ -248,6 +248,9 @@ def round_mantissa(z: np.ndarray, keep_bits: int = 10) -> None:
         Lower numbers will truncate the mantissa more and enable
         more compression.
         Default is 10.
+    chunk_rows : int
+        Number of rows to process at a time to limit memory usage.
+        Default is 1024.
 
     References
     ----------
@@ -261,8 +264,8 @@ def round_mantissa(z: np.ndarray, keep_bits: int = 10) -> None:
     }
     # recurse for complex data
     if np.iscomplexobj(z):
-        round_mantissa(z.real, keep_bits)
-        round_mantissa(z.imag, keep_bits)
+        round_mantissa(z.real, keep_bits, chunk_rows)
+        round_mantissa(z.imag, keep_bits, chunk_rows)
         return
 
     if not z.dtype.kind == "f" or z.dtype.itemsize > 8:
@@ -276,9 +279,10 @@ def round_mantissa(z: np.ndarray, keep_bits: int = 10) -> None:
         return z
     if keep_bits > bits:
         raise ValueError("keep_bits too large for given dtype")
-    b = z.view(a_int_dtype)
     maskbits = bits - keep_bits
     mask = (all_set >> maskbits) << maskbits
     half_quantum1 = (1 << (maskbits - 1)) - 1
-    b += ((b >> maskbits) & 1) + half_quantum1
-    b &= mask
+    for i in range(0, z.shape[0], chunk_rows):
+        b = z[i : i + chunk_rows].view(a_int_dtype)
+        b += ((b >> maskbits) & 1) + half_quantum1
+        b &= mask