From a1101f05a2cef8ef5f478da0c39e0ef4ee304f00 Mon Sep 17 00:00:00 2001 From: daanvaningen Date: Thu, 7 Dec 2023 10:09:28 +0100 Subject: [PATCH] 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