diff --git a/src/2d/shallow/topo_module.f90 b/src/2d/shallow/topo_module.f90 index 5734e53ef..7a85b0be6 100644 --- a/src/2d/shallow/topo_module.f90 +++ b/src/2d/shallow/topo_module.f90 @@ -54,29 +54,29 @@ module topo_module ! Per-file preprocessing parameters — read from the 8 new lines in topo.data ! for all topo file types (2, 3, 4, 5). Type-1 files are skipped. - real(kind=8), allocatable :: tp_crop_extent(:,:) ! (4,n): x1,x2,y1,y2 domain coords; all-zero=no crop - integer, allocatable :: tp_coarsen(:) ! subsampling factor (1=none) - integer, allocatable :: tp_buffer(:) ! buffer cells (grid points) around crop region - real(kind=8), allocatable :: tp_align(:,:) ! (2,n): align offset; all-zero=none - real(kind=8), allocatable :: tp_x_shift(:) ! x coordinate registration offset - real(kind=8), allocatable :: tp_y_shift(:) ! y coordinate registration offset - real(kind=8), allocatable :: tp_z_shift(:) ! bulk datum offset - logical, allocatable :: tp_negate_z(:) ! explicit Z sign flip + real(kind=8), allocatable :: topo_crop_extent(:,:) ! (4,n): x1,x2,y1,y2 domain coords; all-zero=no crop + integer, allocatable :: topo_coarsen(:) ! subsampling factor (1=none) + integer, allocatable :: topo_buffer(:) ! buffer cells (grid points) around crop region + real(kind=8), allocatable :: topo_align(:,:) ! (2,n): align offset; all-zero=none + real(kind=8), allocatable :: topo_x_shift(:) ! x coordinate registration offset + real(kind=8), allocatable :: topo_y_shift(:) ! y coordinate registration offset + real(kind=8), allocatable :: topo_z_shift(:) ! bulk datum offset + logical, allocatable :: topo_negate_z(:) ! explicit Z sign flip ! Full (un-cropped) file sizes and kept-window start indices for types 2/3, ! saved by read_topo_settings so read_topo_file can extract the kept ! (cropped/coarsened) lattice: 0-based indices from the file's south-west - ! corner; kept points are tp_i_start + k*tp_coarsen, k = 0..mxtopo-1. + ! corner; kept points are topo_i_start + k*topo_coarsen, k = 0..mxtopo-1. integer, allocatable :: mxtopo_full(:), mytopo_full(:) - integer, allocatable :: tp_i_start(:), tp_j_start(:) + integer, allocatable :: topo_i_start(:), topo_j_start(:) ! Per-file dtopo preprocessing parameters — read from the 8 preprocessing ! lines in dtopo.data. x_shift, y_shift, z_shift and negate_z are ! implemented for dtopography; the others trigger a stop in ! read_dtopo_settings. - real(kind=8), allocatable :: dtp_x_shift(:), dtp_y_shift(:) - real(kind=8), allocatable :: dtp_z_shift(:) - logical, allocatable :: dtp_negate_z(:) + real(kind=8), allocatable :: dtopo_x_shift(:), dtopo_y_shift(:) + real(kind=8), allocatable :: dtopo_z_shift(:) + logical, allocatable :: dtopo_negate_z(:) ! NetCDF descriptor values for type-4 dtopo files, parsed from the ! key=value block in dtopo.data by read_dtopo_netcdf_descriptor. The @@ -226,26 +226,26 @@ subroutine read_topo_settings(restart,file_name) nc_has_fill = .false. nc_has_crop = .false. - allocate(tp_crop_extent(4, mtopofiles), tp_coarsen(mtopofiles)) - allocate(tp_buffer(mtopofiles), tp_align(2, mtopofiles)) - allocate(tp_x_shift(mtopofiles), tp_y_shift(mtopofiles)) - allocate(tp_z_shift(mtopofiles)) - allocate(tp_negate_z(mtopofiles)) - tp_crop_extent = 0.0d0 - tp_coarsen = 1 - tp_buffer = 0 - tp_align = 0.0d0 - tp_x_shift = 0.0d0 - tp_y_shift = 0.0d0 - tp_z_shift = 0.0d0 - tp_negate_z = .false. + allocate(topo_crop_extent(4, mtopofiles), topo_coarsen(mtopofiles)) + allocate(topo_buffer(mtopofiles), topo_align(2, mtopofiles)) + allocate(topo_x_shift(mtopofiles), topo_y_shift(mtopofiles)) + allocate(topo_z_shift(mtopofiles)) + allocate(topo_negate_z(mtopofiles)) + topo_crop_extent = 0.0d0 + topo_coarsen = 1 + topo_buffer = 0 + topo_align = 0.0d0 + topo_x_shift = 0.0d0 + topo_y_shift = 0.0d0 + topo_z_shift = 0.0d0 + topo_negate_z = .false. allocate(mxtopo_full(mtopofiles), mytopo_full(mtopofiles)) - allocate(tp_i_start(mtopofiles), tp_j_start(mtopofiles)) + allocate(topo_i_start(mtopofiles), topo_j_start(mtopofiles)) mxtopo_full = 0 mytopo_full = 0 - tp_i_start = 0 - tp_j_start = 0 + topo_i_start = 0 + topo_j_start = 0 mtopofiles = mtopofiles - num_dtopo ! decrement after allocates @@ -262,26 +262,26 @@ subroutine read_topo_settings(restart,file_name) ! Read the 8 new preprocessing parameter lines present in ! topo.data for ALL file types (written by TopographyData.write()). - read(iunit,*) tp_crop_extent(1,i), tp_crop_extent(2,i), & - tp_crop_extent(3,i), tp_crop_extent(4,i) - read(iunit,*) tp_coarsen(i) - read(iunit,*) tp_buffer(i) - read(iunit,*) tp_align(1,i), tp_align(2,i) - read(iunit,*) tp_x_shift(i) - read(iunit,*) tp_y_shift(i) - read(iunit,*) tp_z_shift(i) - read(iunit,*) tp_negate_z(i) + read(iunit,*) topo_crop_extent(1,i), topo_crop_extent(2,i), & + topo_crop_extent(3,i), topo_crop_extent(4,i) + read(iunit,*) topo_coarsen(i) + read(iunit,*) topo_buffer(i) + read(iunit,*) topo_align(1,i), topo_align(2,i) + read(iunit,*) topo_x_shift(i) + read(iunit,*) topo_y_shift(i) + read(iunit,*) topo_z_shift(i) + read(iunit,*) topo_negate_z(i) ! Preprocessing is not supported for topo_type=1. if (abs(itopotype(i)) == 1 .and. ( & - any(tp_crop_extent(:,i) /= 0.0d0) .or. & - tp_coarsen(i) > 1 .or. & - tp_buffer(i) /= 0 .or. & - any(tp_align(:,i) /= 0.0d0) .or. & - tp_negate_z(i) .or. & - tp_x_shift(i) /= 0.0d0 .or. & - tp_y_shift(i) /= 0.0d0 .or. & - tp_z_shift(i) /= 0.0d0)) then + any(topo_crop_extent(:,i) /= 0.0d0) .or. & + topo_coarsen(i) > 1 .or. & + topo_buffer(i) /= 0 .or. & + any(topo_align(:,i) /= 0.0d0) .or. & + topo_negate_z(i) .or. & + topo_x_shift(i) /= 0.0d0 .or. & + topo_y_shift(i) /= 0.0d0 .or. & + topo_z_shift(i) /= 0.0d0)) then print *, "ERROR: preprocessing attributes (crop, coarsen, shift) are not" print *, " supported for topo_type=1. Convert to type 2/3/4 first." stop 1 @@ -310,22 +310,22 @@ subroutine read_topo_settings(restart,file_name) mxtopo_full(i) = mxtopo(i) mytopo_full(i) = mytopo(i) - if (any(tp_crop_extent(:,i) /= 0.0d0) .or. & - tp_coarsen(i) > 1) then + if (any(topo_crop_extent(:,i) /= 0.0d0) .or. & + topo_coarsen(i) > 1) then ! Base window: indices of file points inside the ! crop extent (first >= lower edge, last <= upper ! edge, with a tolerance so grid-aligned extents ! stay exact), or the full file if only coarsening. - ! tp_crop_extent is in domain coordinates; subtract - ! tp_x_shift to get file coordinates (the header + ! topo_crop_extent is in domain coordinates; subtract + ! topo_x_shift to get file coordinates (the header ! values are in file coords at this point, before - ! tp_x_shift is applied below). - if (any(tp_crop_extent(:,i) /= 0.0d0)) then - x1_c = tp_crop_extent(1,i) - tp_x_shift(i) - x2_c = tp_crop_extent(2,i) - tp_x_shift(i) - y1_c = tp_crop_extent(3,i) - tp_y_shift(i) - y2_c = tp_crop_extent(4,i) - tp_y_shift(i) + ! topo_x_shift is applied below). + if (any(topo_crop_extent(:,i) /= 0.0d0)) then + x1_c = topo_crop_extent(1,i) - topo_x_shift(i) + x2_c = topo_crop_extent(2,i) - topo_x_shift(i) + y1_c = topo_crop_extent(3,i) - topo_y_shift(i) + y2_c = topo_crop_extent(4,i) - topo_y_shift(i) i_start = max(0, ceiling((x1_c - xlowtopo(i)) & / dxtopo(i) - 1.0d-6)) i_end = min(mxtopo(i)-1, & @@ -337,9 +337,9 @@ subroutine read_topo_settings(restart,file_name) floor((y2_c - ylowtopo(i)) & / dytopo(i) + 1.0d-6)) if (i_start > i_end .or. j_start > j_end) then - print *, "ERROR: tp_crop_extent does not overlap topo file:" + print *, "ERROR: topo_crop_extent does not overlap topo file:" print *, " file = ", trim(topofname(i)) - print *, " crop = ", tp_crop_extent(:,i) + print *, " crop = ", topo_crop_extent(:,i) stop 1 end if else @@ -355,43 +355,43 @@ subroutine read_topo_settings(restart,file_name) ! format written by TopographyData.write()). call apply_align_buffer_coarsen(i_start, i_end, & mxtopo(i), xlowtopo(i), dxtopo(i), & - tp_align(1,i) - tp_x_shift(i), & - any(tp_align(:,i) /= 0.0d0), & - max(1, tp_coarsen(i)), tp_buffer(i)) + topo_align(1,i) - topo_x_shift(i), & + any(topo_align(:,i) /= 0.0d0), & + max(1, topo_coarsen(i)), topo_buffer(i)) call apply_align_buffer_coarsen(j_start, j_end, & mytopo(i), ylowtopo(i), dytopo(i), & - tp_align(2,i) - tp_y_shift(i), & - any(tp_align(:,i) /= 0.0d0), & - max(1, tp_coarsen(i)), tp_buffer(i)) + topo_align(2,i) - topo_y_shift(i), & + any(topo_align(:,i) /= 0.0d0), & + max(1, topo_coarsen(i)), topo_buffer(i)) ! Update stored extents/sizes to the kept lattice - ! (still file coords; tp_x_shift is applied below). + ! (still file coords; topo_x_shift is applied below). xlowtopo(i) = xlowtopo(i) + real(i_start,8)*dxtopo(i) xhitopo(i) = xlowtopo(i) + real(i_end-i_start,8)*dxtopo(i) ylowtopo(i) = ylowtopo(i) + real(j_start,8)*dytopo(i) yhitopo(i) = ylowtopo(i) + real(j_end-j_start,8)*dytopo(i) - mxtopo(i) = (i_end - i_start)/max(1, tp_coarsen(i)) + 1 - mytopo(i) = (j_end - j_start)/max(1, tp_coarsen(i)) + 1 - dxtopo(i) = dxtopo(i)*max(1, tp_coarsen(i)) - dytopo(i) = dytopo(i)*max(1, tp_coarsen(i)) - tp_i_start(i) = i_start - tp_j_start(i) = j_start + mxtopo(i) = (i_end - i_start)/max(1, topo_coarsen(i)) + 1 + mytopo(i) = (j_end - j_start)/max(1, topo_coarsen(i)) + 1 + dxtopo(i) = dxtopo(i)*max(1, topo_coarsen(i)) + dytopo(i) = dytopo(i)*max(1, topo_coarsen(i)) + topo_i_start(i) = i_start + topo_j_start(i) = j_start end if end if - ! Apply tp_x_shift to the stored domain extents so that all + ! Apply topo_x_shift to the stored domain extents so that all ! coordinate comparisons (topoarea, rectintegral, ! set_topo_for_dtopo) use the shifted bounds. For types 2/3 ! this is the only place x_shift is applied; for type 4 we ! also shift xlocs inside read_topo_file so the xstart ! lookup against xlowtopo (= xll) still finds the right index. - if (tp_x_shift(i) /= 0.0d0) then - xlowtopo(i) = xlowtopo(i) + tp_x_shift(i) - xhitopo(i) = xhitopo(i) + tp_x_shift(i) + if (topo_x_shift(i) /= 0.0d0) then + xlowtopo(i) = xlowtopo(i) + topo_x_shift(i) + xhitopo(i) = xhitopo(i) + topo_x_shift(i) end if - if (tp_y_shift(i) /= 0.0d0) then - ylowtopo(i) = ylowtopo(i) + tp_y_shift(i) - yhitopo(i) = yhitopo(i) + tp_y_shift(i) + if (topo_y_shift(i) /= 0.0d0) then + ylowtopo(i) = ylowtopo(i) + topo_y_shift(i) + yhitopo(i) = yhitopo(i) + topo_y_shift(i) end if topoID(i) = i @@ -838,8 +838,8 @@ subroutine read_topo_file(mx,my,topo_type,fname,xll,yll,topo,topo_idx) missing = 0 if (present(topo_idx)) then - crop_active = any(tp_crop_extent(:,topo_idx) /= 0.0d0) & - .or. tp_coarsen(topo_idx) > 1 + crop_active = any(topo_crop_extent(:,topo_idx) /= 0.0d0) & + .or. topo_coarsen(topo_idx) > 1 else crop_active = .false. end if @@ -848,8 +848,8 @@ subroutine read_topo_file(mx,my,topo_type,fname,xll,yll,topo,topo_idx) ! ---------------------------------------------------------- ! Crop and/or coarsen active: read the full file into a ! temporary array, then extract the kept lattice computed - ! by read_topo_settings (tp_i_start/tp_j_start, stride - ! tp_coarsen) into the output topo. + ! by read_topo_settings (topo_i_start/topo_j_start, stride + ! topo_coarsen) into the output topo. ! ---------------------------------------------------------- mx_full = mxtopo_full(topo_idx) my_full = mytopo_full(topo_idx) @@ -880,17 +880,17 @@ subroutine read_topo_file(mx,my,topo_type,fname,xll,yll,topo,topo_idx) end select ! Extract the kept mx-by-my lattice from topo_full. The - ! output topo rows run N-to-S while tp_j_start counts from + ! output topo rows run N-to-S while topo_j_start counts from ! the south, so output row jcrop corresponds to the kept - ! south-based row tp_j_start + (my - jcrop)*c, which is + ! south-based row topo_j_start + (my - jcrop)*c, which is ! file row my_full - that (file rows are 1-based from the ! north). - cf = max(1, tp_coarsen(topo_idx)) + cf = max(1, topo_coarsen(topo_idx)) do jcrop = 1, int(my, 8) - j_in_full = int(my_full,8) - (int(tp_j_start(topo_idx),8) & + j_in_full = int(my_full,8) - (int(topo_j_start(topo_idx),8) & + (int(my,8) - jcrop)*int(cf,8)) do icrop = 1, int(mx, 8) - i_in_full = int(tp_i_start(topo_idx),8) & + i_in_full = int(topo_i_start(topo_idx),8) & + (icrop - 1)*int(cf,8) + 1 topo((jcrop-1)*int(mx,8) + icrop) = & topo_full((j_in_full-1)*int(mx_full,8) + i_in_full) @@ -990,7 +990,7 @@ subroutine read_topo_file(mx,my,topo_type,fname,xll,yll,topo,topo_idx) call check_netcdf_error(nf90_get_var(nc_file, y_var_id, ylocs, & start=(/ 1 /), count=(/ my_tot /))) - ! Apply lon_wrap_offset then tp_x_shift to convert file coordinates + ! Apply lon_wrap_offset then topo_x_shift to convert file coordinates ! to domain coordinates, so xstart matches the domain-coord xll ! that was computed by read_topo_header (which also has both ! offsets applied via read_topo_settings). @@ -998,11 +998,11 @@ subroutine read_topo_file(mx,my,topo_type,fname,xll,yll,topo,topo_idx) if (nc_lon_wrap_offset(topo_idx) /= 0.0d0) then xlocs = xlocs + nc_lon_wrap_offset(topo_idx) end if - if (tp_x_shift(topo_idx) /= 0.0d0) then - xlocs = xlocs + tp_x_shift(topo_idx) + if (topo_x_shift(topo_idx) /= 0.0d0) then + xlocs = xlocs + topo_x_shift(topo_idx) end if - if (tp_y_shift(topo_idx) /= 0.0d0) then - ylocs = ylocs + tp_y_shift(topo_idx) + if (topo_y_shift(topo_idx) /= 0.0d0) then + ylocs = ylocs + topo_y_shift(topo_idx) end if end if @@ -1010,7 +1010,7 @@ subroutine read_topo_file(mx,my,topo_type,fname,xll,yll,topo,topo_idx) ! nf90_get_var, so mx/my count kept points and the spacing ! between them is cf times the file spacing. if (present(topo_idx)) then - cf = max(1, tp_coarsen(topo_idx)) + cf = max(1, topo_coarsen(topo_idx)) else cf = 1 end if @@ -1180,16 +1180,16 @@ subroutine read_topo_file(mx,my,topo_type,fname,xll,yll,topo,topo_idx) if (present(topo_idx) .and. abs(topo_type) /= 1) then ! 1. negate_z: explicit Z sign flip, independent of topo_type sign. - if (tp_negate_z(topo_idx)) then + if (topo_negate_z(topo_idx)) then topo(1:mtot) = -topo(1:mtot) end if ! 2. z_shift: bulk datum adjustment. Skip fill/missing cells. - if (tp_z_shift(topo_idx) /= 0.0d0) then + if (topo_z_shift(topo_idx) /= 0.0d0) then do i = 1, mtot if (abs(topo(i) - topo_missing) > & 1.0d-6 * max(abs(topo_missing), 1.0d0)) then - topo(i) = topo(i) + tp_z_shift(topo_idx) + topo(i) = topo(i) + topo_z_shift(topo_idx) end if end do end if @@ -1470,7 +1470,7 @@ subroutine read_topo_header(fname,topo_type,mx,my,xll,yll,xhi,yhi,dx,dy,topo_idx ! ---------------------------------------------------------------- ! Select the subset of the file to load. - ! Priority: nc_crop_bounds (file coords) > tp_crop_extent (domain + ! Priority: nc_crop_bounds (file coords) > topo_crop_extent (domain ! coords, converted to file coords) > AMR domain fallback. ! ---------------------------------------------------------------- nbuf4 = 0 @@ -1482,21 +1482,21 @@ subroutine read_topo_header(fname,topo_type,mx,my,xll,yll,xhi,yhi,dx,dy,topo_idx y_in_dom = (ylocs >= nc_crop_bounds(3, topo_idx)) .and. & (ylocs <= nc_crop_bounds(4, topo_idx)) else if (present(topo_idx) .and. & - any(tp_crop_extent(:,topo_idx) /= 0.0d0)) then - ! tp_crop_extent is in domain coordinates (after nc_lon_wrap_offset - ! and tp_x_shift). Convert to file coordinates for comparison + any(topo_crop_extent(:,topo_idx) /= 0.0d0)) then + ! topo_crop_extent is in domain coordinates (after nc_lon_wrap_offset + ! and topo_x_shift). Convert to file coordinates for comparison ! against raw xlocs (which have neither offset applied yet). - ! tp_buffer is applied in index space below, mirroring + ! topo_buffer is applied in index space below, mirroring ! Python Topography.crop(). - x1_c = tp_crop_extent(1,topo_idx) - nc_lon_wrap_offset(topo_idx) & - - tp_x_shift(topo_idx) - x2_c = tp_crop_extent(2,topo_idx) - nc_lon_wrap_offset(topo_idx) & - - tp_x_shift(topo_idx) - y1_c = tp_crop_extent(3,topo_idx) - tp_y_shift(topo_idx) - y2_c = tp_crop_extent(4,topo_idx) - tp_y_shift(topo_idx) + x1_c = topo_crop_extent(1,topo_idx) - nc_lon_wrap_offset(topo_idx) & + - topo_x_shift(topo_idx) + x2_c = topo_crop_extent(2,topo_idx) - nc_lon_wrap_offset(topo_idx) & + - topo_x_shift(topo_idx) + y1_c = topo_crop_extent(3,topo_idx) - topo_y_shift(topo_idx) + y2_c = topo_crop_extent(4,topo_idx) - topo_y_shift(topo_idx) x_in_dom = (xlocs >= x1_c) .and. (xlocs <= x2_c) y_in_dom = (ylocs >= y1_c) .and. (ylocs <= y2_c) - nbuf4 = tp_buffer(topo_idx) + nbuf4 = topo_buffer(topo_idx) else x_in_dom = (xlocs > (xlower - dx - hxposs(1)*nghost)) .and. & (xlocs < (xupper + dx + hxposs(1)*nghost)) @@ -1529,13 +1529,13 @@ subroutine read_topo_header(fname,topo_type,mx,my,xll,yll,xhi,yhi,dx,dy,topo_idx ! Python Topography.crop(). The masks select contiguous ! runs, so first/last true bound the window. if (present(topo_idx)) then - cf4 = max(1, tp_coarsen(topo_idx)) - has_al4 = any(tp_align(:,topo_idx) /= 0.0d0) + cf4 = max(1, topo_coarsen(topo_idx)) + has_al4 = any(topo_align(:,topo_idx) /= 0.0d0) ! Align targets are in domain coordinates; xlocs now - ! include lon_wrap_offset but not tp_x_shift, so shift the + ! include lon_wrap_offset but not topo_x_shift, so shift the ! x target into the same frame. - al_x4 = tp_align(1,topo_idx) - tp_x_shift(topo_idx) - al_y4 = tp_align(2,topo_idx) - tp_y_shift(topo_idx) + al_x4 = topo_align(1,topo_idx) - topo_x_shift(topo_idx) + al_y4 = topo_align(2,topo_idx) - topo_y_shift(topo_idx) else cf4 = 1 has_al4 = .false. @@ -1666,8 +1666,8 @@ subroutine read_dtopo_settings(file_name) real(kind=8) :: area_i,area_j integer :: i,j ! Unsupported preprocessing values: read, guarded, then discarded - real(kind=8) :: dtp_crop(4), dtp_align_v(2) - integer :: dtp_coarsen_v, dtp_buffer_v + real(kind=8) :: dtopo_crop(4), dtopo_align_v(2) + integer :: dtopo_coarsen_v, dtopo_buffer_v write(GEO_PARM_UNIT,*) ' ' write(GEO_PARM_UNIT,*) '--------------------------------------------' @@ -1698,12 +1698,12 @@ subroutine read_dtopo_settings(file_name) allocate(index0_dtopowork1(num_dtopo),index0_dtopowork2(num_dtopo)) allocate(tdtopo1(num_dtopo),tdtopo2(num_dtopo),taudtopo(num_dtopo)) allocate(mdtopoorder(num_dtopo),topoaltered(num_dtopo)) - allocate(dtp_x_shift(num_dtopo),dtp_y_shift(num_dtopo)) - allocate(dtp_z_shift(num_dtopo),dtp_negate_z(num_dtopo)) - dtp_x_shift = 0.0d0 - dtp_y_shift = 0.0d0 - dtp_z_shift = 0.0d0 - dtp_negate_z = .false. + allocate(dtopo_x_shift(num_dtopo),dtopo_y_shift(num_dtopo)) + allocate(dtopo_z_shift(num_dtopo),dtopo_negate_z(num_dtopo)) + dtopo_x_shift = 0.0d0 + dtopo_y_shift = 0.0d0 + dtopo_z_shift = 0.0d0 + dtopo_negate_z = .false. allocate(dnc_var_name(num_dtopo),dnc_x_name(num_dtopo)) allocate(dnc_y_name(num_dtopo),dnc_dim_order(num_dtopo)) @@ -1728,17 +1728,17 @@ subroutine read_dtopo_settings(file_name) ! topo.data, written by DTopoData.write()). x_shift, y_shift, ! z_shift and negate_z are implemented for dtopography; fail ! loudly on the rest rather than silently ignoring them. - read(iunit,*) dtp_crop(1), dtp_crop(2), dtp_crop(3), dtp_crop(4) - read(iunit,*) dtp_coarsen_v - read(iunit,*) dtp_buffer_v - read(iunit,*) dtp_align_v(1), dtp_align_v(2) - read(iunit,*) dtp_x_shift(i) - read(iunit,*) dtp_y_shift(i) - read(iunit,*) dtp_z_shift(i) - read(iunit,*) dtp_negate_z(i) - - if (any(dtp_crop /= 0.0d0) .or. dtp_coarsen_v > 1 .or. & - dtp_buffer_v /= 0 .or. any(dtp_align_v /= 0.0d0)) then + read(iunit,*) dtopo_crop(1), dtopo_crop(2), dtopo_crop(3), dtopo_crop(4) + read(iunit,*) dtopo_coarsen_v + read(iunit,*) dtopo_buffer_v + read(iunit,*) dtopo_align_v(1), dtopo_align_v(2) + read(iunit,*) dtopo_x_shift(i) + read(iunit,*) dtopo_y_shift(i) + read(iunit,*) dtopo_z_shift(i) + read(iunit,*) dtopo_negate_z(i) + + if (any(dtopo_crop /= 0.0d0) .or. dtopo_coarsen_v > 1 .or. & + dtopo_buffer_v /= 0 .or. any(dtopo_align_v /= 0.0d0)) then print *, "ERROR: only x_shift, y_shift, z_shift and negate_z" print *, " preprocessing are implemented for dtopography files." print *, " Unsupported attribute set for file: ", trim(dtopofname(i)) @@ -1764,13 +1764,13 @@ subroutine read_dtopo_settings(file_name) ! Apply coordinate registration shifts to the stored domain ! extents (domain = file + shift), mirroring the topo handling ! and the Python DTopography.read() x/y_shift application. - if (dtp_x_shift(i) /= 0.0d0) then - xlowdtopo(i) = xlowdtopo(i) + dtp_x_shift(i) - xhidtopo(i) = xhidtopo(i) + dtp_x_shift(i) + if (dtopo_x_shift(i) /= 0.0d0) then + xlowdtopo(i) = xlowdtopo(i) + dtopo_x_shift(i) + xhidtopo(i) = xhidtopo(i) + dtopo_x_shift(i) end if - if (dtp_y_shift(i) /= 0.0d0) then - ylowdtopo(i) = ylowdtopo(i) + dtp_y_shift(i) - yhidtopo(i) = yhidtopo(i) + dtp_y_shift(i) + if (dtopo_y_shift(i) /= 0.0d0) then + ylowdtopo(i) = ylowdtopo(i) + dtopo_y_shift(i) + yhidtopo(i) = yhidtopo(i) + dtopo_y_shift(i) end if enddo @@ -1819,14 +1819,14 @@ subroutine read_dtopo_settings(file_name) ! Apply preprocessing in-memory (original file unchanged), same ! convention as read_topo_file. No fill-cell guard: deformation ! grids have no no-data convention. - if (dtp_negate_z(i)) then + if (dtopo_negate_z(i)) then dtopowork(i0dtopo(i):i0dtopo(i)+mdtopo(i)-1) = & -dtopowork(i0dtopo(i):i0dtopo(i)+mdtopo(i)-1) end if - if (dtp_z_shift(i) /= 0.0d0) then + if (dtopo_z_shift(i) /= 0.0d0) then dtopowork(i0dtopo(i):i0dtopo(i)+mdtopo(i)-1) = & dtopowork(i0dtopo(i):i0dtopo(i)+mdtopo(i)-1) & - + dtp_z_shift(i) + + dtopo_z_shift(i) end if enddo diff --git a/src/python/geoclaw/data.py b/src/python/geoclaw/data.py index 446f5e965..47300c3ed 100755 --- a/src/python/geoclaw/data.py +++ b/src/python/geoclaw/data.py @@ -155,10 +155,13 @@ def write(self,data_source='setrun.py', out_file='refinement.data'): def _write_preprocessing_block(f, t): - """Write the 7 preprocessing-attribute lines for one topo/dtopo file. + """Write the 8 preprocessing-attribute lines for one topo/dtopo file. + + The lines are, in order: crop_extent, coarsen, buffer, align, x_shift, + y_shift, z_shift, negate_z. Shared by TopographyData.write() (topo.data) and DTopoData.write() - (dtopo.data); Fortran reads the same 7 lines in read_topo_settings and + (dtopo.data); Fortran reads the same 8 lines in read_topo_settings and read_dtopo_settings. *t* is a Topography or DTopography object. Float values use repr (shortest round-trip representation) so coordinates @@ -279,7 +282,11 @@ def _normalize_topofiles(self): raw_type = entry.get('topo_type', None) if raw_type is not None: topo.topo_type = int(raw_type) - # 'extent' in the dict spec maps to 'crop_extent' on Topography + # The legacy dict key 'extent' is an alias for 'crop_extent' + # (the requested crop; see Topography "Region terminology"). + # Note: if a spec supplies BOTH 'extent' and 'crop_extent', the + # 'crop_extent' key wins -- it is applied by the loop below, + # which runs after this alias assignment. if 'extent' in entry: topo.crop_extent = entry['extent'] for attr in ('crop_extent', 'coarsen', 'buffer', 'align', @@ -644,7 +651,7 @@ def write(self, data_source='setrun.py', out_file='dtopo.data'): unsupported = [name for name, is_set in ( ("crop_extent", d.crop_extent is not None), ("coarsen", d.coarsen != 1), - ("buffer", d.buffer != 0.0), + ("buffer", d.buffer != 0), ("align", d.align is not None), ) if is_set] if unsupported: @@ -736,7 +743,7 @@ def _data(line): crop = [float(v) for v in _data(lines[i + 2]).split()] d.crop_extent = None if all(v == 0. for v in crop) else crop d.coarsen = int(_data(lines[i + 3])) - d.buffer = float(_data(lines[i + 4])) + d.buffer = int(_data(lines[i + 4])) # grid-point count align = [float(v) for v in _data(lines[i + 5]).split()] d.align = None if all(v == 0. for v in align) else align d.x_shift = float(_data(lines[i + 6])) diff --git a/src/python/geoclaw/dtopotools.py b/src/python/geoclaw/dtopotools.py index 50da0a6ac..cf451a1cd 100644 --- a/src/python/geoclaw/dtopotools.py +++ b/src/python/geoclaw/dtopotools.py @@ -306,7 +306,7 @@ def __init__(self, path=None, dtopo_type=None, time_reference=None): # others raise NotImplementedError if set (see read()). self.crop_extent = None # [x1,x2,y1,y2]; None = full domain self.coarsen = 1 - self.buffer = 0.0 + self.buffer = 0 # grid-point count (see Topography.crop) self.align = None self.x_shift = 0.0 self.y_shift = 0.0 @@ -360,7 +360,7 @@ def read(self, path=None, dtopo_type=None, verbose=False, unsupported = [name for name, is_set in ( ("crop_extent", self.crop_extent is not None), ("coarsen", self.coarsen != 1), - ("buffer", self.buffer != 0.0), + ("buffer", self.buffer != 0), ("align", self.align is not None), ) if is_set] if unsupported: diff --git a/src/python/geoclaw/etopotools.py b/src/python/geoclaw/etopotools.py index 54809a3b6..3fac15e43 100644 --- a/src/python/geoclaw/etopotools.py +++ b/src/python/geoclaw/etopotools.py @@ -1,24 +1,71 @@ """ -Tools to download etopo topography/bathymetry data from NCEI (formerly NGDC). -See http://www.ngdc.noaa.gov/mgg/global/global.html +Tools to download ETOPO topography/bathymetry data from NCEI (formerly NGDC). +See https://www.ncei.noaa.gov/products/etopo-global-relief-model -Note the new etopo1_download_nc is better to use than etopo1_download. +Two entry points are provided: + +- :func:`fetch_etopo` -- the recommended path. Reads an ETOPO netCDF DEM + (ETOPO 2022 by default) into a :class:`~clawpack.geoclaw.topotools.Topography` + via :func:`clawpack.geoclaw.topotools.fetch_remote_topo`. This is the + consistently-available, near-best-available-data source. + +- :func:`etopo1_download` -- legacy. Downloads a topo_type 3 (ASCII) file from + the old NGDC WCS-proxy endpoint. That endpoint is legacy and often flaky; + prefer :func:`fetch_etopo` (or + :func:`clawpack.geoclaw.topotools.fetch_remote_topo`) instead. """ -from __future__ import absolute_import -from __future__ import print_function +def fetch_etopo(name='etopo22_30sec', crop_extent=None, coarsen=1, buffer=0, + align=None, nc_params={}, verbose=False): + r"""Fetch an ETOPO netCDF DEM as a `Topography`. + + Thin convenience wrapper over + :func:`clawpack.geoclaw.topotools.fetch_remote_topo` for the ETOPO netCDF + nicknames in ``topotools.remote_topo_urls`` (e.g. the default + ``'etopo22_30sec'`` = ETOPO 2022 30 arcsecond, or ``'etopo1'``). + + :Input: + + - *name* (str) - nickname (key of ``topotools.remote_topo_urls``) or a URL + to an ETOPO netCDF file. Default ``'etopo22_30sec'``. + - *crop_extent* ([x1, x2, y1, y2] or None) - requested crop in domain + coordinates; ``None`` reads the whole file. + - *coarsen* (int) - coarsening factor (1 = none). + - *buffer* (int) - points to keep outside the crop on each side. + - *align* (tuple) - alignment when coarsening; see ``Topography.crop``. + - *nc_params* (dict) - options forwarded to the ``topo_type=4`` reader + (e.g. ``z_var``, ``assume_units``); see ``Topography.read``. + - *verbose* (bool) - if True, print the resolved source. + + :Output: + + - a :class:`~clawpack.geoclaw.topotools.Topography` object. + """ + + from clawpack.geoclaw import topotools + + return topotools.fetch_remote_topo(name, crop_extent=crop_extent, + coarsen=coarsen, buffer=buffer, + align=align, nc_params=nc_params, + verbose=verbose) -def etopo1_download(xlimits, ylimits, dx=0.0166666666667, dy=None, \ - output_dir='.', file_name=None, force=False, verbose=True, \ - return_topo=False): + +def etopo1_download(xlimits, ylimits, dx=0.0166666666667, dy=None, + output_dir='.', file_name=None, force=False, verbose=True, + return_topo=None): """ - Create a url to download etopo1 topography from NCEI and - save as a topo_type 3 file. Uses the database described at - http://www.ngdc.noaa.gov/mgg/global/global.html + Download etopo1 topography from NCEI and save as a topo_type 3 file, then + return it as a `Topography` object. + + .. note:: + This uses the old NGDC WCS-proxy endpoint, which is legacy and often + flaky. For a consistently-available, modern netCDF source prefer + :func:`fetch_etopo` or + :func:`clawpack.geoclaw.topotools.fetch_remote_topo`. :Inputs: @@ -38,8 +85,14 @@ def etopo1_download(xlimits, ylimits, dx=0.0166666666667, dy=None, \ - *file_name*: name of file, default is constructed from xlimits,ylimits - *force*: if True, download even if the file already exists. - *verbose*: if True, print info from clawpack.clawutil.data.get_remote_file + - *return_topo*: deprecated and ignored; a `Topography` is always returned. - Note: New NGDC format gives cell-registered values, so shift the + :Output: + + - a :class:`~clawpack.geoclaw.topotools.Topography` object read from the + downloaded topo_type 3 file. + + Note: New NGDC format gives cell-registered values, so shift the values `xllcorner` and `yllcorner` to the specified corner. **To do:** Check whether it is possible to specify grid-registered @@ -49,11 +102,18 @@ def etopo1_download(xlimits, ylimits, dx=0.0166666666667, dy=None, \ so add this in too. """ - from clawpack.geoclaw import util, topotools + from clawpack.geoclaw import topotools from clawpack.clawutil.data import get_remote_file import os + import warnings from numpy import round + if return_topo is not None: + warnings.warn( + "etopo1_download's return_topo argument is deprecated and ignored; " + "a Topography object is now always returned.", + DeprecationWarning, stacklevel=2) + format = '&format=aaigrid' # topo_type 3 if dy is None: @@ -103,7 +163,8 @@ def etopo1_download(xlimits, ylimits, dx=0.0166666666667, dy=None, \ x1 = x1 - longitude_shift # shift back before writing header - lines = open(file_path).readlines() + with open(file_path) as f: + lines = f.readlines() if lines[2].split()[0] != 'xllcorner': print("*** Error downloading, check the file!") else: @@ -114,13 +175,10 @@ def etopo1_download(xlimits, ylimits, dx=0.0166666666667, dy=None, \ if 'nodata_value' not in lines[5]: lines = lines[:5] + ['nodata_value -99999\n'] + lines[5:] print("Added nodata_value line") - f = open(file_path,'w') - f.writelines(lines) - f.close() + with open(file_path, 'w') as f: + f.writelines(lines) print("Created file: ",file_path) - if return_topo: - topo = topotools.Topography() - topo.read(file_path, topo_type=3) - return topo - + topo = topotools.Topography() + topo.read(file_path, topo_type=3) + return topo diff --git a/src/python/geoclaw/netcdf_utils.py b/src/python/geoclaw/netcdf_utils.py index 76fa611fa..323558486 100644 --- a/src/python/geoclaw/netcdf_utils.py +++ b/src/python/geoclaw/netcdf_utils.py @@ -38,6 +38,7 @@ import contextlib import dataclasses import importlib.util +import re import warnings from pathlib import Path from typing import Optional @@ -404,7 +405,17 @@ def __init__( path: str | Path, crop_bounds: Optional[tuple[float, float, float, float]] = None, ) -> None: - self.path = Path(path) + # A remote OPeNDAP/THREDDS URL (e.g. "https://.../foo.nc") must reach + # xarray as a string. Wrapping it in pathlib.Path collapses "https://" + # to "https:/" and makes it a *relative* path, which the netCDF4 backend + # then resolves against the cwd -- producing a bogus local-file lookup + # (PR #726). The scheme-anchored regex ignores Windows drive paths + # like "C:\\..." (no "//"). + if isinstance(path, str) and re.match(r"^[A-Za-z][A-Za-z0-9+.-]*://", + path): + self.path = path + else: + self.path = Path(path) self.crop_bounds = crop_bounds # Activate Dask-lazy chunking if dask is available; fall back to # netCDF4 native lazy loading so dask is an optional dependency. diff --git a/src/python/geoclaw/topotools.py b/src/python/geoclaw/topotools.py index d8c13180f..a6343f4ba 100644 --- a/src/python/geoclaw/topotools.py +++ b/src/python/geoclaw/topotools.py @@ -91,43 +91,98 @@ def determine_topo_type(path, default=None): return topo_type -def _netcdf_window_indices(coords, lo, hi, margin, n, stride=1): - r"""Half-open index window ``[i0, i1)`` into 1-D monotonic *coords*. +def _crop_indices(x, y, crop_extent, coarsen, buffer, align): + r"""Half-open index bounds for a crop+coarsen+align window. - Used by :meth:`Topography.read` to push a ``crop_extent`` down to the - NetCDF read so only the needed hyperslab is loaded from disk (rather than - materializing a whole global variable and cropping afterward). + Returns ``(ilower, iupper, jlower, jupper)`` to be sliced as + ``arr[jlower:jupper:coarsen, ilower:iupper:coarsen]`` (and the coordinate + subsets ``x[ilower:iupper:coarsen]``, ``y[jlower:jupper:coarsen]``), or + ``None`` if *crop_extent* does not overlap the data. + + This is the shared coarsen/align index arithmetic used by both + :meth:`Topography.crop` (on the in-memory arrays) and the ``topo_type=4`` + read (on the cheap 1-D NetCDF coordinate arrays), so ASCII and NetCDF reads + of the same data produce identical grids. :Input: - - *coords* (ndarray) - 1-D coordinate array, monotonic ascending **or** - descending (as stored in the file). - - *lo*, *hi* (float) - requested inclusive coordinate bounds. - - *margin* (int) - extra points to keep on each side so that a subsequent - :meth:`Topography.crop` still has every point it needs (its own buffer - plus the ``coarsen`` alignment search). - - *n* (int) - length of *coords*. - - *stride* (int) - read stride; the low index is snapped **down** to a - multiple of *stride* so the strided sub-window shares the same phase as - striding the full array, keeping the sampled grid identical. - - Returns ``(0, n)`` (the full range) when the interval does not overlap - *coords*, mirroring ``crop()``'s fall-back of leaving the array uncropped - when the filter region misses the topography. + - *x*, *y* (ndarray) - 1-D coordinate arrays, **ascending** (precondition; + this routine contains no N->S/E->W flip logic). + - *crop_extent* ([x1, x2, y1, y2]) - requested crop in the same coords. + - *coarsen* (int) - subsampling factor; ``align`` only has effect when + ``coarsen > 1``. + - *buffer* (int) - grid points to keep outside *crop_extent* on each side + (expanded by ``buffer*coarsen`` native points, as in :meth:`crop`). + - *align* ((xalign, yalign) or None) - desired alignment when coarsening; + ``None`` means no phase snap (start at the crop window). + + ``coarsen`` and ``buffer`` are assumed already ``int()``-coerced. + """ + # dx/dy computed the same way as the `delta` property (round to 15 places), + # so the align fractional-offset search matches crop() bit-for-bit. + dx = numpy.round(abs(x[1] - x[0]), 15) + dy = numpy.round(abs(y[1] - y[0]), 15) + dx_new = dx * coarsen + dy_new = dy * coarsen + + # Find indices of the arrays inside crop_extent: + try: + ilower = (x >= crop_extent[0]).nonzero()[0][0] + iupper = (x <= crop_extent[1]).nonzero()[0][-1] + jlower = (y >= crop_extent[2]).nonzero()[0][0] + jupper = (y <= crop_extent[3]).nonzero()[0][-1] + except IndexError: + # crop_extent does not overlap the data + return None + + # Shift indices if needed for alignment (matches crop() lines historically + # at 2085-2099: pick the low index whose coord best lands on `align`). + if (coarsen > 1) and (align is not None): + xs = numpy.array([x[ilower + i] for i in range(coarsen)]) + offsets = (xs - align[0]) / dx_new + offsets_frac = offsets - numpy.round(offsets) + ioffset = numpy.argmin(abs(offsets_frac)) + ilower = ilower + ioffset + iupper = iupper - numpy.remainder(iupper - ilower, coarsen) + + ys = numpy.array([y[jlower + j] for j in range(coarsen)]) + offsets = (ys - align[1]) / dy_new + offsets_frac = offsets - numpy.round(offsets) + joffset = numpy.argmin(abs(offsets_frac)) + jlower = jlower + joffset + jupper = jupper - numpy.remainder(jupper - jlower, coarsen) + + # buffer, checking limits of arrays: + ilower = numpy.maximum(0, ilower - buffer * coarsen) + jlower = numpy.maximum(0, jlower - buffer * coarsen) + iupper = numpy.minimum(len(x) - 1, iupper + buffer * coarsen) + 1 + jupper = numpy.minimum(len(y) - 1, jupper + buffer * coarsen) + 1 + + return int(ilower), int(iupper), int(jlower), int(jupper) + + +def _axis_file_slice(coord_full, descending, lo, hi, step, n): + r"""Map an ascending window ``[lo:hi:step]`` to a positive-stride slice. + + Used by the ``topo_type=4`` read: :func:`_crop_indices` returns bounds into + an *ascending* view of a file axis, but NetCDF/xarray lazy indexing requires + a **positive** step into the *file-order* axis. For an axis stored + descending (e.g. latitude N->S), the ascending sample indices + ``lo, lo+step, ..., lo+(m-1)*step`` map to file indices ``n-1-(that)``, whose + minimum is ``f0``; reading ``slice(f0, f0+m*step, step)`` returns those same + ``m`` samples in file (descending) order, to be flipped to ascending in + memory afterward. + + Returns ``(file_slice, coord_subset, flip)`` where + ``coord_subset == coord_full[file_slice]`` and ``flip`` is True when the + subset (and the corresponding data axis) must be reversed to be ascending. """ - # A monotonic array clipped to [lo, hi] yields a contiguous True block for - # either sort order, so first/last True index bound the window. - mask = (coords >= lo) & (coords <= hi) - idx = numpy.nonzero(mask)[0] - if idx.size == 0: - return 0, n - i0 = max(0, int(idx[0]) - margin) - i1 = min(n, int(idx[-1]) + margin + 1) - # Snap the low index down to a multiple of stride so coords[i0:i1:stride] - # is a phase-aligned subset of coords[::stride]. - i0 -= i0 % stride - if i1 <= i0: - return 0, n - return i0, i1 + if not descending: + sl = slice(lo, hi, step) + return sl, coord_full[sl], False + m = len(range(lo, hi, step)) + f0 = n - 1 - (lo + (m - 1) * step) + sl = slice(f0, f0 + m * step, step) + return sl, coord_full[sl], True def create_topo_func(loc,verbose=False): @@ -313,6 +368,49 @@ def swapheader(inputfile, outputfile): +# Sentinel for the deprecated crop-region kwargs (filter_region=, extent=). +# Their real default (None) is itself a valid user value, so a distinct sentinel +# is needed to tell "not passed" apart from "passed None". +_CROP_EXTENT_UNSET = object() + +# Sentinel for read()'s align= kwarg. align=None is a valid user value ("no +# phase snap"), and callers such as fetch_remote_topo set the self.align +# attribute *before* calling read(); a distinct sentinel lets read() tell "not +# passed" (leave self.align alone) apart from "passed None" (override it). +_ALIGN_UNSET = object() + + +def _resolve_crop_extent(crop_extent, deprecated): + r"""Fold deprecated region kwargs onto ``crop_extent``. + + The requested-crop rectangle has one canonical name, ``crop_extent`` (see the + "Region terminology" section of :class:`Topography`). Older APIs spelled it + ``filter_region`` or ``extent``; this helper maps those onto ``crop_extent``. + + :Input: + - *crop_extent* - the value of the canonical ``crop_extent`` argument. + - *deprecated* (dict) - maps each old kwarg name to the value it was called + with, or ``_CROP_EXTENT_UNSET`` if it was not passed. + + For any old name that WAS passed, emit a ``DeprecationWarning`` and use its + value as ``crop_extent``; raise ``TypeError`` if ``crop_extent`` is also + supplied (ambiguous). + """ + import warnings + for name, value in deprecated.items(): + if value is _CROP_EXTENT_UNSET: + continue + if crop_extent is not None: + raise TypeError( + "Got both 'crop_extent' and the deprecated '%s'; " + "pass only 'crop_extent'." % name) + warnings.warn( + "The '%s' argument is deprecated; use 'crop_extent' instead." % name, + DeprecationWarning, stacklevel=3) + crop_extent = value + return crop_extent + + # ============================================================================== # Topography class # ============================================================================== @@ -339,6 +437,35 @@ class Topography(object): >>> topo_file.read('./topo.tt3', topo_type=3) >>> topo_file.plot() + :Region terminology: + + Several attributes/arguments describe rectangular regions, all ordered + ``[x1, x2, y1, y2]`` (x-pair then y-pair). Two axes distinguish them: + *role* (a derived result vs. a requested crop) and *coordinate frame* + (domain vs. file). + + - ``extent`` -- the *actual* bounds of the loaded data, in DOMAIN + coordinates. A read-only, lazily-computed :func:`property` (a *result*, + not an input); see :attr:`extent`. + - ``crop_extent`` -- the *requested* crop rectangle, in DOMAIN coordinates. + This is the single canonical name for the request: it is both a persisted + attribute (default ``None`` = no crop) and the argument name accepted by + :meth:`read`, :meth:`crop`, :meth:`interp_unstructured`, and + :func:`fetch_remote_topo`. ``read(crop_extent=r)`` is equivalent to + setting the attribute and then reading. Mirrors Fortran + ``topo_crop_extent``. (The older argument spellings ``filter_region`` and + ``extent=`` are deprecated aliases.) + - ``crop_bounds`` -- the same requested crop expressed in FILE coordinates. + Used only in the NetCDF (type-4) layer + (``netcdf_utils.FileMetadata.crop_bounds`` -> Fortran ``nc_crop_bounds``); + converted from ``crop_extent`` by subtracting ``lon_wrap_offset``/ + ``x_shift``. + - ``crop()`` -- the operation that turns a ``crop_extent`` request (plus + ``coarsen``/``buffer``/``align``) into a new cropped object. + + Convention: the ``_extent`` suffix denotes domain coordinates; ``_bounds`` + denotes file coordinates. + """ @property @@ -427,7 +554,15 @@ def Y(self): @property def extent(self): - r"""Extent of the topography.""" + r"""Actual bounds of the loaded data, ordered ``[x1, x2, y1, y2]``. + + This is a derived *result* (the min/max of the loaded ``x``/``y``), in + domain coordinates -- not a requested crop; for the crop request see + ``crop_extent`` and the "Region terminology" section of the class + docstring. Computed lazily and cached in ``_extent``; the cache is + invalidated (set to ``None``) whenever ``x``/``y`` change or ``read()`` + reloads/crops the data. + """ if self._extent is None: self._extent = ( numpy.min(self.x), numpy.max(self.x), numpy.min(self.y), numpy.max(self.y) ) @@ -509,28 +644,15 @@ def __init__(self, path=None, topo_type=None, topo_func=None, self.coordinate_transform = lambda x,y: (x,y) # Preprocessing attributes — applied by read() after data is loaded. - # - # BOUNDS TERMINOLOGY (consistent across Python and Fortran): - # * ``extent`` — the *actual* loaded-data bounds. Read-only - # @property (below), in DOMAIN coordinates. - # * ``crop_extent`` — the *requested* crop region, in DOMAIN - # coordinates. Mirrors Fortran ``tp_crop_extent`` - # and the gridded met module's ``crop_extent``. - # * ``crop_bounds`` — the same region in FILE coordinates, carried - # only in the NetCDF (type-4) descriptor - # (``netcdf_utils.FileMetadata.crop_bounds`` -> - # Fortran ``nc_crop_bounds``). The descriptor - # writer converts crop_extent -> crop_bounds by - # subtracting lon_wrap_offset/x_shift. - # Convention: the ``_extent`` suffix is domain coords; ``_bounds`` is - # file coords. All 4-element vectors are ordered [x1, x2, y1, y2] - # (x-pair then y-pair). 'crop_extent' is named to avoid shadowing the - # 'extent' property above. - # PATH NOTE: 'path' already exists as an instance attribute set above. - # No topo_path alias is needed; callers should use self.path. + # See the "Region terminology" section of the class docstring for the + # extent / crop_extent / crop_bounds glossary and the Python<->Fortran + # name mapping. Convention recap: the ``_extent`` suffix is domain + # coords, ``_bounds`` is file coords, all ordered [x1, x2, y1, y2]. + # PATH NOTE: 'path' already exists as an instance attribute set above; + # no topo_path alias is needed, callers should use self.path. self.crop_extent: list[float] | None = None # [x1,x2,y1,y2]; None=full domain self.coarsen: int = 1 - self.buffer: float = 0.0 + self.buffer: int = 0 self.align = None self.x_shift: float = 0.0 self.y_shift: float = 0.0 @@ -698,8 +820,9 @@ def generate_2d_coordinates(self, mask=False): def read(self, path=None, topo_type=None, unstructured=False, - mask=False, filter_region=None, force=False, stride=[1, 1], - nc_params={}): + mask=False, crop_extent=None, force=False, + coarsen=None, align=_ALIGN_UNSET, buffer=None, stride=None, + nc_params={}, filter_region=_CROP_EXTENT_UNSET): r"""Read in the data from the object's *path* attribute. Stores the resulting data in one of the sets of *x*, *y*, and *z* or @@ -711,10 +834,30 @@ def read(self, path=None, topo_type=None, unstructured=False, - *unstructured* (bool) - default is False for lat-long grids. - *mask* (bool) - whether to store as masked array for missing values (default if False) - - *filter_region* (tuple) - - *stride* (list) - List of strides for the x and y dimensions - respectively. Default is *[1, 1]*. Note that this is only - implemented for NetCDF reading currently. + - *crop_extent* ([x1, x2, y1, y2] or None) - requested crop region in + domain coordinates (see the "Region terminology" section of the + class docstring). Passing it here is equivalent to setting the + ``crop_extent`` attribute before calling ``read()``; the crop is + applied (together with ``coarsen``/``buffer``/``align``) via + :meth:`crop`. Default ``None`` = no crop. The older ``filter_region`` + keyword is a deprecated alias. + - *coarsen* (int) - subsampling factor (1 = no coarsening). Applied + identically for ASCII and NetCDF reads. Passing it here is + equivalent to setting the ``coarsen`` attribute before ``read()``. + See :meth:`crop`. + - *align* ((xalign, yalign) or None) - desired alignment when + coarsening; see :meth:`crop`. ``None`` (the default) means **no + phase snap** -- subsampling starts at the crop window, matching + ASCII/:meth:`crop`. (This differs from the old NetCDF ``stride`` + behavior, which snapped to the file's grid origin.) Pass e.g. + ``align=[integer_lon, integer_lat]`` to lock the coarsened grid to a + fixed lattice regardless of the requested ``crop_extent``. + - *buffer* (int) - grid points to keep outside ``crop_extent`` on each + side; see :meth:`crop`. + - *stride* (list or int) - **Deprecated**: use ``coarsen`` instead. + A NetCDF-only knob that silently did nothing for ASCII reads and used + a different alignment convention. A scalar (or equal-valued list) is + mapped onto ``coarsen``; per-axis striding is no longer supported. - *nc_params* (dict) - options for NetCDF (`topo_type=4`) reading: - `z_var` (str): name of the elevation variable, if it cannot be @@ -730,6 +873,54 @@ def read(self, path=None, topo_type=None, unstructured=False, """ + # A crop_extent passed here is equivalent to setting the attribute first; + # fold the deprecated filter_region alias onto it, then store it so the + # single attribute-driven crop below (and the type-4 pushdown) apply it. + crop_extent = _resolve_crop_extent(crop_extent, + {'filter_region': filter_region}) + if crop_extent is not None: + self.crop_extent = crop_extent + + # Fold coarsen/align/buffer args onto the attributes (mirrors crop_extent + # above): passing them to read() is equivalent to setting the attribute + # first. Sentinels distinguish "not passed" from an explicit value so a + # caller that presets self.align/self.coarsen/self.buffer before read() + # (e.g. fetch_remote_topo) is not silently clobbered. + if coarsen is not None: + self.coarsen = int(coarsen) + if buffer is not None: + self.buffer = int(buffer) + if align is not _ALIGN_UNSET: + self.align = align + + # `stride` is deprecated: a NetCDF-only knob that silently did nothing + # for ASCII reads and used a different alignment convention than + # crop()/coarsen. Map it onto the unified scalar `coarsen`. + if stride is not None: + import warnings + warnings.warn( + "The 'stride' argument to Topography.read() is deprecated; use " + "'coarsen' (a scalar subsampling factor) instead. 'coarsen' is " + "applied identically for ASCII and NetCDF reads.", + DeprecationWarning, + stacklevel=2, + ) + if numpy.ndim(stride) == 0: + _stride = int(stride) + else: + _s = list(stride) + if len(_s) == 0 or any(int(v) != int(_s[0]) for v in _s): + raise ValueError( + "Per-axis stride is no longer supported; 'coarsen' is a " + "single scalar factor. Got stride=%r." % (stride,)) + _stride = int(_s[0]) + if _stride != 1: + if self.coarsen != 1 and self.coarsen != _stride: + raise ValueError( + "Pass either 'stride' or 'coarsen', not both " + "(stride=%r, coarsen=%r)." % (stride, self.coarsen)) + self.coarsen = _stride + if (path is None) and (self.path is None): raise ValueError("*** Need to set path for file to read") @@ -765,17 +956,17 @@ def read(self, path=None, topo_type=None, unstructured=False, points = [] values = [] - # Filter region if requested - if filter_region is not None: + # Filter region if requested (crop_extent, domain coords) + if self.crop_extent is not None: for coordinate in data: - if filter_region[0] <= coordinate[0] <= filter_region[1]: - if filter_region[2] <= coordinate[1] <= filter_region[3]: + if self.crop_extent[0] <= coordinate[0] <= self.crop_extent[1]: + if self.crop_extent[2] <= coordinate[1] <= self.crop_extent[3]: points.append(coordinate[0:2]) values.append(coordinate[2]) if len(points) == 0: raise Exception("No points were found inside requested " \ - + "filter region.") + + "crop_extent region.") # Cast lists as ndarrays self._x = numpy.array(points[:,0]) @@ -804,7 +995,7 @@ def read(self, path=None, topo_type=None, unstructured=False, _preprocessing_requested = ( self.crop_extent is not None or self.coarsen != 1 - or self.buffer != 0.0 + or self.buffer != 0 or self.align is not None or self.x_shift != 0.0 or self.y_shift != 0.0 @@ -917,43 +1108,66 @@ def read(self, path=None, topo_type=None, unstructured=False, # Topography. _da = _da.transpose(_y_name, _x_name) - # Push both `stride` and `crop_extent` down to xarray's lazy - # indexing so the NetCDF backend reads ONLY the requested + # Push the crop+coarsen+align window down to xarray's lazy + # indexing so the NetCDF backend reads ONLY the final # hyperslab. Otherwise `_da.values` materializes the whole # variable (a global DEM is many GB, and CF fill decoding # promotes it to float), which is prohibitively slow and can - # exhaust memory even when the caller asked for a small - # subset. The crop window is computed on the cheap 1-D - # coordinate arrays and expanded by a margin so the post-read - # crop() below still reproduces its exact bounds/buffer/ - # align/coarsen result on the in-memory sub-window. + # exhaust memory even when the caller asked for a small, + # coarsened subset. The window is computed on the cheap 1-D + # coordinate arrays via _crop_indices -- the SAME routine + # crop() uses -- so this read matches an ASCII read + crop() + # of the same data exactly; the post-read crop() is then + # skipped for topo_type 4 (the data is already final). _nx = _lon_full.size _ny = _lat_full.size + _c = max(int(self.coarsen), 1) + + # _crop_indices requires ascending coords; build ascending + # views of the file axes (lon is usually ascending, lat is + # often stored N→S). + _lon_desc = _lon_full[0] > _lon_full[-1] + _lat_desc = not _meta.y_increasing + _lon_asc = _lon_full[::-1] if _lon_desc else _lon_full + _lat_asc = _lat_full[::-1] if _lat_desc else _lat_full + if self.crop_extent is not None: - _x1, _x2, _y1, _y2 = self.crop_extent - _margin = (int(self.buffer) + 1) * max(int(self.coarsen), 1) - _i0, _i1 = _netcdf_window_indices( - _lon_full, _x1, _x2, _margin, _nx, stride[0] - ) - _j0, _j1 = _netcdf_window_indices( - _lat_full, _y1, _y2, _margin, _ny, stride[1] - ) + _ce = list(self.crop_extent) else: - _i0, _i1 = 0, _nx - _j0, _j1 = 0, _ny - - _da = _da.isel({ - _y_name: slice(_j0, _j1, stride[1]), - _x_name: slice(_i0, _i1, stride[0]), - }) + # whole file (mirrors crop()'s crop_extent=self.extent) + _ce = [_lon_asc[0], _lon_asc[-1], + _lat_asc[0], _lat_asc[-1]] + + _idx = _crop_indices(_lon_asc, _lat_asc, _ce, _c, + int(self.buffer), self.align) + if _idx is None: + # crop_extent misses the file: fall back to the full grid + # at native resolution (mirrors crop() returning None -> + # no-op), rather than coarsening the whole file. + _il, _iu, _jl, _ju = 0, _nx, 0, _ny + _step = 1 + else: + _il, _iu, _jl, _ju = _idx + _step = _c + + # Map each ascending [lo:hi:step] window to a positive-stride + # slice into the FILE-order axis; descending axes are flipped + # in memory afterward. + _x_slice, _lon_vals, _flip_x = _axis_file_slice( + _lon_full, _lon_desc, _il, _iu, _step, _nx) + _y_slice, _lat_vals, _flip_y = _axis_file_slice( + _lat_full, _lat_desc, _jl, _ju, _step, _ny) + + _da = _da.isel({_y_name: _y_slice, _x_name: _x_slice}) _z_vals = numpy.asarray(_da.values, dtype=float) - _lon_vals = _lon_full[_i0:_i1:stride[0]] - _lat_vals = _lat_full[_j0:_j1:stride[1]] - # Flip to S→N (y increasing) if file stores N→S - if not _meta.y_increasing: + # Flip descending axes to ascending (S→N, W→E) in memory. + if _flip_y: _lat_vals = _lat_vals[::-1] _z_vals = _z_vals[::-1, :] + if _flip_x: + _lon_vals = _lon_vals[::-1] + _z_vals = _z_vals[:, ::-1] # Apply unit conversion if source is not already meters _contract = _NC_UNITS.get('topo', 'm') @@ -1029,21 +1243,6 @@ def read(self, path=None, topo_type=None, unstructured=False, if mask: self._Z = numpy.ma.masked_invalid(self._Z) - # Perform region filtering by delegating to crop() so the index - # bounds are computed in exactly one place and are inclusive of the - # filter_region edges (the previous inline slice dropped the upper - # edge row/column). - if filter_region is not None: - _filtered = self.crop(filter_region=filter_region) - if _filtered is not None: - self._x = _filtered._x - self._y = _filtered._y - self._Z = _filtered._Z - self._X = None - self._Y = None - self._extent = None - self._delta = None - # --------------------------------------------------------------- # Apply preprocessing attributes in-memory (original file unchanged). # Fortran applies the same attributes independently in read_topo_file @@ -1054,7 +1253,10 @@ def read(self, path=None, topo_type=None, unstructured=False, # 3. x_shift (shift x array; Fortran shifts xlowtopo/xhitopo) # 3b. y_shift (shift y array; Fortran shifts ylowtopo/yhitopo) # 4+5. crop + coarsen via self.crop() (Fortran: crop+buffer done, - # coarsen not yet implemented) + # coarsen not yet implemented). SKIPPED for topo_type 4: + # the NetCDF read already applied crop+coarsen+align+buffer + # via _crop_indices while reading the hyperslab, so running + # crop() again would double-coarsen. # Steps are skipped when the attribute equals its default value. # --------------------------------------------------------------- if self.negate_z: @@ -1068,9 +1270,10 @@ def read(self, path=None, topo_type=None, unstructured=False, if self.y_shift != 0.0: self._y = self._y + self.y_shift self._extent = None - if self.crop_extent is not None or self.coarsen > 1: + if abs(self.topo_type) != 4 \ + and (self.crop_extent is not None or self.coarsen > 1): _cropped = self.crop( - filter_region=self.crop_extent, + crop_extent=self.crop_extent, coarsen=int(self.coarsen), buffer=int(self.buffer), align=self.align, @@ -1632,17 +1835,18 @@ def plot(self, axes=None, contour_levels=None, contour_kwargs={}, return axes - def interp_unstructured(self, fill_topo, extent=None, method='nearest', + def interp_unstructured(self, fill_topo, crop_extent=None, method='nearest', delta=None, delta_limit=20.0, no_data_value=-99999, buffer_length=100.0, proximity_radius=100.0, - resolution_limit=2000): + resolution_limit=2000, + extent=_CROP_EXTENT_UNSET): r"""Interpolate unstructured data on to regular grid. Function to interpolate the unstructured data in the topo object onto a structured grid. Utilizes a bounding box plus a buffer of size - *buffer_length* (meters) containing all data unless *extent is not None* - is *True*. Then uses the fill topography *fill_topo* to fill in the + *buffer_length* (meters) containing all data unless *crop_extent is not + None* is *True*. Then uses the fill topography *fill_topo* to fill in the gaps in the unstructured data. By default this is done by masking the fill data with the extents, the value *no_data_value* and if *proximity_radius* (meters) is not 0, by a radius of *proximity_radius* @@ -1659,8 +1863,10 @@ def interp_unstructured(self, fill_topo, extent=None, method='nearest', :Input: - *fill_topo* (list) - List of Topography objects to use as fill data in the projection. - - *extent* (tuple) - A tuple defining the rectangle of the sub-section. - Must be in the form (x lower,x upper,y lower, y upper). + - *crop_extent* (tuple) - A tuple defining the rectangle of the + sub-section, in the form (x1, x2, y1, y2). Default ``None`` uses the + data bounding box plus *buffer_length*. The older ``extent`` keyword + is a deprecated alias. - *method* (string) - Method used for interpolation, valid methods are found in *scipy.interpolate.griddata*. Default is *nearest*. - *delta* (tuple) - Directly set the grid spacing of the interpolation @@ -1672,7 +1878,7 @@ def interp_unstructured(self, fill_topo, extent=None, method='nearest', - *no_data_value* (float) - Value to use if no data was found to fill in a missing value, ignored if `method = 'nearest'`. Default is `-99999`. - *buffer_length* (float) - Buffer around bounding box, only applicable - when *extent* is None. Default is `100.0` meters. + when *crop_extent* is None. Default is `100.0` meters. - *proximity_radius* (float) - Radius every unstructured data point used to mask the fill data with. Default is `100.0` meters. - *resolution_limit* (int) - Limit the number of grid points in a @@ -1683,6 +1889,10 @@ def interp_unstructured(self, fill_topo, extent=None, method='nearest', """ + crop_extent = _resolve_crop_extent(crop_extent, {'extent': extent}) + # Internal working name for the interpolation output bounding box. + extent = crop_extent + import scipy.interpolate as interpolate from scipy.spatial import cKDTree @@ -1920,18 +2130,26 @@ def smooth_data(self, indices, r=1): self.Z[index[0], index[1]] = summation / num_points - def crop(self, filter_region=None, coarsen=1, buffer=0, align=None): - r"""Crop region to *filter_region* + def crop(self, crop_extent=None, coarsen=1, buffer=0, align=None, + filter_region=_CROP_EXTENT_UNSET): + r"""Crop region to *crop_extent* Create a new Topography object that is identical to this one but cropped - to the region specified by filter_region + to the region specified by *crop_extent* (see the "Region terminology" + section of the class docstring). :Input: - - *filter_region* (tuple): (x1,x2,y1,y2) desired new extent - - *coarsen* (int): coarsening factor (by subsampling) - - *buffer* (int): when possible, have at least this many points - outside the filter_region on each side + - *crop_extent* (tuple): (x1,x2,y1,y2) desired new extent, in domain + coordinates. Default ``None`` crops to the current ``extent`` (so + only *coarsen* has effect). The older ``filter_region`` keyword is + a deprecated alias. + - *coarsen* (int): coarsening factor (by subsampling). Truncated to + an integer via ``int()``. + - *buffer* (int): integer number of grid points to keep on each side + of *crop_extent* (when possible) -- NOT a coordinate distance (cf. + ``interp_unstructured``'s ``buffer_length``, which is in meters). + Truncated to an integer via ``int()``. - *align* (tuple): (xalign,yalign) = desired alignment if coarsening Setting *buffer > 0* may be useful to insure that the @@ -1973,53 +2191,36 @@ def crop(self, filter_region=None, coarsen=1, buffer=0, align=None): leave the resulting topography as unstructured effectively. """ + crop_extent = _resolve_crop_extent(crop_extent, + {'filter_region': filter_region}) + + # buffer and coarsen are integer grid-point counts (they feed the index + # arithmetic below); truncate any float via int() as documented, rather + # than failing later with an opaque "slice indices must be integers". + buffer = int(buffer) + coarsen = int(coarsen) + if self.unstructured: raise NotImplementedError("*** Cannot currently crop unstructured topo") - if filter_region is None: + if crop_extent is None: # only want to coarsen, so this is entire region: - #filter_region = [self.x[0],self.x[-1],self.y[0],self.y[-1]] - filter_region = self.extent + #crop_extent = [self.x[0],self.x[-1],self.y[0],self.y[-1]] + crop_extent = self.extent - xlower,xupper,ylower,yupper = filter_region + xlower,xupper,ylower,yupper = crop_extent dx,dy = self.delta dx_new = dx*coarsen dy_new = dy*coarsen - # Find indices of topo arrays in filter_region: - try: - ilower = (self.x >= filter_region[0]).nonzero()[0][0] - iupper = (self.x <= filter_region[1]).nonzero()[0][-1] - jlower = (self.y >= filter_region[2]).nonzero()[0][0] - jupper = (self.y <= filter_region[3]).nonzero()[0][-1] - except: - print('*** filter_region does not overlap topo') + # Find crop+coarsen+align index window (shared with the topo_type=4 + # read path so ASCII and NetCDF reads of the same data match exactly). + idx = _crop_indices(self.x, self.y, crop_extent, coarsen, buffer, align) + if idx is None: + print('*** crop_extent does not overlap topo') return None - - # shift indices if needed for alignment: - if (coarsen > 1) and (align is not None): - xs = numpy.array([self.x[ilower + i] for i in range(coarsen)]) - offsets = (xs - align[0]) / dx_new - offsets_frac = offsets - numpy.round(offsets) - ioffset = numpy.argmin(abs(offsets_frac)) - ilower = ilower + ioffset - iupper = iupper - numpy.remainder(iupper-ilower, coarsen) - #print(f'+++ shifted ilower by ioffset={ioffset} to {ilower}') - - ys = numpy.array([self.y[jlower + j] for j in range(coarsen)]) - offsets = (ys - align[1]) / dy_new - offsets_frac = offsets - numpy.round(offsets) - joffset = numpy.argmin(abs(offsets_frac)) - jlower = jlower + joffset - jupper = jupper - numpy.remainder(jupper-jlower, coarsen) - #print(f'+++ shifted jlower by joffset={joffset} to {jlower}') - - # buffer, checking limits of arrays: - ilower = numpy.maximum(0, ilower - buffer*coarsen) - jlower = numpy.maximum(0, jlower - buffer*coarsen) - iupper = numpy.minimum(len(self.x)-1, iupper + buffer*coarsen) + 1 - jupper = numpy.minimum(len(self.y)-1, jupper + buffer*coarsen) + 1 + ilower, iupper, jlower, jupper = idx # Create new topography object: newtopo = Topography() @@ -2182,210 +2383,139 @@ def topo_func(x,y): -def read_netcdf(path, zvar=None, extent='all', coarsen=1, return_topo=True, - return_xarray=False, buffer=0, align=None, verbose=False): +def fetch_remote_topo(name_or_url, crop_extent=None, coarsen=1, buffer=0, + align=None, nc_params={}, verbose=False): + r"""Resolve a remote (or local) netCDF DEM into a `Topography`. + + This is the modern one-call "remote DEM -> Topography" path. It resolves a + nickname or URL and reads it through the `topo_type=4` reader + (`Topography.read`, backed by `netcdf_utils.TopoInspector`), so it inherits + that path's unit handling (elevation must be in meters, or supply + `assume_units` via `nc_params`), datum handling, fill->NaN conversion, CF + coordinate/variable detection, and lazy hyperslab windowing. - """ :Input: - - *path* (str) - Path to the file to read, or url to remote file, - or a key into the topotools.remote_topo_urls dictionary. - - *zvar* (str) - variable to read as Z=elevation. - if None, will try 'Band1', 'z', 'elevation'. - - *extent* - [x1,x2,y1,y2] for desired subset, or 'all' for entire file - - *coarsen* (int) - factor to coarsen by, 1 by default. - - *return_topo* (bool) - if True, return a topotools.Topography object. - default is True - - *return_xarray* (bool) - if True, return an xarray.Dataset object. - default is False - - *buffer* (int): when possible, have at least this many points - outside the filter_region on each side - - *align* (tuple): (xalign,yalign) = desired alignment if coarsening - See the doc string for Topography.crop() + - *name_or_url* (str) - a key into `topotools.remote_topo_urls`, or a URL + (OPeNDAP/THREDDS `dodsC` URLs are read by xarray's netCDF4 backend), or a + path to a local netCDF file. + - *crop_extent* ([x1, x2, y1, y2] or None) - requested crop in domain + coordinates; `None` reads the whole file. Only the requested hyperslab + is read from a remote file. + - *coarsen* (int) - factor to coarsen by (1 = no coarsening). + - *buffer* (int) - when possible, keep at least this many points outside + `crop_extent` on each side. + - *align* ((xalign, yalign) or None) - desired alignment when coarsening; + see `Topography.crop`. + - *nc_params* (dict) - options forwarded to the `topo_type=4` reader, e.g. + `z_var` (elevation variable name) or `assume_units` (unit to assume when + the file has no `units` attribute). See `Topography.read`. + - *verbose* (bool) - if True, print the resolved source. :Output: - - topo and/or xarray_ds depending on what was requested. - (either a single object or a tuple of two objects.) - If `return_xarray == True` then `xarray` is used to read the data, - otherwise `netCDF4` is used directly. + - a `topotools.Topography` object. + + Remote-read failures propagate as `OSError`/`RuntimeError` so callers (and + tests marked `@pytest.mark.remote`) can skip when a server is unavailable. Sample usage: from clawpack.geoclaw import topotools - extent = [-126,-122,46,49] - path = 'etopo1' - topo = topotools.read_netcdf(path, extent=extent, coarsen=2, \ - buffer=1, align=(-126,46), verbose=True) + topo = topotools.fetch_remote_topo('etopo22_30sec', + crop_extent=[-126, -122, 46, 49], + coarsen=2, buffer=1, verbose=True) + topo.write('etopo_sample.tt3', topo_type=3) + """ - # results in topo.x = array([-126.03333333, -126., ...]) + # Resolve a nickname; otherwise treat as a URL or local path. + if name_or_url in remote_topo_urls: + url = remote_topo_urls[name_or_url] + else: + url = name_or_url - # to plot: - topo.plot() + if verbose: + print("Will read netCDF data from \n %s" % url) - # to save topofile for input to GeoClaw: - topo.write('etopo_sample_2min.tt3', topo_type=3, Z_format='%.0f') + # Set the preprocessing attributes *before* reading: Topography.__init__ + # reads immediately when constructed with a path, which would default these + # away, so construct empty and read explicitly. + topo = Topography() + topo.crop_extent = crop_extent + topo.coarsen = coarsen + topo.buffer = buffer + topo.align = align - This should give a 2-minute resolution DEM of the Western Washington coast. - Note that etopo1 Z values are integers (vertical resolution is 1 meter) - and using `Z_format='%.0f'` will save as integers to minimize file size. + try: + topo.read(path=url, topo_type=4, nc_params=nc_params) + except (OSError, RuntimeError): + # Remote/OPeNDAP servers are flaky; let callers/tests decide to skip. + raise + except Exception as e: + raise RuntimeError( + "Failed to read remote topo from %s: %s" % (url, e)) from e - Note that the newer etopo 2022 30 arcsecond DEM can be sampled using - path = 'etopo22_30sec', but this topo is aligned differently with e.g. - x = -126. falling half way between points. Also note that Z values in the - newer dataset are no longer integers. - """ + return topo - from numpy import array - import netCDF4 - if return_xarray: - import xarray - # check if path is a key in the remote_topo_urls dictionary: - if path in remote_topo_urls.keys(): - path = remote_topo_urls[path] +def read_netcdf(path, zvar=None, extent='all', coarsen=1, return_topo=True, + return_xarray=False, buffer=0, align=None, verbose=False): - if verbose: - print("Will read netCDF data from \n %s" % path) + r"""Deprecated: read a netCDF DEM into a Topography and/or xarray.Dataset. - assert (type(coarsen) is int) and (coarsen >= 1), \ - '*** coarsen must be a positive integer' + .. deprecated:: + Use :func:`fetch_remote_topo` (or ``Topography.read(topo_type=4)``) + instead. This is now a thin wrapper over :func:`fetch_remote_topo`; the + standalone ``netCDF4``-based reader it used to contain has been removed + in favor of the modern ``topo_type=4`` read path (unit checking, datum, + fill->NaN, CF coordinate/variable detection, lazy hyperslab windowing). - if return_xarray: - f = xarray.open_dataset(path) - else: - f = netCDF4.Dataset(path, 'r') + The legacy signature is preserved: - if 'lon' in f.variables: - x = f.variables['lon'] - elif 'x' in f.variables: - x = f.variables['x'] - else: - print('*** f.variables = ',f.variables) - raise ValueError("*** Unrecognized x, lon in netCDF file") + - *path* (str) - nickname (key of ``remote_topo_urls``), URL, or local file. + - *zvar* (str) - elevation variable name; mapped to ``nc_params['z_var']``. + - *extent* - ``[x1,x2,y1,y2]`` requested crop, or ``'all'`` for whole file. + - *coarsen* (int) - coarsening factor (1 = none). + - *return_topo* (bool) - if True, include a ``Topography`` in the result. + - *return_xarray* (bool) - if True, include an ``xarray.Dataset``. + - *buffer* (int) - points to keep outside the crop on each side. + - *align* (tuple) - alignment when coarsening; see ``Topography.crop``. - if 'lat' in f.variables: - y = f.variables['lat'] - elif 'y' in f.variables: - y = f.variables['y'] - else: - print('*** f.variables = ',f.variables) - raise ValueError("*** Unrecognized y, lat in netCDF file") - - # for selecting subset based on extent, convert to arrays if netCDF4 used: - #if not return_xarray: - - x = array(x) - y = array(y) - - if zvar is None: - if 'Band1' in f.variables: - zvar = 'Band1' - elif 'z' in f.variables: - zvar = 'z' - elif 'elevation' in f.variables: - zvar = 'elevation' - else: - print('*** f.variables = ',f.variables) - raise ValueError("*** Unrecognized zvar in netCDF file") + :Output: + - a ``Topography``, an ``xarray.Dataset``, or a ``(topo, ds)`` tuple, + depending on ``return_topo`` / ``return_xarray`` (unchanged contract). + """ + import warnings + warnings.warn( + "topotools.read_netcdf is deprecated; use " + "topotools.fetch_remote_topo (or Topography.read(topo_type=4)) instead.", + DeprecationWarning, stacklevel=2) - if extent == 'all': - ilower = 0 - iupper = len(x) - 1 - jlower = 0 - jupper = len(y) - 1 - else: - x1,x2,y1,y2 = extent - # find indices of x,y arrays for points lying within extent: - iindex = numpy.where(numpy.logical_and(x >= x1, x <= x2))[0] - jindex = numpy.where(numpy.logical_and(y >= y1, y <= y2))[0] - ilower = iindex[0] - iupper = iindex[-1] - jlower = jindex[0] - jupper = jindex[-1] - - dx_new = coarsen * (x[1] - x[0]) - dy_new = coarsen * (y[1] - y[0]) - - # shift indices if needed for alignment: - if (coarsen > 1) and (align is not None): - xs = numpy.array([x[ilower + i] for i in range(coarsen)]) - offsets = (xs - align[0]) / dx_new - offsets_frac = offsets - numpy.round(offsets) - ioffset = numpy.argmin(abs(offsets_frac)) - ilower = ilower + ioffset - iupper = iupper - numpy.remainder(iupper-ilower, coarsen) - print(f'+++ shifted ilower by ioffset={ioffset} to {ilower}') + assert (type(coarsen) is int) and (coarsen >= 1), \ + '*** coarsen must be a positive integer' - ys = numpy.array([y[jlower + j] for j in range(coarsen)]) - offsets = (ys - align[1]) / dy_new - offsets_frac = offsets - numpy.round(offsets) - joffset = numpy.argmin(abs(offsets_frac)) - jlower = jlower + joffset - jupper = jupper - numpy.remainder(jupper-jlower, coarsen) - print(f'+++ shifted jlower by joffset={joffset} to {jlower}') + # Map the legacy arguments onto the modern helper. + crop_extent = None if (isinstance(extent, str) and extent == 'all') \ + else extent + nc_params = {} + if zvar is not None: + nc_params['z_var'] = zvar - # buffer, checking limits of arrays: - i1 = numpy.maximum(0, ilower - buffer*coarsen) - j1 = numpy.maximum(0, jlower - buffer*coarsen) - i2 = numpy.minimum(len(x)-1, iupper + buffer*coarsen) + 1 - j2 = numpy.minimum(len(y)-1, jupper + buffer*coarsen) + 1 - - xs = x[i1:i2:coarsen] - ys = y[j1:j2:coarsen] - Zs = f.variables[zvar][j1:j2:coarsen, i1:i2:coarsen] - - Zs = array(Zs) - - if 0: - # debugging checks: - xlower,xupper,ylower,yupper = extent - dx_new = xs[1] - xs[0] - dy_new = ys[1] - ys[0] - xlower_outside = (xlower - xs[0]) / dx_new - ylower_outside = (ylower - ys[0]) / dy_new - xupper_outside = (xs[-1] - xupper) / dx_new - yupper_outside = (ys[-1] - yupper) / dy_new - - print(f'+++ fractions of cells outside should be between' \ - + f' {buffer-1} and {buffer} since buffer={buffer}:') - # note: the statement above is not true if filter_region extends - # to or beyond the edges of the original topo self.extent - print(f'+++ xlower_outside={xlower_outside},' \ - + f' xupper_outside={xupper_outside}') - print(f'+++ ylower_outside={ylower_outside},' \ - + f' yupper_outside={yupper_outside}') - - if align is not None: - xalign = (xs[0] - align[0])/dx_new - yalign = (ys[0] - align[1])/dy_new - print(f'+++ x alignment: {xalign} should be integer') - print(f'+++ y alignment: {yalign} should be integer') + topo = fetch_remote_topo(path, crop_extent=crop_extent, coarsen=coarsen, + buffer=buffer, align=align, nc_params=nc_params, + verbose=verbose) - if verbose: - print('Returning a DEM with shape = %s' \ - % str(Zs.shape)) - print('x ranges from %.5f to %.5f with dx = %.8f' \ - % (xs[0], xs[-1], (xs[1]-xs[0]))) - print('y ranges from %.5f to %.5f with dy = %.8f' \ - % (ys[0], ys[-1], (ys[1]-ys[0]))) - if align is not None: - xalign = (xs[0] - align[0])/dx_new - yalign = (ys[0] - align[1])/dy_new - print(f'aligned in x as requested if {xalign} is an integer') - print(f'aligned in y as requested if {yalign} is an integer') output = None - if return_topo: - topo = Topography() - topo.set_xyZ(xs,ys,Zs) output = topo if return_xarray: - # Create a new xarray.Dataset with this subsampled, coarsened data: - dims = (len(xs),len(ys)) - xarray_ds = xarray.Dataset({'z':(dims,Zs)}, coords={'lon':xs, 'lat':ys}) + import xarray + # Rebuild an xarray.Dataset from the resulting Topography so the legacy + # return contract is unchanged. Z has shape (len(y), len(x)). + xarray_ds = xarray.Dataset({'z': (('lat', 'lon'), topo.Z)}, + coords={'lon': topo.x, 'lat': topo.y}) if output is None: output = xarray_ds else: diff --git a/tests/netcdf/test_base_inspector.py b/tests/netcdf/test_base_inspector.py index d287fc948..93f93f33c 100644 --- a/tests/netcdf/test_base_inspector.py +++ b/tests/netcdf/test_base_inspector.py @@ -5,6 +5,7 @@ fill value resolution, crop bound validation, and laziness guarantee. """ import warnings +from pathlib import Path import numpy as np import pytest @@ -297,3 +298,40 @@ def test_coord_variants_smoke(topo_file_factory, coord_kwargs): assert meta.lon_wrap == expected_conv assert meta.y_increasing == expected_increasing + + +# ============================================================ +# URL vs local path handling (PR #726) +# ============================================================ + +def test_url_passed_through_unmangled(monkeypatch): + """ + A remote OPeNDAP/THREDDS URL must reach xarray as a plain string with its + scheme intact -- not wrapped in pathlib.Path, which would collapse + "https://" to "https:/" and turn it into a cwd-relative local lookup. + """ + import clawpack.geoclaw.netcdf_utils as ncu + + captured = {} + + def fake_open_dataset(path, *args, **kwargs): + captured["path"] = path + # We only need to inspect the argument; stop before a real open. + raise RuntimeError("stop before real open") + + monkeypatch.setattr(ncu.xr, "open_dataset", fake_open_dataset) + + url = ("https://www.ngdc.noaa.gov/thredds/dodsC/regional/" + "crescent_city_13_mhw_2010.nc") + with pytest.raises(RuntimeError): + NetCDFInspector(url) + + assert captured["path"] == url # exact string, "https://" intact + assert not isinstance(captured["path"], Path) + + +def test_local_path_still_becomes_path(topo_file_factory): + """A local filesystem path is still wrapped in pathlib.Path.""" + path = topo_file_factory() + with NetCDFInspector(path) as insp: + assert isinstance(insp.path, Path) diff --git a/tests/test_topotools.py b/tests/test_topotools.py index b42d720a8..71b0db1c2 100644 --- a/tests/test_topotools.py +++ b/tests/test_topotools.py @@ -1,5 +1,6 @@ #!/usr/bin/env python +import os import sys from pathlib import Path from urllib.error import URLError @@ -9,6 +10,7 @@ import clawpack.clawutil.data import clawpack.geoclaw.topotools as topotools +import clawpack.geoclaw.etopotools as etopotools # Local test directory and bundled test data testdir = Path(__file__).parent @@ -317,14 +319,13 @@ def test_get_remote_file_remote(tmp_path): # --- ETOPO1 integration tests and helpers --- -def _read_etopo1_topography(coarsen=10, return_xarray=False): - """Read a small ETOPO1 subset for integration testing.""" +def _read_etopo1_topography(coarsen=10): + """Read a small ETOPO1 subset via the modern fetch_remote_topo path.""" try: - return topotools.read_netcdf( + return topotools.fetch_remote_topo( "etopo1", - extent=etopo1_extent, + crop_extent=etopo1_extent, coarsen=coarsen, - return_xarray=return_xarray, verbose=True, ) except (OSError, RuntimeError): @@ -335,7 +336,7 @@ def _read_etopo1_topography(coarsen=10, return_xarray=False): @pytest.mark.netcdf @pytest.mark.remote def test_etopo1_topography(): - """Integration test for reading a remote ETOPO1 subset via topotools.""" + """Cross-reader equivalence: fetch_remote_topo vs the archived ASCII DEM.""" pytest.importorskip("netCDF4") topo1 = _read_etopo1_topography(coarsen=1) @@ -355,23 +356,157 @@ def test_etopo1_topography(): @pytest.mark.python @pytest.mark.netcdf @pytest.mark.remote -def test_etopo1_xarray(): - """Integration test for the xarray-returning ETOPO1 reader path.""" +def test_etopo1_read_netcdf_shim(): + """The deprecated read_netcdf shim still agrees with archived data and + matches the fetch_remote_topo result, while emitting a DeprecationWarning.""" pytest.importorskip("xarray") - topo10, topo10_xarray = _read_etopo1_topography(coarsen=10, return_xarray=True) + try: + with pytest.warns(DeprecationWarning): + topo10, topo10_xarray = topotools.read_netcdf( + "etopo1", + extent=etopo1_extent, + coarsen=10, + return_xarray=True, + verbose=True, + ) + except (OSError, RuntimeError): + pytest.skip("Reading ETOPO1 failed; check whether the remote server is available.") testdata_path = data_dir / "etopo1_10min.asc" topo10input = topotools.Topography() topo10input.read(testdata_path, topo_type=3) + # Shim result agrees with the archived golden ... assert topo10.Z.shape == topo10input.Z.shape assert topo10_xarray["z"].shape == topo10input.Z.shape assert np.allclose(topo10_xarray["z"], topo10input.Z), ( - "topo10_xarray['z'] does not agree with archived data" + "read_netcdf shim does not agree with archived data" + ) + # ... and with the modern helper it now delegates to. + topo10_helper = _read_etopo1_topography(coarsen=10) + assert np.allclose(topo10.Z, topo10_helper.Z), ( + "read_netcdf shim disagrees with fetch_remote_topo" ) +def _fake_fetch_topo(): + """A small synthetic Topography for offline shim/contract tests.""" + x = np.array([0.0, 1.0, 2.0]) + y = np.array([0.0, 1.0]) + Z = np.arange(6, dtype=float).reshape(2, 3) + topo = topotools.Topography() + topo.set_xyZ(x, y, Z) + return topo + + +@pytest.mark.python +def test_read_netcdf_deprecated_shim(monkeypatch): + """read_netcdf warns and maps its legacy args onto fetch_remote_topo.""" + calls = {} + + def fake_fetch(name_or_url, crop_extent=None, coarsen=1, buffer=0, + align=None, nc_params={}, verbose=False): + calls["name_or_url"] = name_or_url + calls["crop_extent"] = crop_extent + calls["coarsen"] = coarsen + calls["nc_params"] = dict(nc_params) + return _fake_fetch_topo() + + monkeypatch.setattr(topotools, "fetch_remote_topo", fake_fetch) + + # extent='all' maps to crop_extent=None and a DeprecationWarning is emitted. + with pytest.warns(DeprecationWarning): + topo = topotools.read_netcdf("etopo1", extent="all") + assert isinstance(topo, topotools.Topography) + assert calls["name_or_url"] == "etopo1" + assert calls["crop_extent"] is None + assert calls["nc_params"] == {} + + # An explicit extent passes through; zvar maps to nc_params['z_var']. + with pytest.warns(DeprecationWarning): + topotools.read_netcdf("etopo1", extent=[-1, 1, -1, 1], zvar="Band1") + assert calls["crop_extent"] == [-1, 1, -1, 1] + assert calls["nc_params"] == {"z_var": "Band1"} + + +@pytest.mark.python +@pytest.mark.netcdf +def test_read_netcdf_return_xarray_contract(monkeypatch): + """return_xarray still yields a (topo, ds) tuple with matching shapes.""" + pytest.importorskip("xarray") + + monkeypatch.setattr(topotools, "fetch_remote_topo", + lambda *a, **k: _fake_fetch_topo()) + + with pytest.warns(DeprecationWarning): + topo, ds = topotools.read_netcdf("etopo1", return_xarray=True) + assert isinstance(topo, topotools.Topography) + assert ds["z"].shape == topo.Z.shape + assert np.allclose(ds["z"], topo.Z) + + +@pytest.mark.python +def test_etopo1_download_returns_topography(monkeypatch, tmp_path): + """etopo1_download returns a Topography without needing the network.""" + # A minimal NGDC-style AAIGrid payload (keyword-first header, no nodata + # line -- etopo1_download inserts one), as the WCS proxy would return. + aaigrid = "\n".join([ + "ncols 3", + "nrows 2", + "xllcorner -125.0", + "yllcorner 48.0", + "cellsize 0.5", + "-1 -2 -3", + "-4 -5 -6", + "", + ]) + + def fake_get_remote_file(url, output_dir=".", file_name=None, + verbose=True, force=False): + with open(os.path.join(output_dir, file_name), "w") as f: + f.write(aaigrid) + + monkeypatch.setattr(clawpack.clawutil.data, "get_remote_file", + fake_get_remote_file) + + topo = etopotools.etopo1_download((-125, -124), (48, 49), + output_dir=str(tmp_path), verbose=False) + assert isinstance(topo, topotools.Topography) + assert topo.Z.shape == (2, 3) + assert np.all(np.isfinite(topo.extent)) + + +@pytest.mark.python +def test_etopo1_download_return_topo_deprecated(monkeypatch, tmp_path): + """The legacy return_topo keyword is accepted but warns and is ignored.""" + aaigrid = "\n".join([ + "ncols 3", + "nrows 2", + "xllcorner -125.0", + "yllcorner 48.0", + "cellsize 0.5", + "-1 -2 -3", + "-4 -5 -6", + "", + ]) + + def fake_get_remote_file(url, output_dir=".", file_name=None, + verbose=True, force=False): + with open(os.path.join(output_dir, file_name), "w") as f: + f.write(aaigrid) + + monkeypatch.setattr(clawpack.clawutil.data, "get_remote_file", + fake_get_remote_file) + + with pytest.warns(DeprecationWarning): + topo = etopotools.etopo1_download((-125, -124), (48, 49), + output_dir=str(tmp_path), + verbose=False, return_topo=False) + # return_topo=False is ignored: a Topography is still returned. + assert isinstance(topo, topotools.Topography) + + def _import_pyplot(): """Import pyplot using a non-interactive backend for test-safe plotting.""" matplotlib = pytest.importorskip("matplotlib") @@ -633,7 +768,7 @@ def test_unstructured_topo(): pytest.importorskip("scipy") fill_topo, topo = _make_unstructured_topo() - topo.interp_unstructured(fill_topo, extent=[0, 1, 0, 1], delta=(1e-2, 1e-2)) + topo.interp_unstructured(fill_topo, crop_extent=[0, 1, 0, 1], delta=(1e-2, 1e-2)) assert not topo.unstructured assert np.isfinite(topo.Z).all() @@ -646,6 +781,90 @@ def test_unstructured_topo(): assert np.allclose(compare_data.Z, topo.Z) +# --- Region-vocabulary unification: crop_extent + deprecated aliases --- + +@pytest.mark.python +def test_crop_crop_extent_and_filter_region_alias(): + """crop(): crop_extent is canonical; filter_region is a deprecated alias + that warns, gives an identical result, and errors if both are supplied.""" + def topo_bowl(x, y): + return 1000.0 * (x**2 + y**2 - 1.0) + topo = topotools.Topography(topo_func=topo_bowl) + topo.x = np.linspace(-1.0, 3.0, 5) + topo.y = np.linspace(0.0, 3.0, 4) + _ = topo.Z # realize lazy Z + + region = [0, 1, 0, 2] + new = topo.crop(crop_extent=region) + with pytest.warns(DeprecationWarning): + old = topo.crop(filter_region=region) + assert np.allclose(new.Z, old.Z) + # positional first arg still binds to crop_extent + assert np.allclose(new.Z, topo.crop(region).Z) + # supplying both is ambiguous + with pytest.raises(TypeError): + topo.crop(crop_extent=region, filter_region=region) + + +@pytest.mark.python +def test_crop_buffer_is_integer_point_count(): + """buffer is an integer grid-point count: a float is truncated via int() + (not a coordinate distance), and buffer>0 keeps extra points each side.""" + def topo_bowl(x, y): + return 1000.0 * (x**2 + y**2 - 1.0) + topo = topotools.Topography(topo_func=topo_bowl) + topo.x = np.linspace(-2.0, 2.0, 9) # dx = 0.5 + topo.y = np.linspace(-2.0, 2.0, 9) + _ = topo.Z + + region = [-0.6, 0.6, -0.6, 0.6] + base = topo.crop(crop_extent=region, buffer=0) + buffered = topo.crop(crop_extent=region, buffer=1) + # A one-point buffer adds a grid point on each side (not a distance). + assert buffered.x.size == base.x.size + 2 + assert buffered.y.size == base.y.size + 2 + + # A float buffer truncates via int() instead of raising a slice error. + truncated = topo.crop(crop_extent=region, buffer=1.5) + assert truncated.x.size == buffered.x.size + assert np.allclose(truncated.Z, buffered.Z) + + +@pytest.mark.python +def test_read_crop_extent_and_filter_region_alias(): + """read(crop_extent=r) equals setting the attribute then reading; the + deprecated filter_region= alias warns and gives the same result.""" + path = data_dir / "etopo1_10min.asc" + region = [-124.8, -124.0, 48.0, 48.5] + + a = topotools.Topography() + a.read(path, topo_type=3, crop_extent=region) + + b = topotools.Topography() + b.crop_extent = region + b.read(path, topo_type=3) + + assert a.Z.shape == b.Z.shape + assert np.allclose(a.Z, b.Z) + + c = topotools.Topography() + with pytest.warns(DeprecationWarning): + c.read(path, topo_type=3, filter_region=region) + assert np.allclose(a.Z, c.Z) + + +@pytest.mark.python +def test_interp_unstructured_extent_alias(): + """interp_unstructured: crop_extent replaces the deprecated extent= kwarg.""" + pytest.importorskip("scipy") + fill_a, topo_a = _make_unstructured_topo() + fill_b, topo_b = _make_unstructured_topo() + topo_a.interp_unstructured(fill_a, crop_extent=[0, 1, 0, 1], delta=(1e-2, 1e-2)) + with pytest.warns(DeprecationWarning): + topo_b.interp_unstructured(fill_b, extent=[0, 1, 0, 1], delta=(1e-2, 1e-2)) + assert np.allclose(topo_a.Z, topo_b.Z) + + def save_unstructured_test_data(output_dir): """Utility function to save unstructured interpolation test data.""" output_dir = Path(output_dir) diff --git a/tests/test_topotools_preprocessing.py b/tests/test_topotools_preprocessing.py index a347f31e3..4cba32032 100644 --- a/tests/test_topotools_preprocessing.py +++ b/tests/test_topotools_preprocessing.py @@ -430,16 +430,19 @@ def test_preprocessing_double_negate_is_identity(tt2_path): # =========================================================================== -# Group 5 — stride × coarsen interaction (NetCDF type 4) +# Group 5 — coarsen for NetCDF type 4 and the deprecated `stride` alias # =========================================================================== @pytest.mark.netcdf -def test_preprocessing_stride_only(nc_topo_path, tmp_path): +def test_preprocessing_stride_deprecated_maps_to_coarsen(nc_topo_path, tmp_path): + """`stride` is deprecated: it warns and maps onto the scalar `coarsen`, so + the result is identical to reading with `coarsen=2`.""" pytest.importorskip("xarray") path, lon_name, lat_name = nc_topo_path t = Topography() - t.read(path, topo_type=4, stride=[2, 2]) + with pytest.warns(DeprecationWarning): + t.read(path, topo_type=4, stride=[2, 2]) assert t.Z.shape == (5, 5) # Values must match stride-2 subsampling of the analytic Z @@ -466,24 +469,39 @@ def test_preprocessing_coarsen_only_netcdf(nc_topo_path, tmp_path): @pytest.mark.netcdf -def test_preprocessing_stride_and_coarsen_compound(nc_topo_path, tmp_path): - """stride=[2,2] then coarsen=2: effective subsampling is stride-4. +def test_preprocessing_coarsen_param_matches_attribute(nc_topo_path): + """Passing coarsen= to read() is equivalent to setting the attribute.""" + pytest.importorskip("xarray") + path, _, _ = nc_topo_path - The shape is floor(10/4) → check against actual output shape. - """ + t_attr = Topography() + t_attr.coarsen = 2 + t_attr.read(path, topo_type=4) + + t_arg = Topography() + t_arg.read(path, topo_type=4, coarsen=2) + + np.testing.assert_array_equal(t_arg.Z, t_attr.Z) + assert t_arg.coarsen == 2 + + +@pytest.mark.netcdf +def test_preprocessing_stride_conflicts_raise(nc_topo_path): + """Per-axis stride is unsupported; conflicting stride+coarsen is an error.""" pytest.importorskip("xarray") - path, lon_name, lat_name = nc_topo_path + path, _, _ = nc_topo_path + # Unequal per-axis stride cannot be expressed by scalar coarsen. t = Topography() - t.coarsen = 2 - t.read(path, topo_type=4, stride=[2, 2]) + with pytest.warns(DeprecationWarning): + with pytest.raises(ValueError): + t.read(path, topo_type=4, stride=[2, 3]) - # After stride=[2,2], array is 5×5. After coarsen=2, it's 3×3 (floor(5/2)+1). - # The exact shape depends on how crop() slices; just check stride-4 values. - Z_orig = _analytic_Z() - for i in range(t.Z.shape[0]): - for j in range(t.Z.shape[1]): - np.testing.assert_allclose(t.Z[i, j], Z_orig[4 * i, 4 * j], rtol=1e-6) + # stride and a conflicting coarsen at once. + t2 = Topography() + with pytest.warns(DeprecationWarning): + with pytest.raises(ValueError): + t2.read(path, topo_type=4, coarsen=3, stride=[2, 2]) # =========================================================================== @@ -1067,7 +1085,7 @@ def test_crop_pushdown_matches_full_read(nc_topo_path): ref = Topography() ref.read(path, topo_type=4) - ref = ref.crop(filter_region=crop) + ref = ref.crop(crop_extent=crop) t = Topography() t.crop_extent = crop @@ -1087,7 +1105,7 @@ def test_crop_pushdown_with_buffer_matches_full_read(nc_topo_path): ref = Topography() ref.read(path, topo_type=4) - ref = ref.crop(filter_region=crop, buffer=1) + ref = ref.crop(crop_extent=crop, buffer=1) t = Topography() t.crop_extent = crop @@ -1108,7 +1126,7 @@ def test_crop_pushdown_with_coarsen_align_matches_full_read(nc_topo_path): ref = Topography() ref.read(path, topo_type=4) - ref = ref.crop(filter_region=crop, coarsen=2, align=(0.0, 0.0)) + ref = ref.crop(crop_extent=crop, coarsen=2, align=(0.0, 0.0)) t = Topography() t.crop_extent = crop @@ -1123,23 +1141,20 @@ def test_crop_pushdown_with_coarsen_align_matches_full_read(nc_topo_path): @pytest.mark.netcdf @pytest.mark.parametrize("buffer", [0, 1]) -def test_crop_pushdown_with_stride_matches_full_read(nc_topo_path, buffer): - """stride + crop_extent: the strided sub-window stays phase-aligned with - striding the full array (the low index is snapped to a multiple of stride), - so the result matches full-strided-then-cropped. buffer=1 forces an odd - raw low index, exercising the snap.""" +def test_crop_pushdown_with_coarsen_buffer_matches_full_read(nc_topo_path, buffer): + """coarsen + crop_extent + buffer survive the pushdown: the read hyperslab + equals a full read then crop(coarsen, buffer). buffer=1 exercises the + buffer expansion of the read window.""" pytest.importorskip("xarray") path, _, _ = nc_topo_path crop = [3.0, 8.0, 3.0, 8.0] ref = Topography() - ref.read(path, topo_type=4, stride=[2, 2]) - ref = ref.crop(filter_region=crop, buffer=buffer) + ref.read(path, topo_type=4) + ref = ref.crop(crop_extent=crop, coarsen=2, buffer=buffer) t = Topography() - t.crop_extent = crop - t.buffer = buffer - t.read(path, topo_type=4, stride=[2, 2]) + t.read(path, topo_type=4, crop_extent=crop, coarsen=2, buffer=buffer) np.testing.assert_array_equal(t.x, ref.x) np.testing.assert_array_equal(t.y, ref.y) @@ -1162,7 +1177,7 @@ def test_crop_pushdown_respects_storage_order(tmp_path, s2n): ref = Topography() ref.read(str(path), topo_type=4) - ref = ref.crop(filter_region=crop) + ref = ref.crop(crop_extent=crop) t = Topography() t.crop_extent = crop @@ -1201,20 +1216,89 @@ def test_crop_pushdown_reads_bounded_window(nc_topo_path, monkeypatch): path, _, _ = nc_topo_path calls = [] - orig = topotools._netcdf_window_indices + orig = topotools._crop_indices - def _spy(coords, lo, hi, margin, n, stride=1): - result = orig(coords, lo, hi, margin, n, stride) - calls.append((result, n)) + def _spy(x, y, crop_extent, coarsen, buffer, align): + result = orig(x, y, crop_extent, coarsen, buffer, align) + calls.append((result, len(x), len(y))) return result - monkeypatch.setattr(topotools, "_netcdf_window_indices", _spy) + monkeypatch.setattr(topotools, "_crop_indices", _spy) t = Topography() t.crop_extent = [3.0, 6.0, 3.0, 6.0] t.read(path, topo_type=4) assert calls, "crop_extent read did not use the window pushdown" - for (i0, i1), n in calls: - assert 0 <= i0 < i1 <= n - assert (i1 - i0) < n, "window spans the whole axis (no pushdown)" + for result, nx, ny in calls: + assert result is not None, "crop_extent overlaps the file but got None" + il, iu, jl, ju = result + assert 0 <= il < iu <= nx + assert 0 <= jl < ju <= ny + assert (iu - il) < nx, "x window spans the whole axis (no pushdown)" + assert (ju - jl) < ny, "y window spans the whole axis (no pushdown)" + + +# =========================================================================== +# Group 12 — ASCII/NetCDF read equivalence for coarsen+align (rjl report) +# +# The bug: reading the same DEM as NetCDF (type 4) vs ASCII (type 3) with the +# same coarsening produced grids misaligned by fractions of a coarse cell, +# because `stride` coarsened+aligned NetCDF only and used a different alignment +# convention than crop(). These tests pin the unified behavior: read() takes +# `coarsen`/`align` and both file types produce identical, lattice-aligned grids. +# =========================================================================== + +@pytest.mark.netcdf +@pytest.mark.parametrize("align", [None, (0.0, 0.0)]) +def test_ascii_netcdf_coarsen_align_identical(nc_topo_path, tmp_path, align): + """Same data read as NetCDF and as ASCII, with the same coarsen+align, must + yield identical grids (the core of rjl's report).""" + pytest.importorskip("xarray") + path, _, _ = nc_topo_path + crop = [1.0, 8.0, 1.0, 8.0] + coarsen = 2 + + tn = Topography() + tn.read(path, topo_type=4, crop_extent=crop, coarsen=coarsen, align=align) + + # Round-trip the same data through ASCII (type 3) and read it back the same way. + full = Topography() + full.read(path, topo_type=4) + asc = tmp_path / "roundtrip.asc" + full.write(str(asc), topo_type=3, Z_format="%.10f") + + ta = Topography() + ta.read(str(asc), topo_type=3, crop_extent=crop, coarsen=coarsen, align=align) + + np.testing.assert_allclose(ta.x, tn.x) + np.testing.assert_allclose(ta.y, tn.y) + np.testing.assert_allclose(ta.Z, tn.Z) + + +def test_coarsen_align_lattice_invariant(tmp_path): + """With a fixed align, shifting crop_extent by whole native cells keeps the + coarsened grid on the align lattice -- (x0-align)/dx_new stays integer -- + rather than drifting by 1/coarsen of a coarse cell (the reported symptom).""" + n = 30 + x = np.arange(n, dtype=float) + y = np.arange(n, dtype=float) + X, Y = np.meshgrid(x, y) + Z = X + 10.0 * Y + t = Topography() + t.set_xyZ(X, Y, Z) + asc = tmp_path / "invariant.asc" + t.write(str(asc), topo_type=3, Z_format="%.6f") + + coarsen = 3 + align = (0.0, 0.0) + for k in range(coarsen): + crop = [6.0 + k, 20.0 + k, 6.0 + k, 20.0 + k] + tk = Topography() + tk.read(str(asc), topo_type=3, crop_extent=crop, + coarsen=coarsen, align=align) + # dx_new = dx*coarsen = 1*3; alignment => (x0-align)/dx_new integer + xphase = (tk.x[0] - align[0]) / coarsen + yphase = (tk.y[0] - align[1]) / coarsen + assert abs(xphase - round(xphase)) < 1e-9, (k, tk.x[0]) + assert abs(yphase - round(yphase)) < 1e-9, (k, tk.y[0])