1 UVOD
Uporaba geografskih informacijskih sistemov (GIS) je
v vzponu. Kot pri drugih področjih se tudi tukaj
čedalje bolj uveljavljajo odprtokodni paketi. V članku
je za osnovo uporabljen odprtokodni paket GRASS
GIS (Geographic Resources Analysis Support System)
[1]. GRASS GIS je odprtokodno okolje, ki na enoten
način omogoča urejanje geografskih podatkov in njihovo
analiziranje, vsebuje pa tudi orodja za izdelavo želenih
storitev. Slaba lastnost posameznih orodij v GRASS
Prejet 12. december, 2011
Odobren 12. marec, 2012
je njihova relativna počasnost izračunavanja. Počasnost
izračunavanja izvira iz dejstva, da se v okoljih GIS
uporabljajo relativno obsežne rastrske mape, na katerih
je potreben izračun ter sekvenčni način računanja. V
sekvenčnem načinu izračunavanja je čas, potreben za
izračun neke funkcije, sorazmeren produktu izračuna
same funkcije nad točko in števila točk (velikostjo
zemljevida).
Boljša možnost je vzporedno (paralelno)
izračunavanje. Pri vzporednem računanju lahko
čas računanja v dobršni meri skrajšamo na račun
delov programa, ki jih lahko prevedemo v vzporedno
računanje. Pri vzporednem računanju moramo odpraviti
še nekaj drugih težav, na katere pri sekvenčnem
računanju ne naletimo, kot na primer sinhronizacija
vzporednih opravil. Ena najbolj popularna rešitev v
sedanjem času za vzporedno računanje so grafične
enote z arhitekturo CUDA (Compute Unified Device
Architecture). Nova arhitektura omogoča vzporedno
procesiranje (lahko tudi negrafičnih matematičnih
problemov) na sami grafični kartici ([2]).
V članku se osredinimo na modul r.los za izračun
optične vidljivosti (LOS - Line of Sight), ki je že
implementiran v odprtokodnem okolju GRASS GIS. Se-
kvenčni modul r.los izračuna polje vidnosti opazovalca.
Glavni vhodni podatek je digitalni zemljevid višine
terena (DEM - digital elevation model). Modulu moramo
podati tudi točko interesa oz. opazovanja (koordinato
x,y), višino opazovanja nad terenom in maksimalno
razdaljo opazovanja oz. računanja (max dist). Sekvenčni
modul r.los je sorazmerno počasen za računanje. Zaradi
20 OSTERMAN
dolgega računanja postane tako rekoč neuporaben nad
maksimalno razdaljo opazovanja več kot 40 km pri
digitalnem zemljevidu ločljivosti 100 x 100 m, saj je
čas računanja več kot 30 s. V tem članku predstavimo
nov modul za računanje vidnosti (LOS - Line of Sight).
Modul pridobi vstavek cuda in se imenuje r.cuda.los.
Uporabljeno je vzporedno računanje na grafični kartici
NVIDIA. Izkaže se, da je vzporedni modul r.cuda.los
sposoben izračunati vidnosti za faktor od 10 do 1000
hitreje od sekvenčnega modula r.los. Glavni časovni
faktor pri uporabi modula ni več računanje, temveč
branje in zapisovanje digitalnega zemljevida iz trdega
diska v sistemski spomin gostitelja ter naprej na globalni
spomin GPE. Velja tudi, da moramo izračunane podatke
zopet prenesti iz globalnega spomina GPE (grafične pro-
cesne enote) v sistemski spomin, iz sistemskega spomina
pa jih je treba zapisati na trdi disk.
Zaradi morebitne uporabe modula r.cuda.los v ra-
dijsko planerske namene (hitra ocenitev, kaj antena
’vidi’) so dodani še parametri za azimut, obseg azimuta
(horizontalni kot sevanja ’antene’), nagib (angl. tilt) in
obseg nagiba (vertikalni kot sevanja ’antene’). Modul
se uspešno uporablja pri začetnem planiranju radijskega
sektorja v mobilni tehnologiji.
Članek je organiziran takole: v poglavju 2 je najprej
na kratko opisana priprava digitalnega zemljevida višine
terena. Poglavje 3 govori o pripravi globalnega spo-
mina. Poglavje 4 razdela geometrijo posamezne rezine
(niti) pri računanju vidnosti. Srce programa, ščepec je
opisan v poglavju 5. Tu je podana tudi okrnjena koda
ščepca. V poglavju 6 so podani rezultati izračuna in pri-
merjava med sekvenčnim modulom r.los ter modulom
r.cuda.los. V poglavju 7 podam razpravo, prednosti in
slabosti nove implementacije. Opisane so tudi mogoče
izboljšave novega modula.
2 PRIPRAVA DIGITALNEGA ZEMLJEVIDA
Digitalni zemljevid višine terena (angl. digital elevation
map - DEM) je shranjen na disku v rastrski obliki. Vse-
bino digitalne mape je treba prebrati v spomin gostitelja
ter nato prenesti na napravo GPE. Na napravi GPE se
izvede vzporedno računanje nad podatki. Po končanem
računanju je treba rezultat prenesti iz naprave GPE v
spomin gostitelja ter nato zapisati rezultat na trdi disk.
Branje in pisanje na trdi disk je izvedeno s klasičnima
C-funkcijama fread() in fwrite(). Prenos podat-
kov z gostitelja na napravo GPE in nazaj je izve-
deno s funkcijama cudaMemcpyHostToDevice()
in cudaMemcpyDeviceToHost().
Izkaže se, da večino računalniškega časa (več kot 90
%) porabijo zgoraj omenjene funkcije. Te funkcije so
sekvenčne in opravljajo prenašanja podatkov z enega
medija na drugega. Vzporedno računanje na podatkih
pa se izvede na napravi GPE.
3 PRIPRAVA SPOMINA
Vhodni digitalni zemljevid ima dimenzijo R ∗ C, pri
čemer je R število vrstic, C pa je število stolpcev. Vsaka
točka zemljevida vsebuje celoštevilčno vrednost (v dveh
ali štirih zlogih), ki pomeni nadmorsko višino. Videz
vhodnega spomina je prikazana na sliki 1.
(xobs,yobs,hobs+hpoi)
TOČKA INTERESA
Točka v vhodnem spominu (int ali short)
vhodni spomin 0
vhodni spomin 1
vhodni spomin 2
vhodni spomin 3
Rezina od točke interesa do robne točke
Točka interesa, lokacija opazovalca
Slika 1: Vhodni spomin
Če je digitalni zemljevid večji od približno polovice
razpoložljivega spomina na grafični kartici, moramo
digitalni zemljevid razdeliti na več pasov. Računanje na
grafični kartici izvedemo najprej s tistim pasom, kjer je
točka, iz katere nas zanima izračun vidnosti. Vrstni red
računanja je v našem primeru: vhodni spomin 2; vhodni
spomin 1; vhodni spomin 0; vhodni spomin 3. Omenjeni
spomin je na kartici tako imenovani globalni spomin.
To pomeni, da lahko do tega spomina katerakoli nit iz
kateregakoli bloka enakovredno dostopa. Prav zato ni
pomembno, kako so niti organizirane v bloke.
Za izhodni digitalni zemljevid (kjer bo shranjen rezul-
tat) predvidimo enako velikost R ∗C. Pred tem ga mo-
ramo zapolniti z vrednostmi null, kar se najhitreje lahko
naredi z vzporednim računanjem. Poleg tega pa moramo
rezervirati še štiri double vektorje {Vn, Vs, Vw, Ve}.
Vw, Ve imata R elementov, Vn, Vs imata C elementov.
V te štiri vektorje shranjujemo začasne vertikalne kote,
enote pa so v radianih. Pred začetkom računanja vse
elemente vseh štirih vektorjev napolnimo z vrednostjo
−π/2. Ta vrednost pomeni pogled v tla s strani opazo-
valca. Pogled naravnost proti horizontu ima vrednost 0,
pogled navpično navzgor pa vrednost π/2.
4 RAČUNANJE VIDNOSTI PO POSAMEZNIH
REZINAH
Vidnost računamo po posameznih rezinah, in sicer vedno
začnemo pri točki, ki nas zanima (Point of Interest), in
se po vnaprej izračunanem koraku premikamo do točke
IMPLEMENTACIJA MODULA R.CUDA.LOS V GRASS GIS 21
(xobs,yobs,hobs+hpoi)
TOČKA INTERESA
Točka v izhodnem spominu (int ali short)
Trenutni naklon v robnem vektorju (double)
izhodni spomin 0
izhodni spomin 1
izhodni spomin 2
izhodni spomin 3
Rezina od točke interesa do robne točke
Vn
Vw
Vs
Ve
Točka interesa, lokacija opazovalca
Slika 2: Izhodni spomin in naklonski vektorji
na robu digitalnega zemljevida (glej sliko 2, rdeče črte
od POI proti robu zemljevida). Za izračun vsake rezine
je določena svoja nit (thread) vzporednega procesorja. V
posamezni rezini za vsako točko izračunamo razdaljo d
od začetne točke (enačba (1)). (xobs, yobs) je koordinata
opazovalca, (xd, yd) je koordinata trenutne točke, ki jo
računamo. d je razdalja med tema dvema točkama.
d =
√
(xobs − xd)2 + (yobs − yd)2 (1)
Za vsako točko rezine izračunamo tudi vertikalni kot
α (enačba (2)), pod katerim bi bila (je) vidna točka (slika
3).
TOČKA
INTERESA
h
p
o
i
nagibni kot
vidno območje
α
razdalja d
(xd,yd,hd - hc)
točka interesa
Slika 3: Rezina
α je kot (angl. tilt), pod katerim opazovalec vidi
točko (xd, yd). hd je nadmorska višina točke (xd, yd).
Nadmorska višina, kjer stoji opazovalec, je označena s
hobs. Višina opazovalca nad tlemi pa je hpoi.
α = arctan
(
(hd − hc)− (hobs + hpoi)
d
)
(2)
Enačba (2) ima korekcijski faktor višine hc. Ta na-
stopa zaradi ukrivljenosti Zemlje (slika 4), izračunamo
ga s pomočjo enačbe (3). V tej enačbi vzamemo za
premer Zemlje povprečno vrednost 6370.997 km in ne
upoštevamo nadmorske višine opazovalca (če vnesemo
povprečno oddaljenost od središče Zemlje za želeno
območje, nastane napaka na četrti decimalki).
(hc + rearth)
2 = d2 + r2earth
hc =
√
d2 + r2earth − rearth
(3)
V element vektorja V [tid], ki pripada končni točki
rezine, vpišemo novo vrednost vertikalnega kota α ta-
krat, ko je nova vrednost večja od že zapisane vrednosti.
Hkrati pa dotično točko v izhodnem digitalnem zemlje-
vidu označimo kot vidno (vpišemo vrednost vertikalnega
kota, pod katerim jo vidimo). Če pa je vertikalni kot α
manjši od vpisane vrednosti v vektorju, ne vpišemo v
dotično točko nič (pred tem je že vpisana vrednost null).
Tako imamo v vidnih točkah vpisane vrednosti kota. V
točkah, ki niso vidne, je vpisan null element. Če izhodni
digitalni zemljevid zapisujemo v celoštevilčnem for-
matu, izračunu prištejemo vrednost π/2 in pomnožimo
z 105.
r
e
a
rth
r
e
a
rth
d
h
c
Slika 4: Popravek višine zaradi ukrivljenosti Zemlje
5 REALIZACIJA PROGRAMA, ŠČEPEC
Poenostavljena izvorna koda ščepca (angl.kernel,
osnovni gradnik vzporednega programa, ki se izvaja na
arhitekturi CUDA [2] ) za izračun vidnosti, je izpisana
spodaj.
Na začetku vsaka nit dobi svoj indeks (tid). Organiza-
cija niti po blokih v tem primeru ni pomembna, ker niti
med seboj neposredno ne sodelujejo. Hkrati pa se za
shranjevanje digitalnega zemljevida uporablja globalni
pomnilnik, do katerega imajo vse niti enakovreden sku-
pen dostop. Slabost globalnega pomnilnika je relativna
počasnost dostopa z latenco nekaj 10 urnih period [2].
Vendar tu ni drugih možnosti zaradi velikosti digitalnih
zemljevidov.
Vse vhodne podatke, ki se med delovanjem ščepca
ne spreminjajo (maksimalna razdalja, omejitev azimuta,
22 OSTERMAN
omejitev naklona ...), vpišemo v tako imenovani kon-
stantni spomin (__constant__).
Izkoriščen je indeks threadIdx.y, ki pomeni tako-
imenovane štiri podniti. Namenjen je za nastavljanje in
izvrševanje računanja režnjev na severno, južno, zaho-
dno in vzhodno stran digitalnega zemljevida. Na začetku
izračunamo korak pomika od točke POI do robne točke
zemljevida (xstep in ystep).
V kodi ni podan izračun predznaka koraka pomika,
odvisen pa je od tega, v katero smer od točke POI se
pri izračunu gibljemo.
Glavna zanka dejansko pelje spremenljivki x in y od
točke POI proti robni točki v smeri x s korakom xstep
in v smeri y s korakom ystep. Znotraj glavne zanke
se izračuna razdalja d, kot α in popravek višine hc
za vsako točko v režnju. V originalni kodi se izračuna
še azimut (kot proti severu), na podlagi katerega lahko
omejimo računanje (če je tako nastavljeno z vhodnimi
parametri). Kot α nazadnje primerjamo z vrednostjo
vpisano v elementu vektorja V [tid]. Če je kot α večji od
kota vpisanega v elementu vektorja V [tid], je trenutna
točka računanja (x, y) vidna in zato vpišemo kot α v
izhodno datoteko.
Število niti v ščepcu je odvisno od velikosti digital-
nega zemljevida. Število niti je N = 2 ∗ C + 2 ∗R, pri
čemer je C število stolpcev in R število vrstic v vho-
dnem digitalnem zemljevidu. Kako so niti organizirane
v bloke, v tem konkretnem primeru ni pomembno, ker
vse niti uporabljajo globalni spomin na GPE.
t i d = t h r e a d I d x . x + b l o c k I d x . x * blockDim . x ;
x s t e p = 1 . 0 ;
y s t e p = 1 . 0 ;
sw i t ch ( t h r e a d I d x . y )
{
case 0 :
i f ( t i d>c o l s ) re turn ;
x end= t i d ;
y end =0;
V=V n ;
x s t e p = f a b s ( ( x end − x obs ) / ( y end − y obs ) ) ;
break ;
case 1 :
i f ( t i d>c o l s ) re turn ;
x end= t i d ;
y end=rows ;
V=V s ;
x s t e p = f a b s ( ( x end − x obs ) / ( y end − y obs ) ) ;
break ;
case 2 :
i f ( t i d>rows ) re turn ;
x end =0;
y end= t i d ;
V=V w;
y s t e p = f a b s ( ( y end − y obs ) / ( x end − x obs ) ) ;
break ;
case 3 :
i f ( t i d>rows ) re turn ;
x end= c o l s ;
y end= t i d ;
V=V e ;
y s t e p = f a b s ( ( y end − y obs ) / ( x end − x obs ) ) ;
break ;
}
t i l t =−PI / 2 . 0 ;
/ / l oop from POI t o edge p o i n t
f o r ( x=x obs , y=y obs ; ( x<=c o l s ) and ( y<=rows ) ;
x+= x s t ep , y+= y s t e p )
{
d= s q r t ( ( x−x obs ) * ( x−x obs )
+( y−y obs ) * ( y−y obs ) ) ;
h c= s q r t ( d* r e s o l *d* r e s o l +6370997 .0*6370997 .0 )
−6370997.0;
h= r e a d i n p u t m a p ( x , y ) ;
h=h−h c ;
a l p h a = a t a n ( ( h−h obs ) / d ) ;
i f ( a l p h a > (V[ t i d ]−0.00005) )
{
w r i t e o u t p u t m a p ( x , y , a l p h a ) ;
i f ( a l p h a > V[ t i d ] ) V[ t i d ]= a l p h a ;
}
/ / POI red c o l o r
i f ( d<2.5)
{
w r i t e o u t p u t m a p ( x , y , PI / 2 ) ;
}
}
s y n c t h r e a d s ( ) ;
6 PRIMERJAVA R.LOS IN R.CUDA.LOS
Modul r.cuda.los poženemo iz terminalne vrstice po-
dobno kot modul r.los. Oba modula sta testirana
na osebnem PC računalniku, na katerem teče Linux.
Grafična kartica je NVIDIA GeForce GTX 560Ti, CPU
je Intel Core(TM) Duo CPU E8200 @ 2.66GHZ, veli-
kost spomina je 2.0 GiB.
Grafična kartica NVIDIA GeForce GTX 560Ti ima
naslednje za nas pomembne lastnosti:
- globalni spomin: 1024 MB
- število jeder: 336
- največje število niti na blok: 1024
- verzija: 2.1
Podrobnejši podatki o omenjeni grafični kartici se
nahajajo na spletni strani [8].
Na sliki št. 5 je prikazan primer izračuna z mo-
dulom r.los na digitalnem zemljevidu z ločljivostjo
25mx 25m. Maksimalna razdalja izračuna je 10 km,
višina opazovalca pa 20 m. Isti izračun je narejen z
novim modulom r.cuda.los (slika št. 6). Obe sliki se po
obliki ne razlikujeta. Namerno sta obarvani z različnimi
barvami, kar je narejeno z modulom r.color.
Na sliki št. 6 tople barve (rumena, oranžna, rdeča)
pomenijo pogled opazovalca nad črto navideznega hori-
zonta (negativni naklon, angl.tilt), hladne barve (svetlo
modra, temno modra) pa pogled opazovalca pod črto
navideznega horizonta (pozitivni naklon, angl.tilt).
V tabeli 1 je prikazana primerjava časa
izvrševanja obeh modulov. Prikazana je kombinacija
s tremi različnimi vhodnimi zemljevidi s tremi
različnimi ločljivostmi (100mx 100m, 25mx 25m
in 12, 5mx 12, 5m). Za vsako ločljivost je vzeta
kombinacija štirih razdalj, do katere računamo
(parameter max dist), in sicer 5 km, 10 km, 20 km
IMPLEMENTACIJA MODULA R.CUDA.LOS V GRASS GIS 23
in 50 km. Čas računanja za posamezne kombinacije je
podan v sekundah.
Modul r.los za računanje na veliko točkah (daljše raz-
dalje oz. srednje razdalje z bolj podrobno ločljivostjo) ni
več primeren, ker je čas računanja praktično neuporaben
(na tabeli označeno s pomišljajem).
Z modulom r.cuda.los lahko računamo vidnost čez
območje celotne Slovenije na zemljevidu velikosti
28161x 17921 točk z ločljivostjo 12, 5mx 12, 5m, čas
računanja je okoli 18 s.
Tabela 1: Primerjava časa računanja modulov r.los in r.cuda.los
Map [m x m] max dist r.los r.cuda.los
100m x 100m
5 km 0.1 s 0.05 s
10 km 0.2 s 0.06 s
20 km 2.2 s 0.09 s
50 km 44 s 0.15 s
25m x 25m
5 km 2.4 s 0.3 s
10 km 30 s 0.3 s
20 km 511 s 0.6 s
50 km - 1.3 s
12.5m x 12.5m
5 km 32 s 0.7 s
10 km 516 s 1.2 s
20 km - 3.1 s
50 km - 6.8 s
Slika 5: Izračun vidnosti z modulom r.los
7 RAZPRAVA IN SKLEP
Kdorkoli se je že kdaj srečal z okolji GIS, ve, da je
priprava geografskih podatkov natančno, vestno, pred-
vsem pa časovno zamudno delo. Vsaka premišljena
poteza pri pripravi podatkov se pozneje obrestuje s
hitrejšo in lažjo obdelavo podatkov. Vendar pa je zelo
pomemben dejavnik tudi samo računanje nad podatki.
Slika 6: Izračun vidnosti z modulom r.cuda.los
Ključno vlogo igra pri tem ustrezna strojna oprema, ki
pa mora biti podprta z optimalno programsko opremo.
Obdelava geografskih podatkov kar kliče po vzpore-
dnem računanju. Če primerjamo sekvenčno in vzpore-
dno obdelavo podatkov GIS, lahko dobimo pri vzporedni
obdelavi za en, dva ali celo tri velikostne razrede hitrejšo
izvedbo!
Članek [2] govori na splošno o arhitekturi CUDA.
Poleg strokovne vrednosti ima ta članek tudi velik
pomen zaradi slovenskega izrazoslovja v tej novi veji
tehnike.
V okolju GRASS-GIS je bilo že nekaj raziskav o
vzporednem računanju. Ena izmed njih je [3], kjer so
vzporedno računanje reševali z grozdom strežnikov, kjer
vsak strežnik naredi eno nalogo.
V članku [4] je opisan splošni način programiranja
GIS v CUDA s primerom na filtru za povprečenje
(angl.meanfilter).
V članku [5] je opisana implementacija transformacij
projekcij digitalnih zemljevidov s pomočjo CUDA.
Članek [6] sicer ne opisuje vzporednega računanja,
vendar se skupina že ukvarja z implementacijo svojih
modulov na vzporedni način računanja.
Pričujoči članek ne obravnava splošnega problema
migracije iz sekvenčnega na vzporedni način računanja,
temveč se omeji na konkretno nalogo: izboljšanje ob-
stoječega modula r.los.
Verjetno bo splošen preskok na vzporedni način
računanja nastal s časom, ko bo več konkretnih modulov
’migriralo’ na vzporedni način računanja.
Za modul r.los je na spletu prosto dostopna izvorna
koda. Kljub temu je modul r.cuda.los popolnoma na
novo narejen. Oba modula izračunata enak rezultat ter
oba uporabljata datotečni način zapisovanja podatkov
GRASS. Klic programa je za oba modula precej podo-
ben. Velika večina argumentov je enakih, pri r.cuda.los
24 OSTERMAN
je nekaj argumentov dodanih za nastavitev dodatne
funkcionalnosti.
Novi modul r.cuda.los je pisan za CPU in GPE.
Prevajalnik za to je nvcc (namesto gcc za GRASS).
GRASS ima izredno učinkovito narejeno ogrodje za
izdelavo modulov. To ogrodje vsebuje tudi GUI za zagon
samega modula s pripadajočim ogrodjem za pomoč
pri uporabi modula. Za zdaj novi modul r.cuda.los ni
združen s tem ogrodjem, zato ga lahko poženemo samo
prek komandne vrstice.
Pri novem modulu se omejimo na dva načina zapisa
vhodnega in izhodnega zemljevida. Modul lahko prebere
nestisnjeni način zapisa v obliki int ali short. Rezultat
(izhodni zemljevid) dobimo v enaki obliki, kot so vho-
dni podatki (vhodni zemljevid). Tu se pokaže slabost
trenutne implementacije, saj v izhodno datoteko vpisu-
jemo samo celoštevilčne vrednosti. Če imamo majhne
vrednosti (kot v našem primeru, kjer zapisujemo kota α
z vrednostmi −π/2 do π/2), mora biti vpisana vrednost
skalirana (pomnožena npr. z 105). Zapis v obliki double
še čaka na implementacijo.
Če je vhodni zemljevid zapisan v stisnjeni obliki, ga
moramo pred tem razširiti z modulom r.compress.
Iz tabele 1 vidimo, kakšna prednost nam prinese
vzporedno računanje. Tabela prikazuje čas izvajanja
modula r.los in modula r.cuda.los glede na vhodne
podatke. Modul r.los postane zelo počasen pri velikih
vrednostih max dist. Vidimo, da postane zelo počasen
pri vrednostih večjih od 50 km za digitalne zemlje-
vide z ločljivostjo 100mx 100m, ter pri vrednostih
max dist večjih od 5 km za zemljevide z ločljivostjo
12.5mx 12.5m.
Modul r.cuda.los ima še kar nekaj možnosti za
izboljšave. Predvsem je treba izboljšati časovno potratno
branje digitalne mape z diska in nazaj. Ena izmed
možnosti je zapis digitalne mape v stisnjeni obliki. Za
enoto GPE bi bilo treba v tem primeru napisati ščepec
za hitro iskanje vrednosti točk iz stisnjene oblike zapisa.
Vsekakor je vzporedno računanje prihodnost za sis-
teme GIS. Izvajanje operacij bo pohitril vsaj za veliko-
stni razred.