From deb8a39b20d58154ff9a4ea7beb3bd45780660b6 Mon Sep 17 00:00:00 2001
From: daanvaningen <daanvaningen@gmail.com>
Date: Wed, 6 Dec 2023 13:44:13 +0100
Subject: [PATCH 1/4] wip

---
 libthreedigrid/_fgrid.pyf                    |  4 +--
 libthreedigrid/cells.f90                     | 26 ++++++++++----------
 libthreedigrid/geo_utils.f90                 |  6 ++---
 libthreedigrid/quadtree.f90                  |  6 ++---
 threedigrid_builder/grid/quadtree.py         |  3 +++
 threedigrid_builder/interface/raster_gdal.py |  9 +++----
 6 files changed, 28 insertions(+), 26 deletions(-)

diff --git a/libthreedigrid/_fgrid.pyf b/libthreedigrid/_fgrid.pyf
index e306e5de..5f98868e 100644
--- a/libthreedigrid/_fgrid.pyf
+++ b/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
+                logical 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
+                logical dimension(:,:),intent(inout) :: area_mask
                 integer dimension(:,:),intent(inout) :: line
                 integer dimension(:,:),intent(inout) :: cross_pix_coords
                 integer intent(in) :: n_line_u
diff --git a/libthreedigrid/cells.f90 b/libthreedigrid/cells.f90
index 3310be25..b5463b38 100644
--- a/libthreedigrid/cells.f90
+++ b/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.
+        logical, 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(:, :)
+        logical, 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(:,:)
+        logical, intent(in) :: area_mask(:,:)
         integer, intent(in) :: quad_idx(:,:)
         integer, intent(in), optional :: nod
         integer, intent(inout), optional :: line(:,:)
@@ -111,7 +111,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
         !!!!! U - direction !!!!!
         !!!!!!!!!!!!!!!!!!!!!!!!!
         if (mn(3) < size(quad_idx, 1)) then
-            if(lg(mn(3)+1, mn(2)) == k .and. any(minval(area_mask(i1:i1+1,j0:j1), 1) > 0)) then
+            if(any(area_mask(i1:i1+1,j0:j1), 1)) then
                 l_u = l_u + 1
                 if (present(line)) then
                     neighbour = quad_idx(mn(3)+1, mn(2)) 
@@ -120,7 +120,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
                     cross_pix_coords(l_u,:) = (/ i1, j0 - 1, i1, j1 /) ! Python indexing so -1 at start slice
                 endif
             else
-                if (lg(mn(3)+1, mn(2)) == k-1 .and. any(minval(area_mask(i1:i1+1,j0:j2), 1) > 0)) then
+                if (lg(mn(3)+1, mn(2)) == k-1 .and. any(area_mask(i1:i1+1,j0:j2), 1)) then
                     l_u = l_u + 1
                     if (present(line)) then
                         neighbour = quad_idx(mn(3)+1, mn(2)) 
@@ -129,7 +129,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
                         cross_pix_coords(l_u,:) = (/ i1, j0 - 1, i1, j2 /) ! Python indexing so -1 at start slice
                     endif
                 endif
-                if (lg(mn(3)+1, mn(4)) == k-1 .and. any(minval(area_mask(i1:i1+1,j3:j1), 1) > 0)) then
+                if (lg(mn(3)+1, mn(4)) == k-1 .and. any(area_mask(i1:i1+1,j3:j1), 1)) then
                     l_u = l_u + 1
                     if (present(line)) then
                         neighbour = quad_idx(mn(3)+1,mn(4))
@@ -142,7 +142,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
         endif
 
         if (mn(1) > 1) then
-            if(lg(mn(1)-1, mn(2)) == k-1 .and. any(minval(area_mask(i0-1:i0,j0:j2), 1) > 0)) then
+            if(lg(mn(1)-1, mn(2)) == k-1 .and. any(area_mask(i0-1:i0,j0:j2), 2)) then
                 l_u = l_u + 1
                 if (present(line)) then
                     neighbour = quad_idx(mn(1)-1, mn(2))
@@ -151,7 +151,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
                     cross_pix_coords(l_u,:) = (/ i0 - 1, j0 - 1, i0 - 1, j2 /)
                 endif
             endif
-            if(lg(mn(1)-1, mn(4)) == k-1 .and. any(minval(area_mask(i0-1:i0,j3:j1), 1) > 0)) then
+            if(lg(mn(1)-1, mn(4)) == k-1 .and. any(area_mask(i0-1:i0,j3:j1), 2)) then
                 l_u = l_u + 1
                 if (present(line)) then
                     neighbour = quad_idx(mn(1)-1, mn(4))
@@ -170,7 +170,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
         !!!!! V - direction !!!!!
         !!!!!!!!!!!!!!!!!!!!!!!!!
         if (mn(4) < size(quad_idx, 2)) then
-            if (lg(mn(1), mn(4)+1) == k .and. any(minval(area_mask(i0:i1,j1:j1+1), 2) > 0)) then
+            if (lg(mn(1), mn(4)+1) == k .and. any(area_mask(i0:i1,j1:j1+1), 2)) then
                 l_v = l_v + 1
                 if (present(line)) then
                     neighbour = quad_idx(mn(1), mn(4)+1)
@@ -179,7 +179,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
                     cross_pix_coords(l_v,:) = (/ i0 - 1, j1, i1, j1 /) ! Python indexing so -1 at start slice
                 endif
             else
-                if (lg(mn(1), mn(4)+1) == k-1 .and. any(minval(area_mask(i0:i2,j1:j1+1), 2) > 0)) then
+                if (lg(mn(1), mn(4)+1) == k-1 .and. any(area_mask(i0:i2,j1:j1+1), 2)) then
                     l_v = l_v + 1
                     if (present(line)) then
                         neighbour = quad_idx(mn(1), mn(4)+1)
@@ -188,7 +188,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
                         cross_pix_coords(l_v,:) = (/ i0 - 1, j1, i2, j1 /) ! Python indexing so -1 at start slice
                     endif
                 endif
-                if (lg(mn(3), mn(4)+1) == k-1 .and. any(minval(area_mask(i3:i1,j1:j1+1), 2) > 0)) then
+                if (lg(mn(3), mn(4)+1) == k-1 .and. any(area_mask(i3:i1,j1:j1+1), 2)) then
                     l_v = l_v + 1
                     if (present(line)) then
                         neighbour = quad_idx(mn(3), mn(4)+1)
@@ -201,7 +201,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
         endif
             
         if (mn(2) > 1) then
-            if(lg(mn(1), mn(2)-1) == k-1 .and. any(minval(area_mask(i0:i2,max(1,j0-1):j0), 2) > 0)) then
+            if(lg(mn(1), mn(2)-1) == k-1 .and. any(area_mask(i0:i2,max(1,j0-1):j0), 2)) then
                 l_v = l_v + 1
                 if (present(line)) then
                     neighbour = quad_idx(mn(1), mn(2)-1)
@@ -210,7 +210,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
                     cross_pix_coords(l_v,:) = (/ i0 - 1, j0 - 1, i2, j0 - 1/) ! Python indexing so -1 at start slice
                 endif
             endif
-            if(lg(mn(3), mn(2)-1) == k-1 .and. any(minval(area_mask(i3:i1,max(1,j0-1):j0), 2) > 0)) then
+            if(lg(mn(3), mn(2)-1) == k-1 .and. any(area_mask(i3:i1,max(1,j0-1):j0), 2)) then
                 l_v = l_v + 1
                 if (present(line)) then
                     neighbour = quad_idx(mn(3), mn(2)-1)
diff --git a/libthreedigrid/geo_utils.f90 b/libthreedigrid/geo_utils.f90
index 61e41cfc..6a6bda00 100644
--- a/libthreedigrid/geo_utils.f90
+++ b/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(:,:)
+        logical, intent(in) :: raster(:,:)
         integer, intent(in) :: i0
         integer, intent(in) :: i1
         integer, intent(in) :: j0
         integer, intent(in) :: j1
-        integer*2, allocatable :: padded_raster(:, :)
+        logical, allocatable :: padded_raster(:, :)
         integer :: i_size, j_size, size_raster_i, size_raster_j
 
         size_raster_i = size(raster, 1)
@@ -92,7 +92,7 @@ function pad_area_mask(raster, i0, i1, j0, j1) result(padded_raster)
         i_size = max(0, i1, size_raster_i)
         j_size = max(0, j1, size_raster_j)
         allocate(padded_raster(1:i_size, 1:j_size))
-        padded_raster = 0
+        padded_raster = .FALSE.
 
         padded_raster(1:size_raster_i, 1:size_raster_j) = raster(1:size_raster_i, 1:size_raster_j)
 
diff --git a/libthreedigrid/quadtree.f90 b/libthreedigrid/quadtree.f90
index 7066f945..480e125e 100644
--- a/libthreedigrid/quadtree.f90
+++ b/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.
+        logical, 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(:,:)
+        logical, 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(:, :)
+        logical, allocatable:: area_mask_padded(:, :)
         integer :: k
         integer :: m,n
         integer :: mn(4)
diff --git a/threedigrid_builder/grid/quadtree.py b/threedigrid_builder/grid/quadtree.py
index 12273bff..e64d221b 100644
--- a/threedigrid_builder/grid/quadtree.py
+++ b/threedigrid_builder/grid/quadtree.py
@@ -107,6 +107,9 @@ def __init__(
         self.n_lines_u = np.array(0, dtype=np.int32, order="F")
         self.n_lines_v = np.array(0, dtype=np.int32, order="F")
 
+        import ipdb
+
+        ipdb.set_trace()
         m_quadtree.make_quadtree(
             self.kmax,
             self.mmax,
diff --git a/threedigrid_builder/interface/raster_gdal.py b/threedigrid_builder/interface/raster_gdal.py
index 691a9063..bc772808 100644
--- a/threedigrid_builder/interface/raster_gdal.py
+++ b/threedigrid_builder/interface/raster_gdal.py
@@ -35,15 +35,14 @@ def __exit__(self, *args, **kwargs):
     def read(self):
         width, height, bbox, mask = self._create_area_arr_from_dem()
 
-        return {
+        data = {
             "pixel_size": self.pixel_size,
             "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=bool, copy=True, order="F"),
         }
+        return data
 
     def _create_area_arr_from_dem(self):
         xpixel, _, xmin, _, ypixel, ymax = self.transform
@@ -59,7 +58,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=bool)
         n_blocks_j = ((width - 1) // size_j) + 1
         n_blocks_i = ((height - 1) // size_i) + 1
 

From be1975394be4cb9b641256c23a03e4a786d36cb8 Mon Sep 17 00:00:00 2001
From: daanvaningen <daanvaningen@gmail.com>
Date: Thu, 7 Dec 2023 09:47:28 +0100
Subject: [PATCH 2/4] bool

---
 CMakeLists.txt                                |  1 +
 libthreedigrid/_fgrid.pyf                     |  4 +--
 libthreedigrid/cells.f90                      | 26 +++++++++----------
 libthreedigrid/geo_utils.f90                  |  4 +--
 libthreedigrid/quadtree.f90                   | 18 +++++++++----
 setup.py                                      |  2 +-
 threedigrid_builder/grid/quadtree.py          |  4 +++
 threedigrid_builder/interface/raster_gdal.py  |  4 ++-
 threedigrid_builder/tests/test_quadtree.py    |  5 ++--
 .../tests/test_raster_interfaces.py           |  2 +-
 10 files changed, 43 insertions(+), 27 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 34e4f6b1..2b300970 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -119,6 +119,7 @@ add_custom_command(
       -m "numpy.f2py"
       -m ${f2py_module_name}
       --lower
+      --verbose
       ${fortran_src_file}
       #only: make_quadtree set_2d_computational_nodes_lines
   WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
diff --git a/libthreedigrid/_fgrid.pyf b/libthreedigrid/_fgrid.pyf
index 5f98868e..6dbec7f3 100644
--- a/libthreedigrid/_fgrid.pyf
+++ b/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
-                logical dimension(:,:),intent(in) :: area_mask
+                logical*1 dimension(:,:),intent(inout) :: 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
-                logical dimension(:,:),intent(inout) :: area_mask
+                logical*1 dimension(:,:),intent(inplace) :: area_mask
                 integer dimension(:,:),intent(inout) :: line
                 integer dimension(:,:),intent(inout) :: cross_pix_coords
                 integer intent(in) :: n_line_u
diff --git a/libthreedigrid/cells.f90 b/libthreedigrid/cells.f90
index b5463b38..47e04c04 100644
--- a/libthreedigrid/cells.f90
+++ b/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
-        logical, intent(inout) :: area_mask(:, :) ! Array with active pixels of model.
+        logical*1, intent(in) :: 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.
-        logical, allocatable :: area_mask_padded(:, :)
+        logical*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
-        logical, intent(in) :: area_mask(:,:)
+        logical*1, intent(in) :: area_mask(:,:)
         integer, intent(in) :: quad_idx(:,:)
         integer, intent(in), optional :: nod
         integer, intent(inout), optional :: line(:,:)
@@ -111,7 +111,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
         !!!!! U - direction !!!!!
         !!!!!!!!!!!!!!!!!!!!!!!!!
         if (mn(3) < size(quad_idx, 1)) then
-            if(any(area_mask(i1:i1+1,j0:j1), 1)) then
+            if(lg(mn(3)+1, mn(2)) == k .and. any(area_mask(i1:i1+1,j0:j1))) then
                 l_u = l_u + 1
                 if (present(line)) then
                     neighbour = quad_idx(mn(3)+1, mn(2)) 
@@ -120,7 +120,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
                     cross_pix_coords(l_u,:) = (/ i1, j0 - 1, i1, j1 /) ! Python indexing so -1 at start slice
                 endif
             else
-                if (lg(mn(3)+1, mn(2)) == k-1 .and. any(area_mask(i1:i1+1,j0:j2), 1)) then
+                if (lg(mn(3)+1, mn(2)) == k-1 .and. any(area_mask(i1:i1+1,j0:j2))) then
                     l_u = l_u + 1
                     if (present(line)) then
                         neighbour = quad_idx(mn(3)+1, mn(2)) 
@@ -129,7 +129,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
                         cross_pix_coords(l_u,:) = (/ i1, j0 - 1, i1, j2 /) ! Python indexing so -1 at start slice
                     endif
                 endif
-                if (lg(mn(3)+1, mn(4)) == k-1 .and. any(area_mask(i1:i1+1,j3:j1), 1)) then
+                if (lg(mn(3)+1, mn(4)) == k-1 .and. any(area_mask(i1:i1+1,j3:j1))) then
                     l_u = l_u + 1
                     if (present(line)) then
                         neighbour = quad_idx(mn(3)+1,mn(4))
@@ -142,7 +142,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
         endif
 
         if (mn(1) > 1) then
-            if(lg(mn(1)-1, mn(2)) == k-1 .and. any(area_mask(i0-1:i0,j0:j2), 2)) then
+            if(lg(mn(1)-1, mn(2)) == k-1 .and. any(area_mask(i0-1:i0,j0:j2))) then
                 l_u = l_u + 1
                 if (present(line)) then
                     neighbour = quad_idx(mn(1)-1, mn(2))
@@ -151,7 +151,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
                     cross_pix_coords(l_u,:) = (/ i0 - 1, j0 - 1, i0 - 1, j2 /)
                 endif
             endif
-            if(lg(mn(1)-1, mn(4)) == k-1 .and. any(area_mask(i0-1:i0,j3:j1), 2)) then
+            if(lg(mn(1)-1, mn(4)) == k-1 .and. any(area_mask(i0-1:i0,j3:j1))) then
                 l_u = l_u + 1
                 if (present(line)) then
                     neighbour = quad_idx(mn(1)-1, mn(4))
@@ -170,7 +170,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
         !!!!! V - direction !!!!!
         !!!!!!!!!!!!!!!!!!!!!!!!!
         if (mn(4) < size(quad_idx, 2)) then
-            if (lg(mn(1), mn(4)+1) == k .and. any(area_mask(i0:i1,j1:j1+1), 2)) then
+            if (lg(mn(1), mn(4)+1) == k .and. any(area_mask(i0:i1,j1:j1+1))) then
                 l_v = l_v + 1
                 if (present(line)) then
                     neighbour = quad_idx(mn(1), mn(4)+1)
@@ -179,7 +179,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
                     cross_pix_coords(l_v,:) = (/ i0 - 1, j1, i1, j1 /) ! Python indexing so -1 at start slice
                 endif
             else
-                if (lg(mn(1), mn(4)+1) == k-1 .and. any(area_mask(i0:i2,j1:j1+1), 2)) then
+                if (lg(mn(1), mn(4)+1) == k-1 .and. any(area_mask(i0:i2,j1:j1+1))) then
                     l_v = l_v + 1
                     if (present(line)) then
                         neighbour = quad_idx(mn(1), mn(4)+1)
@@ -188,7 +188,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
                         cross_pix_coords(l_v,:) = (/ i0 - 1, j1, i2, j1 /) ! Python indexing so -1 at start slice
                     endif
                 endif
-                if (lg(mn(3), mn(4)+1) == k-1 .and. any(area_mask(i3:i1,j1:j1+1), 2)) then
+                if (lg(mn(3), mn(4)+1) == k-1 .and. any(area_mask(i3:i1,j1:j1+1))) then
                     l_v = l_v + 1
                     if (present(line)) then
                         neighbour = quad_idx(mn(3), mn(4)+1)
@@ -201,7 +201,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
         endif
             
         if (mn(2) > 1) then
-            if(lg(mn(1), mn(2)-1) == k-1 .and. any(area_mask(i0:i2,max(1,j0-1):j0), 2)) then
+            if(lg(mn(1), mn(2)-1) == k-1 .and. any(area_mask(i0:i2,max(1,j0-j0)))) then
                 l_v = l_v + 1
                 if (present(line)) then
                     neighbour = quad_idx(mn(1), mn(2)-1)
@@ -210,7 +210,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
                     cross_pix_coords(l_v,:) = (/ i0 - 1, j0 - 1, i2, j0 - 1/) ! Python indexing so -1 at start slice
                 endif
             endif
-            if(lg(mn(3), mn(2)-1) == k-1 .and. any(area_mask(i3:i1,max(1,j0-1):j0), 2)) then
+            if(lg(mn(3), mn(2)-1) == k-1 .and. any(area_mask(i3:i1,max(1,j0-j0)))) then
                 l_v = l_v + 1
                 if (present(line)) then
                     neighbour = quad_idx(mn(3), mn(2)-1)
diff --git a/libthreedigrid/geo_utils.f90 b/libthreedigrid/geo_utils.f90
index 6a6bda00..8c891b73 100644
--- a/libthreedigrid/geo_utils.f90
+++ b/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)
 
-        logical, intent(in) :: raster(:,:)
+        logical*1, intent(in) :: raster(:,:)
         integer, intent(in) :: i0
         integer, intent(in) :: i1
         integer, intent(in) :: j0
         integer, intent(in) :: j1
-        logical, allocatable :: padded_raster(:, :)
+        logical*1, allocatable :: padded_raster(:, :)
         integer :: i_size, j_size, size_raster_i, size_raster_j
 
         size_raster_i = size(raster, 1)
diff --git a/libthreedigrid/quadtree.f90 b/libthreedigrid/quadtree.f90
index 480e125e..7a4050f0 100644
--- a/libthreedigrid/quadtree.f90
+++ b/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
-        logical, intent(in) :: area_mask(:, :) ! Array with active pixels of model.
+        logical*1, intent(inout) :: 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
@@ -25,7 +25,9 @@ subroutine make_quadtree(kmax, mmax, nmax, lgrmin, use_2d_flow, area_mask, lg, q
         integer :: k
         integer :: m, n
 
+        write(*,*) sizeof(area_mask)
         write(*,*) '** INFO: Start making quadtree.'
+        call sleep(60)
         do m=1, mmax(kmax)
             do n=1, nmax(kmax)
                 call divide(kmax, m, n, lg)
@@ -121,11 +123,11 @@ subroutine find_active_2d_comp_cells(&
         integer, intent(in) :: lgrmin
         logical, intent(in) :: use_2d_flow
         integer, intent(inout) :: lg(:,:)
-        logical, intent(in) :: area_mask(:,:)
+        logical*1, intent(in) :: area_mask(:,:)
         integer, intent(inout) :: quad_idx(:,:)
         integer, intent(inout) :: n_line_u
         integer, intent(inout) :: n_line_v
-        logical, allocatable:: area_mask_padded(:, :)
+        logical*1, allocatable:: area_mask_padded(:, :)
         integer :: k
         integer :: m,n
         integer :: mn(4)
@@ -137,8 +139,14 @@ subroutine find_active_2d_comp_cells(&
         n_line_u = 0
         n_line_v = 0
         quad_idx = 0
+
+        write(*,*) "Find active cells"
+        call sleep(60)
         call get_pix_corners(kmax, mmax(kmax), nmax(kmax), lgrmin, i0, i1, j0, j1)
-        area_mask_padded = pad_area_mask(area_mask, i0, i1, j0, j1) 
+        area_mask_padded = pad_area_mask(area_mask, i0, i1, j0, j1)
+        write(*,*) "pad area mask"
+        write(*,*) sizeof(area_mask_padded)
+
         do k=kmax,1,-1
             do m=1,mmax(k)
                 do n=1,nmax(k)
@@ -147,7 +155,7 @@ subroutine find_active_2d_comp_cells(&
                     i1 = min(i1, size(area_mask, 1))
                     j1 = min(j1, size(area_mask, 2))
                     if (all(lg(mn(1):mn(3),mn(2):mn(4)) == k)) then !! TODO: CHECK OF MODEL AREA CHECK IS NECESSARY???
-                        if (all(area_mask_padded(i0:i1, j0:j1) == 0)) then
+                        if (all(area_mask_padded(i0:i1, j0:j1))) then
                             lg(mn(1):mn(3),mn(2):mn(4)) = -99
                         else
                             n_cells = n_cells + 1
diff --git a/setup.py b/setup.py
index d301f23e..0dc8d954 100644
--- a/setup.py
+++ b/setup.py
@@ -34,7 +34,7 @@
 if sys.platform == "win32":
     cmake_args = ["-G", "MSYS Makefiles", "-DCMAKE_GNUtoMS=ON"]
 else:
-    cmake_args = []
+    cmake_args = ["-DF2PY_REPORT_ON_ARRAY_COPY=1", "-DCMAKE_BUILD_TYPE=DEBUG"]
 
 long_description = open("README.rst").read()
 
diff --git a/threedigrid_builder/grid/quadtree.py b/threedigrid_builder/grid/quadtree.py
index e64d221b..d8441a25 100644
--- a/threedigrid_builder/grid/quadtree.py
+++ b/threedigrid_builder/grid/quadtree.py
@@ -107,6 +107,7 @@ def __init__(
         self.n_lines_u = np.array(0, dtype=np.int32, order="F")
         self.n_lines_v = np.array(0, dtype=np.int32, order="F")
 
+        print("make_quadtree")
         import ipdb
 
         ipdb.set_trace()
@@ -197,6 +198,9 @@ def get_nodes_lines(self, area_mask, node_id_counter, line_id_counter):
         line = np.empty((total_lines, 2), dtype=np.int32, order="F")
         cross_pix_coords = np.full((total_lines, 4), -9999, dtype=np.int32, order="F")
 
+        import ipdb
+
+        ipdb.set_trace()
         m_cells.set_2d_computational_nodes_lines(
             np.array([self.origin[0], self.origin[1]], dtype=np.float64),
             self.min_cell_pixels,
diff --git a/threedigrid_builder/interface/raster_gdal.py b/threedigrid_builder/interface/raster_gdal.py
index bc772808..863a3a71 100644
--- a/threedigrid_builder/interface/raster_gdal.py
+++ b/threedigrid_builder/interface/raster_gdal.py
@@ -40,7 +40,9 @@ def read(self):
             "width": width,
             "height": height,
             "bbox": bbox,
-            "area_mask": np.flipud(mask).T.astype(dtype=bool, copy=True, order="F"),
+            "area_mask": np.asfortranarray(
+                np.flipud(mask).T.astype(dtype=np.bool_, copy=True)
+            ),
         }
         return data
 
diff --git a/threedigrid_builder/tests/test_quadtree.py b/threedigrid_builder/tests/test_quadtree.py
index 3036141f..6d61ea00 100644
--- a/threedigrid_builder/tests/test_quadtree.py
+++ b/threedigrid_builder/tests/test_quadtree.py
@@ -14,8 +14,9 @@
 def subgrid_meta():
     width = 20
     height = 16
-    mask = np.ones((width, height), dtype=np.int16, order="F")
-    mask[15:, :] = 0
+    mask = np.array((width, height), dtype=bool, order="F")
+    mask = True
+    mask[15:, :] = False
     return {
         "pixel_size": 0.5,
         "width": width,
diff --git a/threedigrid_builder/tests/test_raster_interfaces.py b/threedigrid_builder/tests/test_raster_interfaces.py
index 0f7b5c00..aeb3023d 100644
--- a/threedigrid_builder/tests/test_raster_interfaces.py
+++ b/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 == bool
         assert np.count_nonzero(result["area_mask"] == 1) == 47180799
         assert result["pixel_size"] == 0.5
         assert result["width"] == 9517

From a1101f05a2cef8ef5f478da0c39e0ef4ee304f00 Mon Sep 17 00:00:00 2001
From: daanvaningen <daanvaningen@gmail.com>
Date: Thu, 7 Dec 2023 10:09:28 +0100
Subject: [PATCH 3/4] int 8

---
 libthreedigrid/_fgrid.pyf                     |  4 +--
 libthreedigrid/cells.f90                      | 26 +++++++++----------
 libthreedigrid/geo_utils.f90                  |  6 ++---
 libthreedigrid/quadtree.f90                   | 18 ++++---------
 threedigrid_builder/grid/quadtree.py          |  7 -----
 threedigrid_builder/interface/raster_gdal.py  |  6 ++---
 threedigrid_builder/tests/test_quadtree.py    |  5 ++--
 .../tests/test_raster_interfaces.py           |  2 +-
 8 files changed, 28 insertions(+), 46 deletions(-)

diff --git a/libthreedigrid/_fgrid.pyf b/libthreedigrid/_fgrid.pyf
index 6dbec7f3..e4622656 100644
--- a/libthreedigrid/_fgrid.pyf
+++ b/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
-                logical*1 dimension(:,:),intent(inout) :: 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
-                logical*1 dimension(:,:),intent(inplace) :: 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
diff --git a/libthreedigrid/cells.f90 b/libthreedigrid/cells.f90
index 47e04c04..f2b377d3 100644
--- a/libthreedigrid/cells.f90
+++ b/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
-        logical*1, intent(in) :: 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.
-        logical*1, 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
-        logical*1, 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(:,:)
@@ -111,7 +111,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
         !!!!! U - direction !!!!!
         !!!!!!!!!!!!!!!!!!!!!!!!!
         if (mn(3) < size(quad_idx, 1)) then
-            if(lg(mn(3)+1, mn(2)) == k .and. any(area_mask(i1:i1+1,j0:j1))) then
+            if(lg(mn(3)+1, mn(2)) == k .and. any(minval(area_mask(i1:i1+1,j0:j1), 1) > 0)) then
                 l_u = l_u + 1
                 if (present(line)) then
                     neighbour = quad_idx(mn(3)+1, mn(2)) 
@@ -120,7 +120,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
                     cross_pix_coords(l_u,:) = (/ i1, j0 - 1, i1, j1 /) ! Python indexing so -1 at start slice
                 endif
             else
-                if (lg(mn(3)+1, mn(2)) == k-1 .and. any(area_mask(i1:i1+1,j0:j2))) then
+                if (lg(mn(3)+1, mn(2)) == k-1 .and. any(minval(area_mask(i1:i1+1,j0:j2), 1) > 0)) then
                     l_u = l_u + 1
                     if (present(line)) then
                         neighbour = quad_idx(mn(3)+1, mn(2)) 
@@ -129,7 +129,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
                         cross_pix_coords(l_u,:) = (/ i1, j0 - 1, i1, j2 /) ! Python indexing so -1 at start slice
                     endif
                 endif
-                if (lg(mn(3)+1, mn(4)) == k-1 .and. any(area_mask(i1:i1+1,j3:j1))) then
+                if (lg(mn(3)+1, mn(4)) == k-1 .and. any(minval(area_mask(i1:i1+1,j3:j1), 1) > 0)) then
                     l_u = l_u + 1
                     if (present(line)) then
                         neighbour = quad_idx(mn(3)+1,mn(4))
@@ -142,7 +142,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
         endif
 
         if (mn(1) > 1) then
-            if(lg(mn(1)-1, mn(2)) == k-1 .and. any(area_mask(i0-1:i0,j0:j2))) then
+            if(lg(mn(1)-1, mn(2)) == k-1 .and. any(minval(area_mask(i0-1:i0,j0:j2), 1) > 0)) then
                 l_u = l_u + 1
                 if (present(line)) then
                     neighbour = quad_idx(mn(1)-1, mn(2))
@@ -151,7 +151,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
                     cross_pix_coords(l_u,:) = (/ i0 - 1, j0 - 1, i0 - 1, j2 /)
                 endif
             endif
-            if(lg(mn(1)-1, mn(4)) == k-1 .and. any(area_mask(i0-1:i0,j3:j1))) then
+            if(lg(mn(1)-1, mn(4)) == k-1 .and. any(minval(area_mask(i0-1:i0,j3:j1), 1) > 0)) then
                 l_u = l_u + 1
                 if (present(line)) then
                     neighbour = quad_idx(mn(1)-1, mn(4))
@@ -170,7 +170,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
         !!!!! V - direction !!!!!
         !!!!!!!!!!!!!!!!!!!!!!!!!
         if (mn(4) < size(quad_idx, 2)) then
-            if (lg(mn(1), mn(4)+1) == k .and. any(area_mask(i0:i1,j1:j1+1))) then
+            if (lg(mn(1), mn(4)+1) == k .and. any(minval(area_mask(i0:i1,j1:j1+1), 2) > 0)) then
                 l_v = l_v + 1
                 if (present(line)) then
                     neighbour = quad_idx(mn(1), mn(4)+1)
@@ -179,7 +179,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
                     cross_pix_coords(l_v,:) = (/ i0 - 1, j1, i1, j1 /) ! Python indexing so -1 at start slice
                 endif
             else
-                if (lg(mn(1), mn(4)+1) == k-1 .and. any(area_mask(i0:i2,j1:j1+1))) then
+                if (lg(mn(1), mn(4)+1) == k-1 .and. any(minval(area_mask(i0:i2,j1:j1+1), 2) > 0)) then
                     l_v = l_v + 1
                     if (present(line)) then
                         neighbour = quad_idx(mn(1), mn(4)+1)
@@ -188,7 +188,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
                         cross_pix_coords(l_v,:) = (/ i0 - 1, j1, i2, j1 /) ! Python indexing so -1 at start slice
                     endif
                 endif
-                if (lg(mn(3), mn(4)+1) == k-1 .and. any(area_mask(i3:i1,j1:j1+1))) then
+                if (lg(mn(3), mn(4)+1) == k-1 .and. any(minval(area_mask(i3:i1,j1:j1+1), 2) > 0)) then
                     l_v = l_v + 1
                     if (present(line)) then
                         neighbour = quad_idx(mn(3), mn(4)+1)
@@ -201,7 +201,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
         endif
             
         if (mn(2) > 1) then
-            if(lg(mn(1), mn(2)-1) == k-1 .and. any(area_mask(i0:i2,max(1,j0-j0)))) then
+            if(lg(mn(1), mn(2)-1) == k-1 .and. any(minval(area_mask(i0:i2,max(1,j0-1):j0), 2) > 0)) then
                 l_v = l_v + 1
                 if (present(line)) then
                     neighbour = quad_idx(mn(1), mn(2)-1)
@@ -210,7 +210,7 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
                     cross_pix_coords(l_v,:) = (/ i0 - 1, j0 - 1, i2, j0 - 1/) ! Python indexing so -1 at start slice
                 endif
             endif
-            if(lg(mn(3), mn(2)-1) == k-1 .and. any(area_mask(i3:i1,max(1,j0-j0)))) then
+            if(lg(mn(3), mn(2)-1) == k-1 .and. any(minval(area_mask(i3:i1,max(1,j0-1):j0), 2) > 0)) then
                 l_v = l_v + 1
                 if (present(line)) then
                     neighbour = quad_idx(mn(3), mn(2)-1)
diff --git a/libthreedigrid/geo_utils.f90 b/libthreedigrid/geo_utils.f90
index 8c891b73..556393f8 100644
--- a/libthreedigrid/geo_utils.f90
+++ b/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)
 
-        logical*1, intent(in) :: raster(:,:)
+        integer*1, intent(in) :: raster(:,:)
         integer, intent(in) :: i0
         integer, intent(in) :: i1
         integer, intent(in) :: j0
         integer, intent(in) :: j1
-        logical*1, 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)
@@ -92,7 +92,7 @@ function pad_area_mask(raster, i0, i1, j0, j1) result(padded_raster)
         i_size = max(0, i1, size_raster_i)
         j_size = max(0, j1, size_raster_j)
         allocate(padded_raster(1:i_size, 1:j_size))
-        padded_raster = .FALSE.
+        padded_raster = 0
 
         padded_raster(1:size_raster_i, 1:size_raster_j) = raster(1:size_raster_i, 1:size_raster_j)
 
diff --git a/libthreedigrid/quadtree.f90 b/libthreedigrid/quadtree.f90
index 7a4050f0..53628a8a 100644
--- a/libthreedigrid/quadtree.f90
+++ b/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
-        logical*1, intent(inout) :: 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
@@ -25,9 +25,7 @@ subroutine make_quadtree(kmax, mmax, nmax, lgrmin, use_2d_flow, area_mask, lg, q
         integer :: k
         integer :: m, n
 
-        write(*,*) sizeof(area_mask)
         write(*,*) '** INFO: Start making quadtree.'
-        call sleep(60)
         do m=1, mmax(kmax)
             do n=1, nmax(kmax)
                 call divide(kmax, m, n, lg)
@@ -123,11 +121,11 @@ subroutine find_active_2d_comp_cells(&
         integer, intent(in) :: lgrmin
         logical, intent(in) :: use_2d_flow
         integer, intent(inout) :: lg(:,:)
-        logical*1, 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
-        logical*1, allocatable:: area_mask_padded(:, :)
+        integer*1, allocatable:: area_mask_padded(:, :)
         integer :: k
         integer :: m,n
         integer :: mn(4)
@@ -139,14 +137,8 @@ subroutine find_active_2d_comp_cells(&
         n_line_u = 0
         n_line_v = 0
         quad_idx = 0
-
-        write(*,*) "Find active cells"
-        call sleep(60)
         call get_pix_corners(kmax, mmax(kmax), nmax(kmax), lgrmin, i0, i1, j0, j1)
-        area_mask_padded = pad_area_mask(area_mask, i0, i1, j0, j1)
-        write(*,*) "pad area mask"
-        write(*,*) sizeof(area_mask_padded)
-
+        area_mask_padded = pad_area_mask(area_mask, i0, i1, j0, j1) 
         do k=kmax,1,-1
             do m=1,mmax(k)
                 do n=1,nmax(k)
@@ -155,7 +147,7 @@ subroutine find_active_2d_comp_cells(&
                     i1 = min(i1, size(area_mask, 1))
                     j1 = min(j1, size(area_mask, 2))
                     if (all(lg(mn(1):mn(3),mn(2):mn(4)) == k)) then !! TODO: CHECK OF MODEL AREA CHECK IS NECESSARY???
-                        if (all(area_mask_padded(i0:i1, j0:j1))) then
+                        if (all(area_mask_padded(i0:i1, j0:j1) == 0)) then
                             lg(mn(1):mn(3),mn(2):mn(4)) = -99
                         else
                             n_cells = n_cells + 1
diff --git a/threedigrid_builder/grid/quadtree.py b/threedigrid_builder/grid/quadtree.py
index d8441a25..12273bff 100644
--- a/threedigrid_builder/grid/quadtree.py
+++ b/threedigrid_builder/grid/quadtree.py
@@ -107,10 +107,6 @@ def __init__(
         self.n_lines_u = np.array(0, dtype=np.int32, order="F")
         self.n_lines_v = np.array(0, dtype=np.int32, order="F")
 
-        print("make_quadtree")
-        import ipdb
-
-        ipdb.set_trace()
         m_quadtree.make_quadtree(
             self.kmax,
             self.mmax,
@@ -198,9 +194,6 @@ def get_nodes_lines(self, area_mask, node_id_counter, line_id_counter):
         line = np.empty((total_lines, 2), dtype=np.int32, order="F")
         cross_pix_coords = np.full((total_lines, 4), -9999, dtype=np.int32, order="F")
 
-        import ipdb
-
-        ipdb.set_trace()
         m_cells.set_2d_computational_nodes_lines(
             np.array([self.origin[0], self.origin[1]], dtype=np.float64),
             self.min_cell_pixels,
diff --git a/threedigrid_builder/interface/raster_gdal.py b/threedigrid_builder/interface/raster_gdal.py
index 863a3a71..1878e1ca 100644
--- a/threedigrid_builder/interface/raster_gdal.py
+++ b/threedigrid_builder/interface/raster_gdal.py
@@ -40,9 +40,7 @@ def read(self):
             "width": width,
             "height": height,
             "bbox": bbox,
-            "area_mask": np.asfortranarray(
-                np.flipud(mask).T.astype(dtype=np.bool_, copy=True)
-            ),
+            "area_mask": np.flipud(mask).T.astype(dtype=np.int8, copy=True, order="F"),
         }
         return data
 
@@ -60,7 +58,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=bool)
+        mask = np.zeros((height, width), dtype=np.int8)
         n_blocks_j = ((width - 1) // size_j) + 1
         n_blocks_i = ((height - 1) // size_i) + 1
 
diff --git a/threedigrid_builder/tests/test_quadtree.py b/threedigrid_builder/tests/test_quadtree.py
index 6d61ea00..e7fd6cba 100644
--- a/threedigrid_builder/tests/test_quadtree.py
+++ b/threedigrid_builder/tests/test_quadtree.py
@@ -14,9 +14,8 @@
 def subgrid_meta():
     width = 20
     height = 16
-    mask = np.array((width, height), dtype=bool, order="F")
-    mask = True
-    mask[15:, :] = False
+    mask = np.ones((width, height), dtype=np.int8, order="F")
+    mask[15:, :] = 0
     return {
         "pixel_size": 0.5,
         "width": width,
diff --git a/threedigrid_builder/tests/test_raster_interfaces.py b/threedigrid_builder/tests/test_raster_interfaces.py
index aeb3023d..2ebb761b 100644
--- a/threedigrid_builder/tests/test_raster_interfaces.py
+++ b/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 == bool
+        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

From 3112f214847481bc6aad5eff0e504ee282cfda07 Mon Sep 17 00:00:00 2001
From: daanvaningen <daanvaningen@gmail.com>
Date: Thu, 7 Dec 2023 10:23:01 +0100
Subject: [PATCH 4/4] cleanup

---
 CHANGES.rst                                  | 2 ++
 CMakeLists.txt                               | 1 -
 setup.py                                     | 2 +-
 threedigrid_builder/interface/raster_gdal.py | 3 +--
 4 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index 2a534354..b209139a 100644
--- a/CHANGES.rst
+++ b/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)
 -------------------
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2b300970..34e4f6b1 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -119,7 +119,6 @@ add_custom_command(
       -m "numpy.f2py"
       -m ${f2py_module_name}
       --lower
-      --verbose
       ${fortran_src_file}
       #only: make_quadtree set_2d_computational_nodes_lines
   WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
diff --git a/setup.py b/setup.py
index 0dc8d954..d301f23e 100644
--- a/setup.py
+++ b/setup.py
@@ -34,7 +34,7 @@
 if sys.platform == "win32":
     cmake_args = ["-G", "MSYS Makefiles", "-DCMAKE_GNUtoMS=ON"]
 else:
-    cmake_args = ["-DF2PY_REPORT_ON_ARRAY_COPY=1", "-DCMAKE_BUILD_TYPE=DEBUG"]
+    cmake_args = []
 
 long_description = open("README.rst").read()
 
diff --git a/threedigrid_builder/interface/raster_gdal.py b/threedigrid_builder/interface/raster_gdal.py
index 1878e1ca..45388bd1 100644
--- a/threedigrid_builder/interface/raster_gdal.py
+++ b/threedigrid_builder/interface/raster_gdal.py
@@ -35,14 +35,13 @@ def __exit__(self, *args, **kwargs):
     def read(self):
         width, height, bbox, mask = self._create_area_arr_from_dem()
 
-        data = {
+        return {
             "pixel_size": self.pixel_size,
             "width": width,
             "height": height,
             "bbox": bbox,
             "area_mask": np.flipud(mask).T.astype(dtype=np.int8, copy=True, order="F"),
         }
-        return data
 
     def _create_area_arr_from_dem(self):
         xpixel, _, xmin, _, ypixel, ymax = self.transform