Skip to content

Commit

Permalink
int 8
Browse files Browse the repository at this point in the history
  • Loading branch information
daanvaningen committed Dec 7, 2023
1 parent be19753 commit a1101f0
Show file tree
Hide file tree
Showing 8 changed files with 28 additions and 46 deletions.
4 changes: 2 additions & 2 deletions libthreedigrid/_fgrid.pyf
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
26 changes: 13 additions & 13 deletions libthreedigrid/cells.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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(:,:)
Expand All @@ -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))
Expand All @@ -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))
Expand All @@ -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))
Expand All @@ -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))
Expand All @@ -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))
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand Down
6 changes: 3 additions & 3 deletions libthreedigrid/geo_utils.f90
Original file line number Diff line number Diff line change
Expand Up @@ -79,20 +79,20 @@ 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)
size_raster_j = size(raster, 2)
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)

Expand Down
18 changes: 5 additions & 13 deletions libthreedigrid/quadtree.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand All @@ -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
Expand Down
7 changes: 0 additions & 7 deletions threedigrid_builder/grid/quadtree.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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,
Expand Down
6 changes: 2 additions & 4 deletions threedigrid_builder/interface/raster_gdal.py
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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

Expand Down
5 changes: 2 additions & 3 deletions threedigrid_builder/tests/test_quadtree.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
2 changes: 1 addition & 1 deletion threedigrid_builder/tests/test_raster_interfaces.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

0 comments on commit a1101f0

Please sign in to comment.