From Newsgroup: comp.lang.fortran
On 2007-04-25 13:51, smurray444 wrote:
Dear all,
I have a table of data values (6509 in each dimension) which
represents
global point values at a resolution of 1.5-degrees... I want to
perform a
linear interpolation to change the resolution of this data to 5-
minutes.
This will require 18 new values to be created above, below and to both
sides
of each existing data value in the current table (except in instances
where
a value is not surrounded by four values, eg. at the edges of the
table,
where only values available should be used).
If anyone has any helpful tips, or could create some Fortran code to
achieve
this, then I'd be very grateful.
Many thanks
smurray444
Time goes, time comes.
Soon enough it will be 20 years since this post was first made.
smurray444, if you're still out there...
Better late than never. Here is the Fortran code you asked for.
The approach is bilinear interpolation. For each point in your fine
grid, you find which four surrounding coarse points enclose it, then
blend their values weighted by distance. At the edges, the clamping
logic degrades gracefully to linear (two points) rather than
extrapolating outward, which is exactly the behaviour you described.
The scale factor works out cleanly: 1.5 degrees is 90 arcminutes, 90 / 5
= 18, so 17 new points are inserted between each pair of originals,
giving an output grid of (6509-1)*18 + 1 = 117145 points per dimension.
Fair warning: at float64 that is roughly 109 GB of memory. You will want
to process it in row-strips and write to a direct-access file rather
than holding it all at once.
Compile with:
gfortran -O2 -o interpolate interpolate.f90
The code is below. Of course, edit to your own needs, whether they
passed or not. It has been tested, every original grid point is
reproduced exactly in the output to within 1e-10.
Hope the project turned out well, whatever it was for.
```f90
program interpolate
implicit none
! --- Parameters ---
integer, parameter :: COARSE_ROWS = 10
integer, parameter :: COARSE_COLS = 10
integer, parameter :: SCALE = 18 ! 1.5deg / (5min=1/12deg) = 18
integer, parameter :: FINE_ROWS = (COARSE_ROWS - 1) * SCALE + 1
integer, parameter :: FINE_COLS = (COARSE_COLS - 1) * SCALE + 1
! --- Arrays ---
real(8) :: coarse(COARSE_ROWS, COARSE_COLS)
real(8) :: fine(FINE_ROWS, FINE_COLS)
! --- Local variables ---
integer :: i, j, r0, c0, r1, c1
real(8) :: dr, dc, f00, f10, f01, f11
real(8) :: row_frac, col_frac
! --- Fill coarse grid with some fake "terrain" data, for your
problem, implement this yourself fully! ---
! We use a simple formula so results are predictable and checkable:
! value = sin(i) * cos(j) * 1000, giving a smooth hilly surface
call fill_coarse(coarse, COARSE_ROWS, COARSE_COLS)
write(*,'(A)') "=== INPUT: Coarse Grid (10x10) ==="
call print_grid(coarse, COARSE_ROWS, COARSE_COLS)
! --- Bilinear interpolation ---
! For each point in the fine grid, find which coarse cell it falls in,
! then blend the four surrounding coarse values by distance.
do i = 1, FINE_ROWS
do j = 1, FINE_COLS
! Map fine index back to fractional coarse index (0-based math)
row_frac = real(i - 1, 8) / real(SCALE, 8) ! e.g. 0.0,
0.0556, ...
col_frac = real(j - 1, 8) / real(SCALE, 8)
! Lower-left coarse cell (1-based indexing for Fortran arrays)
r0 = floor(row_frac) + 1
c0 = floor(col_frac) + 1
! Clamp upper neighbours to array bounds (handles
bottom/right edges)
r1 = min(r0 + 1, COARSE_ROWS)
c1 = min(c0 + 1, COARSE_COLS)
! Fractional offset within the coarse cell (0.0 to 1.0)
dr = row_frac - real(r0 - 1, 8)
dc = col_frac - real(c0 - 1, 8)
! The four surrounding coarse values
f00 = coarse(r0, c0) ! top-left
f10 = coarse(r1, c0) ! bottom-left
f01 = coarse(r0, c1) ! top-right
f11 = coarse(r1, c1) ! bottom-right
! Bilinear blend
fine(i, j) = f00 * (1.0d0 - dr) * (1.0d0 - dc) &
+ f10 * dr * (1.0d0 - dc) &
+ f01 * (1.0d0 - dr) * dc &
+ f11 * dr * dc
end do
end do
write(*,'(A)') ""
write(*,'(A,I0,A,I0,A)') "=== OUTPUT: Fine Grid (", FINE_ROWS, "x", FINE_COLS, ") - top-left 5x5 shown ==="
call print_subgrid(fine, FINE_ROWS, FINE_COLS, 5, 5)
! --- Verify: original coarse points must be exactly preserved ---
write(*,'(A)') ""
write(*,'(A)') "=== VERIFICATION: Checking original points are
preserved ==="
call verify(coarse, fine, COARSE_ROWS, COARSE_COLS, SCALE)
! --- Write full fine grid to file ---
call write_output(fine, FINE_ROWS, FINE_COLS)
write(*,'(A)') ""
write(*,'(A)') "Done! Full output written to: fine_grid.txt"
end program interpolate
! This will fill the coarse grid with a smooth synthetic surface
subroutine fill_coarse(grid, nrows, ncols)
implicit none
integer, intent(in) :: nrows, ncols
real(8), intent(out) :: grid(nrows, ncols)
integer :: i, j
real(8), parameter :: PI = 3.14159265358979d0
do i = 1, nrows
do j = 1, ncols
grid(i, j) = sin(real(i, 8) * PI / real(nrows, 8)) &
* cos(real(j, 8) * PI / real(ncols, 8)) &
* 1000.0d0
end do
end do
end subroutine fill_coarse
! This prints a full grid nicely
subroutine print_grid(grid, nrows, ncols)
implicit none
integer, intent(in) :: nrows, ncols
real(8), intent(in) :: grid(nrows, ncols)
integer :: i, j
do i = 1, nrows
do j = 1, ncols
write(*, '(F8.2, " ")', advance='no') grid(i, j)
end do
write(*,*)
end do
end subroutine print_grid
! We will print just the top-left corner of a large grid
subroutine print_subgrid(grid, nrows, ncols, show_rows, show_cols)
implicit none
integer, intent(in) :: nrows, ncols, show_rows, show_cols
real(8), intent(in) :: grid(nrows, ncols)
integer :: i, j
do i = 1, min(show_rows, nrows)
do j = 1, min(show_cols, ncols)
write(*, '(F8.2, " ")', advance='no') grid(i, j)
end do
write(*,*)
end do
end subroutine print_subgrid
! Here, we check that every original coarse point survived interpolation.
subroutine verify(coarse, fine, nrows, ncols, scale)
implicit none
integer, intent(in) :: nrows, ncols, scale
real(8), intent(in) :: coarse(nrows, ncols)
real(8), intent(in) :: fine((nrows-1)*scale+1, (ncols-1)*scale+1)
integer :: i, j
real(8) :: diff
real(8), parameter :: TOL = 1.0d-10
logical :: ok
ok = .true.
do i = 1, nrows
do j = 1, ncols
diff = abs(fine((i-1)*scale + 1, (j-1)*scale + 1) -
coarse(i, j))
if (diff > TOL) then
write(*,'(A,I0,A,I0,A,E12.4)') &
" MISMATCH at coarse(", i, ",", j, ") diff=", diff
ok = .false.
end if
end do
end do
if (ok) then
write(*,'(A)') " All original grid points exactly preserved."
end if
end subroutine verify
! Here, we write the full fine grid to a file.
subroutine write_output(grid, nrows, ncols)
implicit none
integer, intent(in) :: nrows, ncols
real(8), intent(in) :: grid(nrows, ncols)
integer :: i, j
open(unit=10, file='fine_grid.txt', status='replace')
do i = 1, nrows
do j = 1, ncols
write(10, '(F10.4, " ")', advance='no') grid(i, j)
end do
write(10, *)
end do
close(10)
end subroutine write_output
```
--- Synchronet 3.21d-Linux NewsLink 1.2