Fortran for Geospatial work

Help
1

Hello All,

I am new to this forum but not very new to Fortran and did use Fortran a lot during my PhD and post-doc working on legacy codes in hydraulics, hydrology and CFD. However, these days I am with a consulting company working with flood models where things are primarily done in Python and C# and I often get laughed at :slight_smile: . But I still prefer to do things with Fortran as that is the language I am most comfortable in.

flowdir
flowdir5047Ă—4781 492 KB

I am writing here as I am stuck with a problem, which I want to solve in Fortran, and could use some help. I have a DEM (digital elevation model) and flow direction file for a catchment which is 220 km2. The resolution of the file is 5 m and every valid pixel flows into one of the eight neighbouring pixel and finally finds its way out at the catchment at the southern end. Time to one of the eight neighbouring pixel changes every hour but is readily computable, I want to compute the total time to the outlet at the southern end for all the pixels.

As of now in my code, I start with any valid pixel and keep on traversing down the flow direction map and adding the individual time to next pixel for each cell until I reach the outlet. As the time-to-next pixel changes every hour this is the most time consuming part (3-5 mts) of the code which makes things worse if I want to conduct simulation for events spanning couple of days. I can certainly not store the complete drainage path for all the cells because I run out of memory when I am working with over 8 to 10 million grid cells.

If this problem is very trivial in someone’s head please do put in a line or two and I would be extremely grateful.

Kind regards
Sumit

3 Likes
2

I have never worked with digital drainage networks. Isn’t there some readily available method in the literature?

If I understand correctly, you have a 2D array which stores the elevation (let’s assume elevation is > 0, so inactive cells are represented with -1):

type :: elevation_map
   real, allocatable :: elevation(:,:)  ! each element stores the elevation 
                                        ! of a pixel is a square 5 m wide
end type

Next, you have a flow map, whereby each pixel has an associated vector [vx,vy] where vx,vy take on the values -1,0,1, depending into which neighbor pixel they drain. The duration for this drainage is also available. The drainage durations vary in time.

type :: flow_map
   integer, allocatable :: flow_vector(:,:,:)   ! 2*nx*ny
   real, allocatable :: drainage_time(:,:)
end type

Starting from each pixel, you can calculate the total drainage time, by “travelling” down the flow map, and summing the duration, in pseudo-code:

t = 0
loop:
t = t + 1 
read flow_map

do for each pixel (x,y):
  if elevation(x,y) > 0:
    total_time = 0
    p = [x,y]
    do until reach outlet:
       p = p + flow_vector(:,p%x,p%y)
       total_time = total_time + drainage_time(p%x,p%y)
    end do
    ! p is now outlet position
    total_drainage_time(x,y) = total_time
  end if
end do

if (t > tmax) exit

end loop

The total drainage time calculation is repeated in a time loop, with a timestep of 1h for a couple of days or circa O(100) steps? If I understand this correctly, cells might “share” the same drainage path, which could be used to speed up the traversal and updates?

3 Likes
3

Hello Ivan,

First of all thanks a lot for your kind reply, extremely appreciated. I think, I was trying to do something similar, however, the difference in my approach is that I stored all the valid elements in a 1D array thinking it will make the code faster. But when I start with a pixel and follow the flow direction map, I always come out of the loop for mapping row and column to the pertinent 1d element. At any rate, I am extremely thankful for another pair of eyes to take a look at it.

Best wishes,
Sumit

1 Like