просто имитируем физический процесс. Сначала в каждой точке смотрим куда из нее потечет вода - в самую нижнюю из соседних точек. Если все соседние точки не ниже, то ищем ближайшую точку, которая ниже и сливаем все туда. Потом по расставленным градиентам запускаем итерации - везде разливаем одинаковое количество воды и крутим в цикле, пока все не стечет к минимуму (глобальному). Вот собственно исходник:
#include <stdio.h>
#include <math.h>

#define LAT 54.0
FILE *hgt, *res, *pic, *picwater;
int x,y,yi,xi,H,W;
int hmin=255*255;
int hmax=0;
int r,g,b;
//поле water - текущее количество воды в точке,
//поле allwater - суммарное количество воды, прошедшее через точку (если умножить на площадь
//"квадратика", то получится площадь водосбора
struct high {
int h;
long water, allwater;
int nx,ny;
} highes[1203][1203];
unsigned char string[255];
double U;
unsigned char move;
int numsteps;

typedef unsigned long DWORD; // Двойное слово - 32 бита (разряда)
typedef unsigned short WORD; // Слово - 16 бит (разрядов)
typedef signed long LONG;
typedef unsigned short UINT;


// Заголовок файла
typedef struct tagBITMAPFILEHEADER
{
UINT bfType; // 'BM' = 4D42h
DWORD bfSize;
UINT bfReserved1;
UINT bfReserved2;
DWORD bfOffBits; // Смещение к растру
} tBFH;

// Заголовок Bitmap
typedef struct tagBITMAPINFOHEADER
{
DWORD biSize;
LONG biWidth;
LONG biHeight;
WORD biPlanes;
WORD biBitCount;
DWORD biCompression;
DWORD biSizeImage;
LONG biXPelsPerMeter;
LONG biYPelsPerMeter;
DWORD biClrUsed;
DWORD biClrImportant;
} tBIH;

tBFH BFH;
tBIH BIH;
double dx2, dy2;
double dmin;


void main(){

dx2 = 30.93*cos(LAT / 180.0 * M_PI);
dy2 = 30.93;
dx2*=dx2; dy2*=dy2; /* в квадрате... */

hgt=fopen("n54e110.hgt","rb");
res=fopen("res.txt","wt");
pic=fopen("pic.bmp","wb");
picwater=fopen("picwater.bmp","wb");
//инициализируем массив
for (y=1;y<1202;y++)
for (x=1;x<1202;x++){
H=getc(hgt)<<8;
H+=getc(hgt);
if (hmin>H) hmin=H;
if (hmax<H) hmax=H;
highes[x][y].h=H;
highes[x][y].water=1;
highes[x][y].allwater=0;
highes[x][y].nx=0;
highes[x][y].ny=0;
}
for (x=0;x<1203;x++){
highes[x][0].h=32000;
highes[x][0].water=0;
highes[x][0].allwater=0;
highes[x][0].nx=0;
highes[x][0].ny=0;
highes[x][1202].h=32000;
highes[x][1202].water=0;
highes[x][1202].allwater=0;
highes[x][1202].nx=0;
highes[x][1202].ny=0;
highes[0][x].h=32000;
highes[0][x].water=0;
highes[0][x].allwater=0;
highes[0][x].nx=0;
highes[0][x].ny=0;
highes[1202][x].h=32000;
highes[1202][x].water=0;
highes[1202][x].allwater=0;
highes[1202][x].nx=0;
highes[1202][x].ny=0;
}
fprintf(res,"min: %d, max: %d\n",hmin,hmax);

//заполняем заголовки для писания битмапов
BFH.bfType=0x4d42;
BFH.bfSize=sizeof(BFH)+sizeof(BIH)+1200ul*1200ul*3ul;
BFH.bfReserved1=0;
BFH.bfReserved2=0;
BFH.bfOffBits=sizeof(BFH)+sizeof(BIH);
BIH.biSize=sizeof(BIH);
BIH.biWidth=1200;
BIH.biHeight=1200;
BIH.biPlanes=1;
BIH.biBitCount=24;
BIH.biCompression=0;
BIH.biSizeImage=0;
BIH.biXPelsPerMeter=11200ul;
BIH.biYPelsPerMeter=11200ul;
BIH.biClrUsed=0;
BIH.biClrImportant=0;

fwrite(&BFH, sizeof(BFH), 1, pic);
fwrite(&BIH, sizeof(BIH), 1, pic);
fwrite(&BFH, sizeof(BFH), 1, picwater);
fwrite(&BIH, sizeof(BIH), 1, picwater);


Продолжение следует...

//рисуем стандартную картинку (с) slazav
for (y=1200;y>=1;y--){
for (x=1;x<1201;x++){
/* высота в метрах: синий - 500, зеленый - 1500, красный 2500*/
H=highes[x][y].h;
if (H < 0) {r = 0; g= 0; b = 255;}
else if (H < 750) {r = 0; g= (H*256)/750; b = 255;}
else if (H < 1500) {r = 0; g= 255; b = 255-((H-750)*256)/750;}
else if (H < 2250) {r = ((H-1500)*256)/750; g= 255; b = 0;}
else if (H < 3000) {r = 255; g= 255-((H-2250)*256)/750; b = 0;}
else if (H < 10000){r = 255; g= 0; b = 0;}
else {r = 255; g= 255; b = 255;}
/* угол наклона в градусах */

U = sqrt((highes[x][y].h-highes[x+1][y].h)*(highes[x][y].h-highes[x+1][y].h)/dx2
+ (highes[x][y].h-highes[x][y+1].h)*(highes[x][y].h-highes[x][y+1].h)/dy2);
U = atan(U)*180.0/M_PI;

fputc((char)(b*(1-U/90.0)),pic);
fputc((char)(g*(1-U/90.0)),pic);
fputc((char)(r*(1-U/90.0)),pic);
}
}
//считаем куда потечет вода (все тривиально - где ниже - туда и течет)
for(y=1;y<1202;y++){
for(x=1;x<1202;x++){
//всю вот эту идиоццкую конструкцию из ифов и особенно вложенный двойной цикл в конце
//надо заменить на проверку по заранее просчитанному шаблону точек, упорядоченных по расстоянию
if(highes[x-1][y-1].h<highes[x+highes[x][y].nx][y+highes[x][y].ny].h){highes[x][y].nx=-1;highes[x][y].ny=-1;}
if(highes[x][y-1].h<highes[x+highes[x][y].nx][y+highes[x][y].ny].h){highes[x][y].nx=0;highes[x][y].ny=-1;}
if(highes[x+1][y-1].h<highes[x+highes[x][y].nx][y+highes[x][y].ny].h){highes[x][y].nx=1;highes[x][y].ny=-1;}
if(highes[x-1][y].h<highes[x+highes[x][y].nx][y+highes[x][y].ny].h){highes[x][y].nx=-1;highes[x][y].ny=0;}
if(highes[x+1][y-1].h<highes[x+highes[x][y].nx][y+highes[x][y].ny].h){highes[x][y].nx=1;highes[x][y].ny=0;}
if(highes[x+1][y-1].h<highes[x+highes[x][y].nx][y+highes[x][y].ny].h){highes[x][y].nx=1;highes[x][y].ny=-1;}
if(highes[x+1][y].h<highes[x+highes[x][y].nx][y+highes[x][y].ny].h){highes[x][y].nx=1;highes[x][y].ny=0;}
if(highes[x+1][y+1].h<highes[x+highes[x][y].nx][y+highes[x][y].ny].h){highes[x][y].nx=1;highes[x][y].ny=1;}
if(!(highes[x][y].nx||highes[x][y].ny)){
dmin=1.4*1200;
for(yi=1;yi<1202;yi++){
for(xi=1;xi<1202;xi++){
if((highes[xi][yi].h<highes[x][y].h)&&(sqrt((x-xi)*(x-xi)+(y-yi)*(y-yi))<dmin)){
highes[x][y].nx=xi-x;
highes[x][y].ny=yi-y;
dmin=sqrt((x-xi)*(x-xi)+(y-yi)*(y-yi));
}
}
}
}
}
}
//и все потекло... крутим цикл, пока все не стечет
move=1;
numsteps=0;
while(move){
move=0;
numsteps++;
for(y=1;y<1202;y++){
for(x=1;x<1202;x++){
if(highes[x][y].water&&(highes[x][y].nx||highes[x][y].ny)){
move=1;
highes[x+highes[x][y].nx][y+highes[x][y].ny].water+=highes[x][y].water;
highes[x+highes[x][y].nx][y+highes[x][y].ny].allwater+=highes[x][y].water;
highes[x][y].water=0;
}
}
}
}
fprintf(res,"Number of steps: %d\n",numsteps);
//рисуем выходной битмап
for (y=1200;y>=1;y--){
for (x=1;x<1201;x++){
W=highes[x][y].allwater;
if(W>1000){
r=g=b=255;
}
else {
r=g=b=W*256.0/1000.0;
}
if(highes[x][y].water){
g=b=0;
}
/* угол наклона в градусах */
fputc((char)b,picwater);
fputc((char)g,picwater);
fputc((char)r,picwater);
}
}

}