verkkoalgoritmi
Ari Jolma
ajolma at water.hut.fi
Fri Nov 19 03:50:25 CST 1999
Tämä ei ole varsinaisesti Perl asia vaan yleisempi mutta silti. Puuhailen
Perl-moduulin
kimpussa, jonka on tarkoitus olla a) yleinen rastereiden (tai mitä termiä
haluaakin käyttää,
kyse on matriiseista) käpistelytyökalu, b) yleinen GIS työkalu ja c)
yleinen hydrologinen
GIS työkalu. Nämä rasterit ovat isoja (esim. 3000x3000) ja siksi algoritmit
on pääasiassa C:llä ja Perlillä on tehty näppärä oliopohjainen ja
operaattorit määrittelevä API.
Mutta siis siihen algoritmiin, joka tuottaa ongelmia koska ne
implementaatiot jotka olen
saanut aikaan ovat hitaita - luokkaa >10 tuntia isolle rasterille ajettuna
eli esim. 10 x nopeutus auttaisi tosi paljon (ajelen näitä omilla
PIII:lla). Mulla on siis FDG (Flow Direction Grid), jonka jokaisessa
pisteessä (i,j, jotka kuten matriiseissa) on kokonaisluku 0..9, 0
tarkoittaa, että solu on "pit" eli kuoppa, 1 tarkoittaa, että vesi virtaa
tästä pohjoiseen eli soluun (i-1,j), 2 tarkoittaa, että vesi virtaa
koilliseen eli soluun (i-1,j+1), jne. Ideana olisi poistaa painanteet eli
alueet, joissa vesi valuu soluun, joka ei ole gridin reunalla. Olen tehnyt
sen, niin, että (i) etsin alueen reunalta alimman solun (ii) muutan polun
siitä pitiin
kulkemaan toisinpäin ja (iii) sen itsensä purkautumaan alueen ulkopuolelle
alimpaan
naapuriinsa.
Olen sillä tavalla hölmö, että minusta näitä on hauska väsätä itse ja tällä
kertaa se voi olla
syynä siihen että tulos on noin hidas. Ehkäpä algoritmikirjallisuudessa
olisi parempi ratkaisu, ja joku voisi viitata minut siihen tai antaa
mielestään nopean koodin niin testaan.
Nopea ratkaisu voisi olla sellainen, joka kävisi koko alueen läpi ja hakisi
reunapisteet
siinä samalla.
Vielä sen verran että FDG lasketaan DEMistä, joka on korkeusmalli mutta
tällä seikalla ei ole tämän algoritmin kanssa tekemistä.
Toivottavasti en pilaa nyt kenenkään viikonloppua ;-)
Ari
tässä mun tämänhetkinen versio:
int fixpit(GRID *fdg, GRID *dem, GRID *mark, int i, int j) {
int i0,j0,i2,j2,ilbp,jlbp,ilbpn,jlbpn,d;
REAL lbp,lbpn;
#ifdef DEBUG_ENABLED
if (debug) fprintf(stderr,"fix pit at %i,%i\n",i,j);
#endif
memset(mark->sdata,0,mark->M*mark->N*sizeof(INTEGER));
markbasin(fdg, mark, i, j);
#ifdef DEBUG_ENABLED
if (debug>1) fprintf(stderr,"basin marked %i,%i\n",i,j);
#endif
/* walk the border, start from top and walk clockwise */
i2 = i;
j2 = j;
while (i2 > 0 AND sgd(mark,i2-1,j2)) {
#ifdef DEBUG_ENABLED
if (debug>3) fprintf(stderr,"up to border, at cell %i,%i\n",i2,j2);
#endif
i2--;
}
/* d is the direction where we go */
d = 3;
#ifdef DEBUG_ENABLED
if (debug>3) fprintf(stderr,"start walking, at cell %i,%i, going
%i\n",i2,j2,d);
#endif
ilbp = i0 = i2;
jlbp = j0 = j2;
lbp = fgd(dem,i2,j2);
do {
/* test clockwise, the first dir to test depends on our current
direction */
/* test direction is d2 */
int d2 = -1;
int done = -1;
#ifdef DEBUG_ENABLED
if (debug>3) fprintf(stderr,"walking, at cell %i,%i, going %i, current
low is %i,%i %f\n",i2,j2,d,ilbp,jlbp,lbp);
#endif
switch (d) {
sagain:
case 1:
case 2:
{
if (i2 > 0 AND j2 > 0 AND sgd(mark,i2-1,j2-1)) {
d2 = 8;
i2--;
j2--;
break;
}
if (i2 > 0 AND sgd(mark,i2-1,j2)) {
d2 = 1;
i2--;
break;
}
}
case 3:
case 4:
{
if (i2 > 0 AND j2 < mark->N-1 AND sgd(mark,i2-1,j2+1)) {
d2 = 2;
i2--;
j2++;
break;
}
if (j2 < mark->N-1 AND sgd(mark,i2,j2+1)) {
d2 = 3;
j2++;
break;
}
}
case 5:
case 6:
{
if (i2 < mark->M-1 AND j2 < mark->N-1 AND sgd(mark,i2+1,j2+1)) {
d2 = 4;
i2++;
j2++;
break;
}
if (i2 < mark->M-1 AND sgd(mark,i2+1,j2)) {
d2 = 5;
i2++;
break;
}
}
case 7:
case 8:
{
if (i2 < mark->M-1 AND j2 > 0 AND sgd(mark,i2+1,j2-1)) {
d2 = 6;
i2++;
j2--;
break;
}
if (j2 > 0 AND sgd(mark,i2,j2-1)) {
d2 = 7;
j2--;
break;
}
}
}
/* if new direction is not found, try once more */
done++;
#ifdef DEBUG_ENABLED
if (debug>3) fprintf(stderr,"do I have the new dir (%i)? %i\n",d2,done);
#endif
if (d2 == -1 AND !done) goto sagain;
d = d2;
if (fgd(dem,i2,j2) < lbp) {
lbp = fgd(dem,i2,j2);
ilbp = i2;
jlbp = j2;
}
} while (!(i2 == i0 AND j2 == j0));
#ifdef DEBUG_ENABLED
if (debug>1) fprintf(stderr,"lowest border point is at %i,%i\n",ilbp,jlbp);
#endif
/*
if the lowest border point is at the border of the whole grid then
that will become a pit, else search for the lowest neighbor not in the
depression
*/
ilbpn = -1;
jlbpn = -1;
if (ilbp > 0 AND ilbp < mark->M-1 AND jlbp > 0 AND jlbp < mark->N-1) {
lbpn = 99999999.9;
i2 = ilbp;
j2 = jlbp;
i2--;
if (i2 >= 0 AND !sgd(mark,i2,j2) AND fgd(dem,i2,j2) < lbpn) {
lbpn = fgd(dem,i2,j2);
ilbpn = i2;
jlbpn = j2;
}
j2++;
if (i2 >= 0 AND j2 < fdg->N AND !sgd(mark,i2,j2) AND fgd(dem,i2,j2) <
lbpn) {
lbpn = fgd(dem,i2,j2);
ilbpn = i2;
jlbpn = j2;
}
i2++;
if (j2 < fdg->N AND !sgd(mark,i2,j2) AND fgd(dem,i2,j2) < lbpn) {
lbpn = fgd(dem,i2,j2);
ilbpn = i2;
jlbpn = j2;
}
i2++;
if (i2 < fdg->M AND j2 < fdg->N AND !sgd(mark,i2,j2) AND
fgd(dem,i2,j2) < lbpn) {
lbpn = fgd(dem,i2,j2);
ilbpn = i2;
jlbpn = j2;
}
j2--;
if (i2 < fdg->M AND !sgd(mark,i2,j2) AND fgd(dem,i2,j2) < lbpn) {
lbpn = fgd(dem,i2,j2);
ilbpn = i2;
jlbpn = j2;
}
j2--;
if (i2 < fdg->M AND j2 >= 0 AND !sgd(mark,i2,j2) AND fgd(dem,i2,j2) <
lbpn) {
lbpn = fgd(dem,i2,j2);
ilbpn = i2;
jlbpn = j2;
}
i2--;
if (j2 >= 0 AND !sgd(mark,i2,j2) AND fgd(dem,i2,j2) < lbpn) {
lbpn = fgd(dem,i2,j2);
ilbpn = i2;
jlbpn = j2;
}
i2--;
if (i2 >= 0 AND j2 >= 0 AND !sgd(mark,i2,j2) AND fgd(dem,i2,j2) < lbpn) {
lbpn = fgd(dem,i2,j2);
ilbpn = i2;
jlbpn = j2;
}
if (ilbpn == -1) {
fprintf(stderr,"fixpit: FAILED: the lowest border point of a
depression at point (%i,%i) does not have a neigbor outside the
depression\n",i,j);
return 0;
}
#ifdef DEBUG_ENABLED
if (debug>1) fprintf(stderr,"escape from the depression:
%i,%i,%i,%i\n",ilbp,jlbp,ilbpn,jlbpn);
#endif
}
if (ilbp != i OR jlbp != j) {
/* invert the path from the lowest border point to the pit point */
int im = ilbp;
int jm = jlbp;
int da = 0;
int db = sgd(fdg,im,jm);
int tmp;
i2 = ilbp;
j2 = jlbp;
do {
movetodir(i2, j2, db);
inversedir(da, tmp);
sgd(fdg,im,jm) = tmp;
im = i2;
jm = j2;
da = db;
db = sgd(fdg,im,jm);
} while (db);
inversedir(da, tmp);
sgd(fdg,im,jm) = tmp;
}
/* a new pit point at the border of the grid or escape from the
depression */
if (ilbpn == -1) sgd(fdg,ilbp,jlbp) = 0;
else {
if (jlbp == jlbpn) {
if (ilbpn-ilbp == 1) sgd(fdg,ilbp,jlbp) = 5;
else if (ilbp-ilbpn == 1) sgd(fdg,ilbp,jlbp) = 1;
} else if (jlbp-jlbpn == 1) {
if (ilbpn == ilbp) sgd(fdg,ilbp,jlbp) = 7;
else if (ilbpn-ilbp == 1) sgd(fdg,ilbp,jlbp) = 6;
else if (ilbp-ilbpn == 1) sgd(fdg,ilbp,jlbp) = 8;
} else if (jlbpn-jlbp == 1) {
if (ilbpn == ilbp) sgd(fdg,ilbp,jlbp) = 3;
else if (ilbpn-ilbp == 1) sgd(fdg,ilbp,jlbp) = 4;
else if (ilbp-ilbpn == 1) sgd(fdg,ilbp,jlbp) = 2;
}
}
#ifdef DEBUG_ENABLED
if (debug>1) fprintf(stderr,"pit fixed, hopefully\n");
#endif
return 1;
}
More information about the Helsinki-pm
mailing list