Merge pull request #344 from nens/daan-dem-mask-boolean

DEM mask to lower datatype size to reduce memory usage

CHANGES.rst

@@ -6,6 +6,8 @@ Changelog of threedigrid-builder
 
 - Add manhole_indicator field to gridadmin for future export.
 
+- Reduce memory use by reducing DEM mask datatype to int8.
+
 
 1.12.1 (2023-08-14)
 -------------------

libthreedigrid/_fgrid.pyf

@@ -11,7 +11,7 @@ python module _fgrid ! in
                 integer dimension(:),intent(in) :: nmax
                 integer intent(in) :: lgrmin
                 integer intent(in) :: use_2d_flow
-                integer*2 dimension(:,:),intent(in) :: area_mask
+                integer*1 dimension(:,:),intent(in) :: area_mask
                 integer dimension(:,:),intent(inout) :: lg
                 integer dimension(:,:),intent(inout) :: quad_idx
                 integer intent(inout) :: n_cells
@@ -36,7 +36,7 @@ python module _fgrid ! in
                 double precision dimension(:,:),intent(inout) :: bounds
                 double precision dimension(:,:),intent(inout) :: coords
                 integer dimension(:,:),intent(inout) :: pixel_coords
-                integer*2 dimension(:,:),intent(inout) :: area_mask
+                integer*1 dimension(:,:),intent(inout) :: area_mask
                 integer dimension(:,:),intent(inout) :: line
                 integer dimension(:,:),intent(inout) :: cross_pix_coords
                 integer intent(in) :: n_line_u

libthreedigrid/cells.f90

@@ -24,12 +24,12 @@ subroutine set_2d_computational_nodes_lines(origin, lgrmin, kmax, mmax, nmax, dx
         double precision, intent(inout) :: bounds(:, :) ! Bbox of comp cell
         double precision, intent(inout) :: coords(:, :) ! Cell center coordinates
         integer, intent(inout) :: pixel_coords(:, :) ! pixel bbox of comp cell
-        integer*2, intent(inout) :: area_mask(:, :) ! Array with active pixels of model.
+        integer*1, intent(inout) :: area_mask(:, :) ! Array with active pixels of model.
         integer, intent(inout) :: line(:, :) ! Array with connecting nodes of line.
         integer, intent(inout) :: cross_pix_coords(:, :) ! Array pixel indices of line interface
         integer, intent(in) :: n_line_u ! Number of active u-dir lines.
         integer, intent(in) :: n_line_v  ! Number of active v-dir lines.
-        integer*2, allocatable :: area_mask_padded(:, :)
+        integer*1, allocatable :: area_mask_padded(:, :)
         integer :: nod
         integer :: k
         integer :: i0, i1, j0, j1
@@ -95,7 +95,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
         integer, intent(in) :: mn(4)
         integer, intent(in) :: lg(:,:)
         integer, intent(in) :: lgrmin
-        integer*2, intent(in) :: area_mask(:,:)
+        integer*1, intent(in) :: area_mask(:,:)
         integer, intent(in) :: quad_idx(:,:)
         integer, intent(in), optional :: nod
         integer, intent(inout), optional :: line(:,:)

libthreedigrid/geo_utils.f90

@@ -79,12 +79,12 @@ end subroutine crop_pix_coords_to_raster
 
     function pad_area_mask(raster, i0, i1, j0, j1) result(padded_raster)
 
-        integer*2, intent(in) :: raster(:,:)
+        integer*1, intent(in) :: raster(:,:)
         integer, intent(in) :: i0
         integer, intent(in) :: i1
         integer, intent(in) :: j0
         integer, intent(in) :: j1
-        integer*2, allocatable :: padded_raster(:, :)
+        integer*1, allocatable :: padded_raster(:, :)
         integer :: i_size, j_size, size_raster_i, size_raster_j
 
         size_raster_i = size(raster, 1)

libthreedigrid/quadtree.f90

@@ -16,7 +16,7 @@ subroutine make_quadtree(kmax, mmax, nmax, lgrmin, use_2d_flow, area_mask, lg, q
         integer, intent(in) :: nmax(:) ! Y Dimension of each refinement level
         integer, intent(in) :: lgrmin ! Number of pixels in cell of smallest refinement level
         integer, intent(in) :: use_2d_flow  ! Whether to add flowlines
-        integer*2, intent(in) :: area_mask(:, :) ! Array with active pixels of model.
+        integer*1, intent(in) :: area_mask(:, :) ! Array with active pixels of model.
         integer, intent(inout) :: lg(:, :) ! Array with all refinement levels.
         integer, intent(inout) :: quad_idx(:, :) ! Array with idx of cell at lg refinement locations
         integer, intent(inout) :: n_cells ! counter for active cells
@@ -121,11 +121,11 @@ subroutine find_active_2d_comp_cells(&
         integer, intent(in) :: lgrmin
         logical, intent(in) :: use_2d_flow
         integer, intent(inout) :: lg(:,:)
-        integer*2, intent(in) :: area_mask(:,:)
+        integer*1, intent(in) :: area_mask(:,:)
         integer, intent(inout) :: quad_idx(:,:)
         integer, intent(inout) :: n_line_u
         integer, intent(inout) :: n_line_v
-        integer*2, allocatable:: area_mask_padded(:, :)
+        integer*1, allocatable:: area_mask_padded(:, :)
         integer :: k
         integer :: m,n
         integer :: mn(4)

threedigrid_builder/interface/raster_gdal.py

@@ -40,9 +40,7 @@ def read(self):
             "width": width,
             "height": height,
             "bbox": bbox,
-            "area_mask": np.flipud(mask).T.astype(
-                dtype=np.int16, copy=False, order="F"
-            ),
+            "area_mask": np.flipud(mask).T.astype(dtype=np.int8, copy=True, order="F"),
         }
 
     def _create_area_arr_from_dem(self):
@@ -59,7 +57,7 @@ def _create_area_arr_from_dem(self):
         band = self._dataset.GetRasterBand(1)
         size_j, size_i = band.GetBlockSize()
         nodata = band.GetNoDataValue()
-        mask = np.zeros((height, width), dtype=np.int16)
+        mask = np.zeros((height, width), dtype=np.int8)
         n_blocks_j = ((width - 1) // size_j) + 1
         n_blocks_i = ((height - 1) // size_i) + 1
 

threedigrid_builder/tests/test_quadtree.py

@@ -14,7 +14,7 @@
 def subgrid_meta():
     width = 20
     height = 16
-    mask = np.ones((width, height), dtype=np.int16, order="F")
+    mask = np.ones((width, height), dtype=np.int8, order="F")
     mask[15:, :] = 0
     return {
         "pixel_size": 0.5,

threedigrid_builder/tests/test_raster_interfaces.py

@@ -29,7 +29,7 @@ def test_read(dem_path):
     with GDALInterface(dem_path) as dem:
         result = dem.read()
         assert result["area_mask"].shape == (9517, 9726)
-        assert result["area_mask"].dtype == np.int16
+        assert result["area_mask"].dtype == np.int8
         assert np.count_nonzero(result["area_mask"] == 1) == 47180799
         assert result["pixel_size"] == 0.5
         assert result["width"] == 9517