12 Geovizualizacija u R jeziku

U ovom poglavlju ukratko će biti prikazani primeri geovizualizacije u R jeziku (R Core Team 2018). Prvi deo će u najkraćem prestaviti R jezik i okruženje, sa osnovnim strukturama podataka, potom će biti predstavljene prostorne strukture podataka u R jeziku i na kraju će biti opisani primeri vizualizacije, kako u okviru R okruženja, tako i na virtuelnim globusima i Web kartama. Svi primeri su dati sa kodom i podacima tako da ih čitalac može sam reprodukovati ili prilagoditi svojim potrebama.

12.1 R jezik

R je nastao kao sistem za statističke proračune i grafike, trenutno je jedan od najmoćnijih alata za analizu podataka. R je interpreter koji je besplatan i softver je otvorenog koda. Pored ostalog, R je programski jezik sa velikim potencijalom za kreiranje grafika i vizualizacije, takođe, R nudi interfejs ka drugim programskim jezicima i GIS aplikacijama. R je vrlo popularan za analizu podataka u mnogim oblastima i to u: statistici, geoinformatici, geografiji, poljoprivredi, ekologiji, bionformatici i mnogim drugim disciplinama.

U ovom delu neće biti detaljno prezentovan uvod u R jezik se preporučuje neki od sledećih slobodno dostupnih materijala, posebno ako čitalac nije upoznat sa niti jednim programskim jezikom:

http://cran.r-project.org/doc/contrib/Paradis-rdebuts_en.pdf
http://www.cyclismo.org/tutorial/R/
https://www.stat.auckland.ac.nz/~ihaka/120/notes.html
http://www.springer.com/us/book/9780387938363
http://www.nostarch.com/artofr.htm
http://adv-r.had.co.nz/ (za napredne korisnike)
https://geocompr.robinlovelace.net/ (za napredne korisnike)

12.1.1 Instalacija R-a

Osnovna instalacija R interpretera sadrži set paketa koje omogućavaju linearnu algebru, deskriptivnu statistiku, kreiranje grafika i drugo. Instalacija zavisi od operativnog sistema, a može se preuzeti sa stranice: https://cran.r-project.org/.

Paketi su skupovi funkcija, dokumentacionih fajlova i podataka uvezanih zajedno. Pakete prave R korisnici i stručnjaci iz razvojnog tima. Trenutno ima preko 12000 paketa. Na primer, rgdal paket je set klasa i metoda za čitanje i pisanje prostornih podataka koristeći GDAL funkcionalnosti.

Korisnici uglavnom instaliraju set paketa koji oni koriste za specifične potrebe. Na primer, set paketa koji se odnose na prostorne podatke je dat na url adresi: https://cran.r-project.org/web/views/Spatial.html, na ovoj stranici se može naći spisak paketa sa kratkim opisom, svi ovi paketi su pridruženi grupi - view Spatial.

Ako korisnik želi da instalira sve pakete namenjene radu sa prostornim podacima (koji su grupi - view Spatial) neophodno je da zada sledeće komande.

install.packages("ctv")
library("ctv")
install.views("Spatial")

Prva komanda je za instalaciju paketa koji se zove ctv (Cran Task View), potom se učita paket ctv u radno okruženje koristeći naredbu library. Na kraju koristeći komandu install.views pokreće se instalacija svih paketa koji su klasifikovani u ovu grupu.

Najpopularnije okruženje za R je RStudio (Slika 12.1), koji nudi integrisano okruženje za R interpreter, editor za pisanje koda, prozor za pregled podataka, grafika, pomoć za funkcije i još mnogo toga. Instalacija je dostupna na url adresi: https://www.rstudio.com/products/rstudio/download/.

Slika 12.1: Izgled RStudio okruženja

12.1.2 Korišćenje ugrađenjih funkcija i operatora u R okruženju

Ako korisnik pokrene R konzolu direktno (Slika 12.2), pokrenuće se interpreter koji očekuje komande, koje se direktno izvršavaju. Najprostiji pristup je koristiti R kao običan kalkulator.

Slika 12.2: R interpreter; R konzola

U Tabeli 12.1 dati su najvažniji matematički i logički operatori u R jeziku.

Tabela 12.1: R operatori sa opisom
Operator	Opis	Operator	Opis
<-	Dodeljivanje
+	Suma	&	i
-	Razlika	<	Manje
*	Množenje	>	Veće
/	Deljenje	<=	Manje ili =
^	Eksponent	>=	Veće ili =
%%	Moduo	!	Ne
*%%**	Matrično množenje	!=	Različito
%/%	Celobrojno deljenje	==	Isto

Kad se dodeljuje promenljiva u R-u, tada se dodeljuje simbol (uobičajeno skraćenica ili ime ili bilo koji niz karaktera) za vrednost koja postoji u R okruženju. Promenljive u R-u su osetljive na velika i mala slova (eng. CASE SENSITIVE), to znači da promenljiva x nije ista kao promenljiva X, ili xX nije isto što xx. U R se dodeljivanje može vršiti koristeći operator <- ili =.

x <- c(1,2,3,4,8,11,18)
mean(x)
## [1] 6.714286
? mean
help(mean)

U ovom primeru, promenljivoj x je dodeljen je vektor od 7 elemenata koristeći komandu c (eng. combine). Potom je pozvana ugrađena funkcija mean za proračun srednje vrednosti, i dat je primer kako da se dođe do opisa funkcije i njenih argumenata koristeći ugrađene opise funkcija (help za funkcije). Slično, možemo promenljivoj x dodeliti drugi vektor, pa potom napraviti histogram podataka, Slika 12.3.

x <- c(12, 15, 13, 20, 14, 16, 10, 10, 8,15)
hist(x)

Slika 12.3: Primer histograma.

Ako želimo da konstruišemo promenljivu x, tako da sadrži niz od 1 do 10 sa korakom 0.5, onda možemo koristiti komandu seq. Potom koristeći matematičke operatore možemo dobiti vrednosti promenljive y (y=x²-5x+1) i na kraju napraviti vizualizaciju rezultata na grafiku koristeći funkciju plot (eng. xy scatter plot, Slika 12.4).

x <- seq(1, 10,by=0.5)
y <- x^2- 5*x +1
plot(x,y)

Slika 12.4: Primer xy plota.

Slično, možemo napraviti vizualizaciju trigonometrijskih funkcija koristeći komandu curve, Slika 12.5.

curve (expr = sin, from = 0, to=2*pi)
curve (expr = cos, from = 0, to=2*pi)

Slika 12.5: Sinusna i kosinusna funkcija.

12.2 Osnovne strukture podataka u R-u

Tri osnovna tipa podataka u R-u su: numerički, karakterni i logički tip (eng. numeric, character i logical). Ovi podaci mogu da budu deo različitih kompleksnih objekata (klasa), kao što su objekti kojima su predstavljeni prostorni podaci.

a<-c(1,2,3)
mode(a)
## [1] "numeric"
b<-c("a","b","c")
mode(b)
## [1] "character"
C<-c(T,F,T)
mode(C)
## [1] "logical"

Bilo koji tip podataka (eng. numeric, character, logical) može sadržati NA vrednost (eng. not available). Ta vrednost ukazuje na nedostajući podatak. Mnoge funkcije ne izvršavaju obradu nad podacima dok im se kroz argument ne ukaže da se zanemare NA vrednosti. Često je taj argument u obliku na.rm=T, koji je opcioni argument za mnoge matematičke i statističke funkcije.

max(c(NA, 4, 7))
## [1] NA
max(c(NA, 4, 7),na.rm=T)
## [1] 7
NA | TRUE
## [1] TRUE
NA & TRUE
## [1] NA

Primer manipulacije podacima koji sadrže NA dat je sa ugrađenom funkcijom za ekstrakciju maksimalne vrednosti. Takođe prikazan je i odnos NA podatka logičkim tipom podataka. NULL, za razliku od NA, ne ukazuje na podatak koji nedostaje, već ukazuje na podatak koji ne postoji. Na primer ako hoćemo da definišemo prazan vektor koristili bi vrednost NULL.

Vektori u R okruženju su jedan od osnovnih tipova podataka. U opštem slučaju bilo koji podatak zapisan u R-u direktno je vektor dužine 1. Vektor je niz podataka istog tipa. Detaljnije o tipovima i karakteristikama vektorskih podataka u R-u, videti na http://adv-r.had.co.nz/Data-structures.html#vectors.

U R jeziku simbol # se koristi za komentar, tako da sve što počinje sa # se ne interpretira u R interpreteru.

a<-c(11,27,38)
length(a) # dužina vektora
## [1] 3
summary(a) # sumarna statistika
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 11.00 19.00 27.00 25.33 32.50 38.00
a[2]*8 # 2. element vektora pomnožen sa 8
## [1] 216
as.character(a) # konverzija u tekstualni vektor
## [1] "11" "27" "38"

U ovom primeru kreiran je vektor a sa 3 elementa, prvi element je 11, drugi je 27, treći je 38. Ako želimo da pristupimo drugom elementu vektora, koriste se srednje zagrade, u ovom slučaju kod bi bio sledeći a[2]. Indeksiranje vektora u R jeziku počinje od jedan, dok u nekim drugim jezicima kao što je Java, JavaScript i Python indeksiranje je od nule.

Matrica je pravougaona tabela podataka, gde su svi podaci istog tipa i svaka kolona iste dužine, najčešće to su numerički podaci.

A<-cbind(c(11,2,3),c(8,4,2), c(88,33,11))
## A
## [,1] [,2] [,3]
## [1,] 11 8 88
## [2,] 2 4 33
## [3,] 3 2 11
A<-rbind(c(11,8,88),c(2,4,33),c(3,2,11) )
## A
## [,1] [,2] [,3]
## [1,] 11 8 88
## [2,] 2 4 33
## [3,] 3 2 11

Data su dva primera kreiranja matrice, prvi kad se matrica kreira kolonu po kolonu koristeći komandu cbind i drugi kad se matrica kreira vrstu po vrstu koristeći komandu rbind. Matrica se može kreirati i koristeći komandu matrix.

x.mat <- matrix(c(3, -1, 2, 2, 0, 3, -3, 6), ncol = 2)
## x.mat
## [,1] [,2]
## [1,] 3 0
## [2,] -1 3
## [3,] 2 -3
## [4,] 2 6
x.mat[3,2] # element matrice iz 3. vrste i 2. kolone
## [1] -3
# indeksi matrice mogu dobiti i imena koja nisu samo numerička
dimnames(x.mat) <- list(c("R1","R2","R3","R4"), c("C1","C2"))
## x.mat
## C1 C2
## R1 3 0
## R2 -1 3
## R3 2 -3
## R4 2 6
x.mat["R3","C2"]
## [1] -3

Nizovi se najčešće koriste kao višedimenzionalne matrice ili vektori.

h<-1:24
Z <- array(h, dim=c(3,4,2))
## Z
## , , 1
## [,1] [,2] [,3] [,4]
## [1,] 1 4 7 10
## [2,] 2 5 8 11
## [3,] 3 6 9 1

## , , 2
## [,1] [,2] [,3] [,4]
## [1,] 13 16 19 22
## [2,] 14 17 20 23
## [3,] 15 18 21 24

12.2.1 Klasa data.frame

Najčešće korišćena klasa podataka u R-u je data.frame. Objekat tipa data.frame se može posmatrati kao pravougaona tabela, gde svaka kolona može biti međusobno različitog tipa, ali sve kolone su iste dužine. Unutar svake kolone pojedinačno su podaci istog tipa, jer je praktično jedna kolona jedan vektor. Ova struktura podataka najviše podseća na tabelu iz baze podataka.

Dat je primer data.frame objekta koji je nastao kao kompozicija 3 vektora i to tipa: karakter, numeričkog i logičkog.

student = c("Nikola", "Jelena", "Marko")
god_stud = c(3, 3, 1)
ranije_progr = c(TRUE, FALSE, TRUE)
df = data.frame(student , god_stud, ranije_progr)
## df
## student god_stud ranije_progr
## 1 Nikola 3 TRUE
## 2 Jelena 3 FALSE
## 3 Marko 1 TRUE
fix(df)

Korišćenjem funkcije data.frame je napravljen objekat tipa data.frame i sačuvan je u promenljivoj df. Pozivom komande str može se videti struktura objekta df. On sadrži tri kolone (atributa) koje imaju nazive student, god_stud i ranije_prog.

str(df)
## 'data.frame': 3 obs. of 3 variables:
## $ student : Factor w/ 3 levels "Jelena","Marko",..: 3 1 2
## $ god_stud : num 3 3 1
## $ ranije_progr: logi TRUE FALSE TRUE

Možemo videti da je kolona (atribut objekta df) koja je nastala od vektora student promenila klasu iz tipa karakter u tip faktor (eng. factor). To je uobičajen način za čuvanje kategorijskih podataka u R okruženju, a funkcija data.frame ima podrazumevani argument da podatke tipa karakter prebacuje u kategorijske. Time se proširuje mogućnost korišćenja statističkih analiza nad podacima u odnosu na to kad su podaci tipa karakter. U osnovi za svaki unikatni element vektora koji je tipa faktor se pridružuje broj, gde je svaki broj oznaka kategorije.

Postoji nekoliko načina da selektujemo celu kolonu iz data.frame objekta, sledeći primeri to pokazuju.

df$student # selekcija po imenu atributa sa operatorom $
## [1] Nikola Jelena Marko
## Levels: Jelena Marko Nikola
df[,1] # selekcija po red. Br. atr., kao indeks kod matrice
## [1] Nikola Jelena Marko
## Levels: Jelena Marko Nikola
df[,'god_stud'] # selekcija po imenu atributa kao indeksom
## [1] 3 3 1

Prva selekcija je koristeći operator $ i ime kolone, druga je koristeći redni broj kolone (slično kao kod matrice) i treća koristeći ime kolone u okviru srednjih zagrada (kao kod matrice kod koje su definasana dinamička imena kolona). Za selekcije pojedinačnog elementa, na primer za element druge vrste i treće kolone koristi se isti način selekcije kao da je u pitanju matrica (df [2,3]).

Objekat tipa df omogućava i atributne selekcije. Kod prostornih vektorskih podataka atributi su smešteni u data.frame objekte, tako da ovaj način atributne selekcije važi i za prostorne podatke.

Ako želimo da selektujemo sve vrste za koje je ispunjen uslov da je ime studenta Nikola, napravićemo sledeću sintaksu za upit po vrstama.

df[df$student=='Nikola',]
## student god_stud ranije_progr
## Nikola 3 TRUE

Slično možemo zadati upit da nam se prikažu svi podaci za studente koji su starija godina studija od prve.

df[df$god_stud>1,]
## student god_stud ranije_progr
## 1 Nikola 3 TRUE
## 2 Jelena 3 FALSE

Ili možemo zadati upit po više kriterijuma, da su na višoj godini od prve i da su ranije pohađali programiranje.

df[df$god_stud&gt;1 & df$ranije_progr==T,]
## student god_stud ranije_progr
## 1 Nikola 3 TRUE

12.2.2 Liste

Lista je uredjen skup podataka koji ne moraju biti istog tipa i iste dužine, tako da predstavlja vrlo fleksibilnu klasu za skladištenje podataka. Na primer lista =[vektor, data.frame, lista].

kurs.l <- list(ime="R kurs", br_ucesnika=3, podaci=df)
## $ime
## [1] "R kurs"
## $br_ucesnika
## [1] 3
## $podaci
## student god_stud ranije_progr
## 1 Nikola 3 TRUE
## 2 Jelena 3 FALSE
## 3 Marko 1 TRUE

Kod liste za selekciju po elementu koriste se dvostruke uglaste zagrade ili operator $.

# Selekcije
kurs.l$ime # po elementu liste koji se zove ime
kurs.l[[3]] # selekcija 3. elementa liste
kurs.l[[3]]$student
# selekcija 3. elementa liste koji je tipa data.frame,
# plus dodatna selekcija po koloni student

12.2.3 Objekti i klase u R jeziku

Osim postojećih klasa, bilo da su u pitanju osnovne ili klase definisane u nekom od paketa, korisnik može da koristi i sopstvene klase i metode. Postoje dve vrste klasa u R-u i to:

S3 klase – stare, neformalne klase,
S4 klase – formalne klase,
R5 klase - R Reference Class.

Primer kreiranja klase i metode za vizualizaciju nad klasom deca, koristeći stari S3 sistem.

# kreiranje klase deca
deca <- function(ID, p, u, v, t){
x <- list(ID=ID,
Pol=p,
Podaci=data.frame(Uzrast.m=u, Visina.cm=v, Tezina.kg=t)
)
class(x) <- 'deca'
invisible(x)
}
# Kreiranje metode nad klasom deca
plot.deca <- function(object){
Podaci <- object$Podaci
par(mfrow=c(1,2))
plot(Podaci$Uzrast, Podaci$Visina.cm, type="b", pch=19,col="red",
xlab="Uzrast (mjeseci)", ylab="Visina (cm)", main="Odnos visine i uzrasta")
plot(Podaci$Uzrast, Podaci$Tezina.kg, type="b", pch=19,col="blue",
mtext(paste("ID",object$ID,",",object$Pol,sep=""), side=3,outer=TRUE,line=-1.5,cex=1.5)
}

U prvom delu je kreirana klasa deca, gde je koristeći klasu liste definisano kako su struktuirani podaci. Potom je kreirana metoda koja će se primenjivati nad klasom deca za vizualizaciju podataka. Praktično metoda kreira dva grafika od kojih jedan predstavlja odnos uzrasta i visine, a drugi odnos uzrasta i telesne mase deteta, Slika 12.6.

Kreiranje istance i primena metode

uzrast = c(0,3,6,12,18,24,30,36)
decaci.t =c(3.53,6.39,8.16,10.46,11.80,12.74,13.56,14.33)
decaci.v = c(49.99,62.08,67.86,76.11,82.41,87.66,92.13,95.45)
# Instanca klase (kreiranje obj. x klase deca)
x <- deca(1, "M", uzrast, decaci.v, decaci.t)
str(x)
## List of 3
## $ ID : num 1
## $ Pol : chr "M"
## $ Podaci:'data.frame': 8 obs. of 3 variables:
## ..$ Uzrast.m : num [1:8] 0 3 6 12 18 24 30 36
## ..$ Visina.cm: num [1:8] 50 62.1 67.9 76.1 82.4 ...
## ..$ Tezina.kg: num [1:8] 3.53 6.39 8.16 10.46 11.8 ...
## - attr(*, "class")= chr "deca"
## class(x)
## [1] "deca"
## 
## plot(x)

Slika 12.6: Metoda plot za klasu deca

Primer kreiranje objekta koristeći S4 klase i metode dat je bez detaljnog objašnjena.

setClass("deca",representation(ID="numeric",
Pol="character",Podaci = "data.frame"))
# Kreiranje metoda za S4 klasu deca
setGeneric("plot")
setMethod ("plot" , "deca",
function (x , xlab ="Uzrast (mjeseci)", ylab="Visina (cm)" , main="Odnos visine i uzrasta", pch=19, ...){
Podaci=x@Podaci
par(mfrow=c(1,2))
plot(Podaci[,1], Podaci[,3], type="b", pch=19,col="red",
xlab="Uzrast (meseci)", ylab="Visina (cm)", main="Odnos visine i uzrasta", ...)
plot(Podaci[,1], Podaci[,2], type="b", pch=19,col="blue",xlab="Uzrast (mjeseci)",
     ylab="Težina (kg)", main="Odnos težine i uzrasta", ...)
mtext(paste("ID",x@ID,",",x@Pol,sep=""), side=3,outer=TRUE,line=-1.5,cex=1.5)
})
# Kreiranje instance objekta deca
uzrast = c(0,3,6,12,18,24,30,36)
decaci.t =c(3.53,6.39,8.16,10.46,11.80,12.74,13.56,14.33)
decaci.v = c(49.99,62.08,67.86,76.11,82.41,87.66,92.13,95.45)
x <- new("deca", ID=1, Pol="M", Podaci=data.frame(uzrast, decaci.t, decaci.v))
# Struktura objekta
str(x)
## Formal class 'deca' [package ".GlobalEnv"] with 3 slots
## ..@ ID : num 1
## ..@ Pol : chr "M"
## ..@ Podaci:'data.frame': 8 obs. of 3 variables:
## .. ..$ uzrast : num [1:8] 0 3 6 12 18 24 30 36
## .. ..$ decaci.t: num [1:8] 3.53 6.39 8.16 10.46 11.8 ...
## .. ..$ decaci.v: num [1:8] 50 62.1 67.9 76.1 82.4 ...
plot(x)

R5 (RC - Reference classes) je slična S4 klasi, ali sa pridruženim okruženjem, vrlo sličnog pristupa kao kod drugih programskih jezika: Python, Ruby, Java i C# .

Tabela 12.2: Funkcije za ispitivanje strukture podataka.
`class`	Klasa objekta. (vector, matrix, function, logical, list, … )
`str`	Struktura podataka
`mode`	Tip podataka. (Numeric, character, logical, …)
`typeof`	Tip koji koristi R da bi skladištio podatak u memoriji (double, integer, character, logical, …)
`is.function`	Da li je funkcija(TRUE if function)
`is.na`	Da li je NA (TRUE if missing)
`names`	Imena pridružena objektu
`dimnames`	Imena za indekse kod vektora, matrica i redova
`slotNames`	Imena slotova – delova objekta (npr SP podaci)
`attributes`	Imena i klase atributa kod objekta…

12.3 Neke funkcije za uvoz i izvoz podataka u R

U R-u postoji više funkcija za čitanje i pisanje podataka iz fajlova i obratno kao što su:

readLines, writeLines;
scan, write;
read.csv, write.csv;
read.csv2, write.csv2;
read.delim, write.delim;
read.delim2, write.delim2;

Takođe, postoji niz paketa koji omogućavaju direktnu komunikaciju sa bazama podataka kao što su:

RODBC,
RMySQL,
RpostGIS,
RRasdaman,
RPostgreSQL,
ROracle,
RJDBC, i drugi.

12.4 Kontrola toka i petlje

Slično kao i u drugim jezicima, za kontrolu toka u R-u koristi se komanda if, koja generalno usmerava interpreter na to šta se izvršava ako je neki logički uslov zadovoljen, a šta ako nije. Funkcionisanje if komande prezentovano je u obliku koda.

if (logički uslov) {
izrazi (statements)
}
else {
alternativni izrazi
}
# { } su opcione u slučaju samo jednog izraza
# za jednostavne slučajeve
ifelse (logički uslov, da izraz, ne izraz)

Funkcionisanje if komande ilustrovano je i sledećim primerom.

x<- runif(100, min=0, max=5)
x=mean(x)
if(x>2.5){
cat('x je vece od 2.5 \n') # \n novi red
Y=2*x
} else {
cat('x je manje od 2.5 \n')
Y=3*x
}
# primena ifelse komande za proste slučajeve kontrole toka
Y<-ifelse(x>2.5, 2*x, 3*x)

Prvo je kreiran vektor čiji su elementi generisani primenom uniformne raspodele sa minimalnom vrednošću 0 i maksimalnom 5. Potom je promenljiva x postala srednja vrednost gore pomenutog vektora. Zatim je definisan logički uslov i set komandi koje se izvršavaju kad je uslov zadovoljen i kad uslov nije zadovoljen.

Petlje u R-u su implementirane slično kao i u drugim jezicima. Dat je primer for i while komande bez dodatnog objašnjenja.

for(i in 1:10) {
print(i*i+2*i)
}
i<-1
while(i<=10) {
print(i*i)
i<-i+sqrt(i)
}
# Pogledati: repeat, break, next

Funkcije se kreiraju koristeći komandu function. Funkcijama kreiramo neku proceduru nad podacima, vrši se obrada ulaznih podataka i dobija rezultat. Dat je primer kreiranja funkcije koja izračunava geometrijsku sredinu.

# kreiranje funkcije
gs <- function(a,b) {
result <- sqrt(a*b)
return(result)
}
gs(7,88) # poziv funkcije
## [1] 24.8193472919817

Funkcije lapply, sapply, apply su ugrađene funkcije u R-u koje olakšavaju manipulaciju nad vektorima, listama i matricama. Izbegava se kontrola elementa koristeći petlje, u ovom slučaju umesto petlji se koriste funkcije u okviru koje pored podataka kao argumenata, uzima se i funkcija nad podacima kao agrument. Najčešće se ubrzava proces obrade podataka, ali ne uvek.

A<-cbind(c(11,2,3) , c(8,4,2) , c(88,33,11) )
A
## [,1] [,2] [,3]
## [1,] 11 8 88
## [2,] 2 4 33
## [3,] 3 2 11
apply(A, 1, sum) # suma po vrstama
## [1] 107 39 16
apply(A, 2, sum) # suma po kolonama
## [1] 16 14 132
apply(A, 2, sort) # sortiranje po kolonama
## [,1] [,2] [,3]
## [1,] 2 2 11
## [2,] 3 4 33
## [3,] 11 8 88

U predhodnom primeru je dato kako komandu apply možemo koristiti da na jednostavan način obrađujemo delove matrice po vrstama i kolonama, a da pri tom ne koristimo for petlju.

Sledeći primer ilustruje korišćenje lapply i sapply funkcija koje se primenjuju nad listama ili vektorima. Funkcija lapply uvek kao odgovor daje listu, dok sapply kao odgovor daje vektor.

Slika 12.7: Ilustruja korišćenja lapply i sapply

U prethodnom primeru primenjen je proračun srednje vrednosti nad svakim elementom liste. Kod primene lapply funkcije kao odgovor je dobijena lista sa rezultatima, a sapply daje vektor. Slično je ilustrovano i sa vektorom txt i funkcijom koja broji karaktere, iako je ulazni podatak bio vektor izlaz lapply funkcije je lista, a sapply je opet vektor.

12.5 Prostorni podaci u R-u

12.5.1 Paket sp

Ranije su mnogi paketi sadržali svoju reprezentaciju prostornih podataka. Takav pristup je otežavao razmenu podataka između paketa, kao i sa eksternim aplikacijama. Paket sp (Pebesma and Bivand 2018) je standardizovao reprezantaciju prostornih podataka u R-u i omogućio jednostavnu konverziju u druge tipove unutar R-a ili ka spoljnim aplikacijama. Paket sp obezbeđuje reprezentaciju i vektorskih i rasterskih prostornih podataka.

Sve sp klase prostornih podataka su prikazane u sledećem primeru.

library(sp)
getClass("Spatial")
## Class "Spatial" [package "sp"]
## 
## Slots:
## 
##   Name:         bbox proj4string
## Class:      matrix         CRS
## 
## Known Subclasses:
##   Class "SpatialPoints", directly
## Class "SpatialMultiPoints", directly
## Class "SpatialGrid", directly
## Class "SpatialLines", directly
## Class "SpatialPolygons", directly
## Class "SpatialPointsDataFrame", by class "SpatialPoints", distance 2
## Class "SpatialPixels", by class "SpatialPoints", distance 2
## Class "SpatialMultiPointsDataFrame", by class "SpatialMultiPoints", distance 2
## Class "SpatialGridDataFrame", by class "SpatialGrid", distance 2
## Class "SpatialLinesDataFrame", by class "SpatialLines", distance 2
## Class "SpatialPixelsDataFrame", by class "SpatialPoints", distance 3
## Class "SpatialPolygonsDataFrame", by class "SpatialPolygons", distance 2

12.5.1.1 Klasa SpatialPointsDataFrame

Klasa Spatial je osnovna klasa koja je sadržana u svim drugim podklasama. Spatial klasa sadrži dve osobine (u R-u se osobina naziva slot) i to bbox - prostorni obuhvat i proj4string slot koji sadrži opis koordinatnog referentnog sistema u PROJ.4 notaciji. Na primer, SpatialPoints pored ova dva slota sadrži i treći slot koji predstavlja matricu koordinata. Klasa SpatialPointsDataFrame pored pomenuta 3 slota koja sadrži podklasa SpatialPoints sadži i slot data, gde su čuvaju atributni podaci kao data.frame objekat. Slot coords.nrs se popunjava samo ako su tačkasti podaci nastali iz tabele (data.frame objekta). Tada se u tom slotu čuva infomacija o poziciji x i y koordinate u tabeli, pa obrnuta konverzija iz SpatialPointsDataFrame u data.frame daje istu ulaznu tabelu.

Slika 12.8: Struktura SpatialPointsDataFrame podklase i njen odnos sa SpatialPoints podklasom i Spatial klasom.

Primer kreiranja SpatialPointsDataFrame objekta iz data.frame objekta demonstriran je koristeći set podataka meuse, koji je u obliku tabele podataka (data.frame objekat), a podaci su sadržani u paketu sp. Više o samim podacima može se naći pozivajući komandu ? meuse u komandnom prozoru. Podaci sadrže koordinate za 155 tačaka i podatke o uzorkovanim karakteristikama zemljišta u tim tačkama.

library(sp)
data(meuse)
str(meuse)
## 'data.frame':    155 obs. of  14 variables:
##   $ x      : num  181072 181025 181165 181298 181307 ...
## $ y      : num  333611 333558 333537 333484 333330 ...
## $ cadmium: num  11.7 8.6 6.5 2.6 2.8 3 3.2 2.8 2.4 1.6 ...
## $ copper : num  85 81 68 81 48 61 31 29 37 24 ...
## $ lead   : num  299 277 199 116 117 137 132 150 133 80 ...
## $ zinc   : num  1022 1141 640 257 269 ...
## $ elev   : num  7.91 6.98 7.8 7.66 7.48 ...
## $ dist   : num  0.00136 0.01222 0.10303 0.19009 0.27709 ...
## $ om     : num  13.6 14 13 8 8.7 7.8 9.2 9.5 10.6 6.3 ...
## $ ffreq  : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
## $ soil   : Factor w/ 3 levels "1","2","3": 1 1 1 2 2 2 2 1 1 2 ...
## $ lime   : Factor w/ 2 levels "0","1": 2 2 2 1 1 1 1 1 1 1 ...
## $ landuse: Factor w/ 15 levels "Aa","Ab","Ag",..: 4 4 4 11 4 11 4 2 2 15 ...
## $ dist.m : num  50 30 150 270 380 470 240 120 240 420 ...
# Konvertovanje tabelarnih podataka u prostorne
coordinates(meuse) <- ~x+y
proj4string(meuse) <- CRS("+init=epsg:28992")

Komandom coordinates se konvertuje objekat data.frame u SpatialPointsDataFrame, neophodno je ukazati na kolone koje predstavljaju koordinate (<-~x+y, u ovom slučaju kolone se nazivaju x i y). Tako se dobija prostorni objekat, kome se izmenom proj4string slota može zadati odgovarajući koordinatni referentni sistem. Za ove podatke to je stereografska projekcija za Holandiju (više na http://epsg.io/28992).

Struktura podataka sa svim slotovima sa detaljama do drugog niovoa je data za novoformirani objekat meuse:

str(meuse, max.level=2)
## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
## ..@ data       :'data.frame':    155 obs. of  12 variables:
##   ..@ coords.nrs : int [1:2] 1 2
## ..@ coords     : num [1:155, 1:2] 181072 181025 181165 181298 181307 ...
## .. ..- attr(*, "dimnames")=List of 2
## ..@ bbox       : num [1:2, 1:2] 178605 329714 181390 333611
## .. ..- attr(*, "dimnames")=List of 2
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot

Selekcija atributnih podataka je ista kao i kod data.frame objekta, sledećim primerom selektujemo prve tri tačke sa 3 atributa.

meuse@data[1:3,c("zinc","dist.m","elev")]
## zinc dist.m  elev
## 1 1022     50 7.909
## 2 1141     30 6.983
## 3  640    150 7.800

Vizualizacija prostornih podataka može se kreirati koristeći standardne funkcije kao što je plot i druge ili metode koje su napravljene za sp klase kao što su spplot ili bubble. Metod spplot je standardan metod za vizualizaciju svih tipova vektorskih i rasterskih podataka tipa sp, dok je bubble metod za geovizualizaciju tačkastih podataka koristeći metod proporcionalnih simbola u zavisnosti od vrednosti atributa koji se vizualizuje. U narednom primeru bubble metod je korišćen za izradu karte koristeći simbole proporcionalne atributnom podatku o koncentraciji cinka u zemljištu (atribut zinc, dat je ppm jedinicama), argumentom main zadat je naslov, Slika 12.10.

spplot(meuse,
       "zinc",
       key.space= list(x=0.2,y=0.9, corner=c(0,1),
       scales=list(draw=T),
       col.regions=rainbow(5),
       main="Cink")

Slika 12.9: Metoda spplot i metoda boja.

bubble(meuse, "zinc", main="Cink" ,maxsize = 1.5)

Slika 12.10: Metoda bubble: proporcionalni simboli.

Metod spplot se najčešće koristi za izradu karata gde se vizualizuju atributni podaci metodom boja, Slika 12.9. Argumentom col.regions zadata je paleta boja, koristeći ugrađenu paletu rainbow, definisan je položaj legende unutar plota (key.space), zadat naslov (main), iscrtavanje koordinatnih osa (scales).

Geovizualizaciju je moguće praviti i koristeći standardne plot funkcionalnosti. Sledeći primer ilustruje izradu karte koristeći prikaz proporcionalnim simbolima, i standardan plot pristup, , Slika 12.11.

plot(meuse, pch = 1, cex = sqrt(meuse$zinc)/12, axes = TRUE)
v = c(100,200,400,800,1600)
legend("topleft", legend = v, pch = 1, pt.cex = sqrt(v)/12)

Argument pch koristi se za definisanje oblika simbola, u ovom slučaju izabran je kružni simbol (1-kružnica). Veličina simbola je proporcionalna korenu koncentracije cinka u tački podeljenom sa brojem 12 (argument cex), takođe je zadato iscrtavanje koordinatnih osa (argument axes=TRUE). Zatim je definisana legenda koristeći komandu legend i odgovarajuće argumente kao u komandi plot.

Slika 12.11: Metoda plot: proporcionalni simboli.

12.5.1.2 Klasa SpatialLinesDataFrame

Klasa SpatialLines je nešto komlikovanija nego klasa koja se odnosi na tačke, u suštini može da podrži geometriju tipa linija ili više linija. Sledeći kod definiše kreiranje jednog objekta tipa SpatialLines. Kod je dat bez detaljnog objašnjenja.

l1 = cbind(c(1, 2, 3), c(3, 2, 2))
l1a = cbind(l1[, 1] + 0.05, l1[, 2] + 0.05)
l2 = cbind(c(1, 2, 3), c(1, 1.5, 1))
Sl1 = Line(l1)
Sl1a = Line(l1a)
Sl2 = Line(l2)
S1 = Lines(list(Sl1, Sl1a), ID = "a")
S2 = Lines(list(Sl2), ID = "b")
Sl=SpatialLines(list(S1,S2))
summary(Sl)
## Object of class #SpatialLines
## Coordinates:
##   min  max
## x   1 3.05
## y   1 3.05
## Is projected: NA 
## proj4string : [NA]

Objektu tipa SpatialLines možemo dodati i atributne podatke koji su objekat tipa data.frame, tako se dobija objekat SpatialLinesDataFrame, koji predstavlja vektorske podatke linijskog tipa.

df = data.frame(ime_linije= c("prva lin.", "druga lin."))
Sldf = SpatialLinesDataFrame(Sl, data = df, mach.ID=F)
plot(Sl, col = c("red", "blue"),axes=T,type="b",lwd=2)
grid(nx=3, ny=3, lwd=1,lty=2,col="green")

Slika 12.12: Metoda plot linije.

df = data.frame(ime_linije= c("prva lin.", "druga lin."))
Sldf = SpatialLinesDataFrame(Sl, data = df, mach.ID=F)
plot(Sl, col = c("red", "blue"),axes=T,type="b",lwd=2)
grid(nx=3, ny=3, lwd=1,lty=2,col="green")
spplot(Sldf, col.regions=c('red','blue'),scales=list(draw=T))

Slika 12.13: Metoda spplot linije.

Objekat Sldf je kartiran koristeći spplot (Slika 12.12) i plot metod. Kod plot metoda je definisan tip plota, tip linije, boje linije, gridne linije, tj. koordinatne linije, iscrtavanje koordinatnih osa, Slika 12.13.

Dat je još jedan primer linijskih podataka koristeći paket geosphere (R. J. Hijmans 2017 a) kojim je generisana geodetska linija (komanda greatCircle), koja sadrži približne koordinate Beograda i Njujorka. Sračunate su i dužine geodetske linije (komanda distGeo) i loksodrome (komanda distRhumb) od Beograda do Njujorka. Iscrtane su granice država, potom geodetska linija, a na kraju su prikazane pozicije Beograda i Njujorka kao crveni simboli u obliku zvezde, (Slika 12.14). U ovim primerima linijski podaci su dati kao standardne matrice, a ne strukturirani prostorni podaci.

library(geosphere)
NY <- c(-73.78, 40.63)
BG <- c(20.44, 44.78)
# Dužina geodetske linije i loksodrome
distGeo(BG, NY); distRhumb(BG, NY)
## 7271839 7716831
gc <- greatCircle(BG,NY)
data(wrld)
plot(wrld, type='l')
lines(gc, lwd=2, col='blue')
points(rbind(BG,NY), pch='*', cex=3, col='red')

Slika 12.14: Metoda plot geodetska linija i tačke.

12.5.1.3 Klasa SpatialPolygonsDataFrame

Poligonski podaci sa atributima su podaci iz klase SpatialPolygonsDataFrame. Ovde se neće objašnjavati struktura sp objekta poligonskog tipa. U sledećem primeru atributni poligonski podaci su učitani direktno iz ESRI Shape fajla, koristeći R paket (Bivand, Keitt, and Rowlingson 2018). Paket rgdal direktno konvertuje kako rasterske tako i vektorske podatke u odgovarajuće sp objekte. Takođe, ako su u fajlu sadržani podaci o koordinatnom referentnom sistemu, ti podaci se upisuju u odgovarajući sp slot. Paket rgdal može se koristiti i za eksportovanje sp objekata u GDAL/OGR formate. Detaljni opis paketa i funkcija može se videti u manualu rgdal paketa.

Sledeći kod učitava rgdal paket, potom se preuzima fajl (komanda download.file) koristeći Web Feature Service (WFS) url upit (sačuvan u promenljivoj url) i na kraju koristeći funkcionalnosti rgdal paketa (komanda ogrInfo) zahteva prikaz osnovnih informacija o preuzetom fajlu (poage.gml).

library(rgdal)
url=paste('http://osgl.grf.bg.ac.rs/geoserver/osgl_4/wfs',
          '?service=WFS&version=1.0.0&',
          'request=GetFeature&typeName=osgl_4:POP_age',sep='')
download.file(url, 'poage.gml')
ogrInfo('poage.gml')
## Source: "poage.gml", layer: "POP_age"
## Driver: GML; number of rows: 197
## Feature type: wkbPolygon with 2 dimensions
## Extent: (18.81499 41.8521) - (23.00637 46.19006)
## Number of fields: 14
## name type length typeName
## 1 gml_id 4 0 String
## 2 ID_Municip 0 0 Integer
## ...

Potom se učitavaju podaci koristeći komandu readOGR, koja ima jedan obavezan argument, a to je ime fajla i putanja do fajla (ako je fajl u radnom direktorijumu putanja nije neophodna). U ovoj komandi definisani su i dodatni argumenti i to ime lejera i opcija da se ne konvertuju atributi tipa string u tip factor. Potom je prikazana struktura podataka sačuvanih u objektu pa. Podaci predstavljaju granice opština u Srbiji sa pridruženim atributnim podacima iz podataka Popisa stanovništva, Slika 12.15. Atribut TOTAL_POP koji predstavlja ukupan broj stanovnika po opštini prema podacima iz Popisa u Srbiji 2011, konvertovan je u numerički tip.

pa <- readOGR(dsn ='poage.gml', layer = 'POP_age', stringsAsFactors=FALSE )
str(pa, max.level=2)
## Formal class 'SpatialPolygonsDataFrame'
## ..@ data :'data.frame': 197 obs. of 14 variables:
## ..@ polygons :List of 197
## ..@ plotOrder : int [1:197] 143 10 20 171 170 168 121 169 152 3 ...
## ..@ bbox : num [1:2, 1:2] 18.8 41.9 23 46.2
## .. ..- attr(*, "dimnames")=List of 2
## ..@ proj4string:Formal class 'CRS'
pa$TOTAL_POP <-as.numeric(pa$TOTAL_POP)
spplot(pa, zcol = 'TOTAL_POP', scales=list(draw=T), 
       main="Broj stanovnika po opštinama popis 2011")

Slika 12.15: Metoda spplot: horoplet karta.

Predhodna karta dobijena je koristeći spplot metod. Dobijena je horoplet karta broja stanovnika po opštinama. Podaci su dati u WGS84 sistemu, a poznato je da nije uobičajeno da se statistički podaci prikazuju koristeći elipsoidne koordinate, već se kartiraju u ekvivalentnim projekcijama.

Koristeći podatke o popisu po opštinama (objekat sačuvan u promenljivoj pa) napravićemo kartu gustine stanovništva. Prvo transformišemo prostorne podatke iz jednog koordinatnog sistema u drugi koristeći komandu spTransform. Transformacija je iz WGS84 u Lambertovu azimutnu ekvivalentnu projekciju sa dodirnom paralelom 44 stepena severno, centralnim meridijanom 21 stepen istočno, sa zadatim koordinatma koordinatnog početka u ravni i datumom WGS84.

### dodeljivanje CRS #######################################
proj4string(pa) <- CRS("+init=epsg:4326")
# Transformacija iz WGS u LAEA
pa <- spTransform(pa, CRS('+proj=laea +lat_0=44 +lon_0=21 +x_0=200000 +y_0=200000 +datum=WGS84'))

Koristeći R paket rgeos (Bivand and Rundel 2017), koji je interfejs ka GEOS biblioteci otvorenog koda, sračunate su površine (komanda gArea) svakog poligona, tj. opštine, redom (komanda sapply). Površine su podeljene sa milion da bi bile izražene u km² i dodeljene objektu pa kao nov atribut, a zatim sračunata gustina po km².

### računanje površine u LAEA ###############################
library(rgeos)
pa$povrsina <-sapply(1:nrow(pa), function(i) gArea(pa[i,])/1000000)
pa$gustina <- pa$TOTAL_POP/pa$povrsina

Potom su definsani intervali u zavisnosti od vrednosti gustine po opštinama. R paket classInt (Bivand 2017), sadrži niz metoda za klasifikaciju atributa za potrebe geovizualizacije. U narednom primeru korišćen je kvantil metod za klasifikaciju podataka u zavisnosti od vrednosti gustine. Kvantil intervali dele atributne podatke na jednake grupe (isti broj opština se nalazi u svakoj grupi), opseg podele varira, tj nije konstantna širina opsega. Za izbor boja je korišćen paket RColorBrewer (Neuwirth 2014), koji obezbeđuje predefinisane palete boja za geovizualizaciju sugerisane od strane američke profesorice kartografije Cynthia Brewer (Brewer, Hatchard, and Harrower 2003), http://colorbrewer2.org.

######### intervali i boje #############################
library(classInt)
library(RColorBrewer)
plotvar <- pa$gustina
nclr <- 8
plotclr <- brewer.pal(8,"YlOrRd")
class <- classIntervals(plotvar, nclr, style="quantile", dataPrecision=1)
colcode <- findColours(class, plotclr)

Nakon kvantil klasifikacije podataka o gustini stanovništva po opštinama generisana je karta koristeći standardne plot metode. Prvo je generisana koordinta mreža u projekciji, koristeći paket graticule (Sumner 2016). Potom je generisana horoplet karta sa iscrtanom koordinatnom mrežom, naslovom, legendom u skladu sa klasifikacijom gustine kvantil metodom i izabranim bojama, Slika 12.16.

######### kreiranje karte ##############################
####### gridne linije u projekciji #########
library(graticule)
lons <- 18:23 ; lats <- 42:46
xl <- range(lons) + c(-0.4, 0.4)
yl <- range(lats) + c(-0.4, 0.4)
grat <- graticule(lons, lats, 
                  proj = proj4string(pa),
                  xlim = xl, ylim = yl) # SpatialLinesDataFrame
labs <- graticule_labels(lons, lats, xline = min(xl),
                         yline = max(yl), 
                         proj = proj4string(pa)) # SpatialPointsDataFrame
#######################
op <- par(mar = c(0,0,2,2))
plot(extent(grat) + c(4, 2) * 1e5, asp = 1, type = "n", axes = FALSE, xlab = "", ylab = "")
plot(pa, add = TRUE)
plot(pa, col=colcode, add=TRUE)
title(main="Gustina stanovnika po opštinama",
sub="Klase su u kvantilima")
plot(grat, add = TRUE, lty = 5, col = rgb(0, 0, 0, 0.8))
text(subset(labs, labs$islon), lab = parse(text = labs$lab[labs$islon]), pos = 3)
text(subset(labs, !labs$islon), lab = parse(text = labs$lab[!labs$islon]), pos = 2)
legend("right", legend=names(attr(colcode, "table")),
fill=attr(colcode, "palette"), cex=0.8, bty="n")

Slika 12.16: Metoda plot: horoplet karta.

Kod horoplet karata je veoma značajan način klasifikacije atributnih podataka (takođe važi i za druge primere vizualizacije kvantitativnih veličina). Sledeći kod demonstrira tri metode klasifikacije podataka za promenljivu koja predstavljaja udeo odraslih stanovnika u odnosu na ukupnu populaciju. Pored pomenute kvantil metode, korišćena je i metoda jednakih intervala (prekida) i metod Dženksovih (eng Jenks) intervala (prekida). Kod metoda jednakih intervala opseg podataka se deli na jednake delove, zavisno od raspodele podataka vrlo često se dešava da je broj entiteta u intervalima vrlo neravnomerno raspoređen, ovaj primer isto to jasno demonstrira. Dženksovi intervali se formiraju tako da su klase dobijene metodom optimizacije, i to sa uslovom da je što manja varijansa unutar klasa, a što veća među različitim klasama.

library(classInt)
library(RColorBrewer)
pal = brewer.pal(7,"Oranges")
pa$ADULTS <- as.numeric(pa$ADULTS)
pa$deoStarijih <- 100 * pa$ADULTS/pa$TOTAL_POP
brks.qt = classIntervals(pa$deoStarijih, n = 7, style = "quantile", dataPrecision = 1)
brks.jk = classIntervals(pa$deoStarijih, n = 7, style = "jenks", dataPrecision = 1)
brks.eq = classIntervals(pa$deoStarijih, n = 7, style = "equal", dataPrecision = 1)

spplot(pa, "deoStarijih",at=brks.eq$brks, col.regions=pal, col="transparent",
main = list(label="Jednaki intervali"))
plot(brks.eq,pal=pal,main="Jednaki intervali")

spplot(pa, "deoStarijih",at=brks.qt$brks,col.regions=pal, col="transparent",
main = list(label="Kvantil intervali"))
plot(brks.qt,pal=pal,main="Kvantil intervali")

spplot(pa, "deoStarijih",at=brks.jk$brks,col.regions=pal, col="transparent",
main = list(label="Dženksovi intervali"))
plot(brks.jk,pal=pal,main="Dženksovi intervali")

Gore navedeni kod pravi podelu na klase koristeći 3 metoda za izbor intervala, potom se kreiraju horoplet karte koristeći spplot metod. Uz svaku horoplet kartu kreiran je i grafik gustine raspodele u odnosu na definisane klase (Slike: 12.17, 12.18, 12.19).

Slika 12.17: Jednaki intervali.

# Jednaki intervali
# interval , broj opština u intervalu
[69.2,72.3)  2
[72.3,75.5)  2
[75.5,78.7)  2
[78.7,81.9)  38
[81.9,85.1)  109
[85.1,88.3)  13
[88.3,91.5]  2
# Kvantil intervali
# interval , broj opština u intervalu
[69.1,81.3) 24
[81.3,81.9) 24
[81.9,82.4) 24 
[82.4,83)  24
[83,83.6)  24
[83.6,84.6) 24
[84.6,91.5] 24
# Dženksovi intervali
# interval , broj opština u intervalu
[69.1,70.8] 2
(70.8,77.6] 4
(77.6,81.4] 23
(81.4,82.7] 58
(82.7,84] 42
(84,85.8] 31
(85.8,91.5] 8

Slika 12.18: Kvantil intervali.

Slika 12.19: Dženksovi intervali.

12.5.1.4 Klasa SpatialPixelsDataFrame i SpatialGridDataFrame

Paket sp čuva rasterske podatke u dve klase. Tipični rasterski podaci mogu imati n atributa, i oni su smešteni u strukturu SpatialGridDataFrame. Ova klasa je struktuirana slično kao i rasterski podaci u bilo kom formatu, georeferencija je obezbeđena topologijom grida i koordinatama početnog piksela. Druga sp klasa za rasterske podatake je SpatialPixelsDataFrame, ona je vrlo slična tačkastim podacima gde se čuvaju koordinate svakog piksela, moguća je direktna konverzija u tip SpatialPointsDataFrame i u SpatialGridDataFrame.

Primer kreiranja objekta SpatialPixelsDataFrame iz tabele podataka meuse.grid (podaci su deo sp paketa) je dat u narednom primeru.

data(meuse.grid)
coordinates(meuse.grid) = \~x+y
proj4string(meuse.grid) <- CRS("+init=epsg:28992")
gridded(meuse.grid) = TRUE
str(meuse.grid, max.level = 2)
## Formal class 'SpatialPixelsDataFrame' [package "sp"] with 7 slots
## ..@ data :'data.frame': 3103 obs. of 5 variables:
## ..@ coords.nrs : num(0)
## ..@ grid :Formal class 'GridTopology' [package "sp"] with 3 slots
## ..@ grid.index : int [1:3103] 69 146 147 148 223 224 225 226 300 301 ...
## ..@ coords : num [1:3103, 1:2] 181180 181140 181180 181220 181100 ...
## .. ..- attr(*, "dimnames")=List of 2
## ..@ bbox : num [1:2, 1:2] 178440 329600 181560 333760
## .. ..- attr(*, "dimnames")=List of 2
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1
spplot(meuse.grid, zcol='dist')

Slika 12.20: Metoda spplot za grid.

Slično kao i kod vektorskih podataka mogu se koristiti ili generičke plot funkcije (u ovom slučaju funkcija image) ili spplot funkcije za vizualizaciju rasterskih podataka, Slika 12.20. Prilikom učitavanje rasterskih podataka koristeći rgdal paket kreira se objekat iz klase SpatialGridDataFrame. Primer vizualizacije digitalnog modela visina za deo grada dat je u narednom primeru koristeći image (Slika 12.21) i spplot metod, Slika 12.22.

gr <- readGDAL("grid1.TIF")
image(gr, col = grey(1:99/100),axes=T)

Slika 12.21: Metoda image za grid.

gr <- readGDAL("grid1.TIF")
# generički
image(gr, col = grey(1:99/100),axes=T)
# ili spplot
spplot(gr, main="DEM dela grada", scales = list(draw = T),
col.regions=terrain.colors(20) )

Slika 12.22: Metoda spplot: DEM.

12.5.2 Paket sf

Paket sf (Pebesma 2018, Lovelace, Nowosad, and Muenchow (2018)) predstavlja nove klase za vektorske podatke koje se zasnivaju na OGC otvorenom standardu (eng. simple feature). Reprezentacija geometrijskih entiteta je podržana za: tačke, linije, poligone, više tačaka, više linija, više poligona i geometrijsku kolekciju. Vektorski podaci su organizovani kao data.frame objekti sa proširenjem da se geometrija čuva u obliku WKT ili WKB formata.

Izrada globalne karte očekivanog životnog veka po državama koristeći podatke sadržane u paketu spData (Bivand, Nowosad, and Lovelace 2018) data je u narednom primeru. Metodom st_transform se vrši transformacija iz jednog koordinatnog referentnog sistema u drugi, ovaj metod se koristi ako se manipuliše objektima iz klasa paketa sf. Ovde su podaci transformisani u pseudocilindričnu ekvivalentnu projekciju. Metod plot se odnosi i na sf objekte. Ovim metodom lako se dodaju neki elementi karte, kao na primer koordinatna mreža, koja se dodaje koristeći argument graticule = TRUE, Slika 12.23.

library(sf)
library(spData)
data(world)
world_mollweide = st_transform(world, crs = "+proj=moll")
plot(world_mollweide["lifeExp"], graticule = TRUE, 
     main = "Očekivani životni vek za rođene 2014.")

Slika 12.23: Metoda sf::plot: horoplet karta.

Konverzija meuse podataka iz sp u sf objekat, kao i izrada pregledne karte po svim atributima (Slika 12.24) data je u sledećem primeru:

library(sf)
library(sp)
data(meuse)
coordinates(meuse) = ~x+y
m.sf = st_as_sf(meuse)
str(m.sf)
## $ cadmium : num 11.7 8.6 6.5 2.6 2.8 3 3.2 2.8 2.4 1.6 ...
## $ copper : num 85 81 68 81 48 61 31 29 37 24 ...
## $ lead : num 299 277 199 116 117 137 132 150 133 80 ...
## $ zinc : num 1022 1141 640 257 269 ...
## $ elev : num 7.91 6.98 7.8 7.66 7.48 ...
## $ dist : num 0.00136 0.01222 0.10303 0.19009 0.27709 ...
## $ om : num 13.6 14 13 8 8.7 7.8 9.2 9.5 10.6 6.3 ...
## $ ffreq : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
## $ soil : Factor w/ 3 levels "1","2","3": 1 1 1 2 2 2 2 1 1 2 ...
## $ lime : Factor w/ 2 levels "0","1": 2 2 2 1 1 1 1 1 1 1 ...
## $ landuse : Factor w/ 15 levels "Aa","Ab","Ag",..: 4 4 4 11 4 11 4 2 2 15 ...
## $ dist.m : num 50 30 150 270 380 470 240 120 240 420 ...
## $ geometry:sfc_POINT of length 155; first list element: Classes 'XY', 'POINT', 'sfg' num [1:2] 181072 333611
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "names")= chr "cadmium" "copper" "lead" "zinc"
opar = par(mar=rep(0,4))
plot(m.sf)

Slika 12.24: Metoda sf::plot: multi-panel plot.

12.5.3 Raster

Paket sp podržava rasterske podatke sa dve klase SpatialGridDataFrame i SpatialPixelsDataFrame. Paket raster (R. J. Hijmans 2017 b) sadrži više klasa za rasterske podatke od kojih su tri glavne RasterLayer, RasterBrick i RasterStack. Raster paket sadrži funkcionalnosti za kreiranje, čitanje, pisanje i manipulaciju nad rasterskim podacima, kao što je na primer rasterska algebra. Isto tako paket sadrži i funkcionalnosti za vizualizaciju rasterskih podataka.

12.5.3.1 RasterLayer

Klasa RasterLayer predstavlja rasterski podatak sa jednim atributom. Objekat tipa RasterLayer sadrži pored atributnih podataka i sve parametre koje se tipično čuvaju u zaglavlju (hederu) rasterskog fajla. To su: broj vrsta i kolona, prostorni obuhvat podataka i parametri koordinatnog referentnog sistema. Zavisno od veličine fajla svi podaci mogu biti u RAM memoriji ili im se pristupa sa diska računara, što nije slučaj sa sp klasama kod kojih su podaci uvek u RAM memoriji.

Naredni primer pokazuje kreiranje rasterskog objekta tipa RasterLayer, koji predstavlja rastersku kartu prosečnih godišnjih temperatura u Srbiji. Podaci su preuzeti koristeći Web Coverage Service (WCS) url zahtev. Na kraju su podaci vizualizovani koristeći metod plot koji je deo paketa raster, Slika 12.25.

library(raster)
library(RColorBrewer)
url =paste('http://osgl.grf.bg.ac.rs/geoserver/osgl_2/ows?',
           'service=WCS&version=2.0.1&',
           'request=GetCoverage&coverageId=osgl_2__tac_temp_1961_2010&format=geotiff',sep='')
download.file(url, "temp.tif")
temp= raster("temp.tif")
nlayers(temp)
## [1] 1
# koordinatni ref. sis. (CRS)
crs(temp)
# CRS arguments:
#   +proj=utm +zone=34 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
# Broj piksela
ncell(temp)
## [1] 159390
# broj vrsta kolona i atributa 
dim(temp)
## [1] 483 330   1
# prostorna rezolucija
res(temp)
## [1] 1000 1000
plot(temp, col=brewer.pal(5,"OrRd"), 
     main= "Prosečne godišnje temperature 1961-2010" )

Slika 12.25: Metoda raster::plot godišnje temperature.

12.5.3.2 RasterBrick

RasterBrick je klasa koja podržava više atributa u okviru rasterskih podataka, tipično to važi kad se koristi satelitski snimak sa više kanala, na primer multispektralni snimak iz Sentinel 2 misije ili Landsat 8. Ono što je specifično za RasterBrick je da on može da se referiše samo na jedan fajl koji je sačuvan na disku, na primer GeoTIFF fajl u kome su sačuvani kanali nekog multispektralnog satelitskog snimka.

Demonstracija korišćenja RasterBrick objekta je data nad satelitskim podacima Sentinel 2 misije za deo područja oko Opštine Ruma. Podaci sadrže kanale od 2 do 11 i svi kanali su resamplovani na rezoluciju od 10 m. U Tabeli 12.3 je dat opis svih kanala sa njihovim osnovnim karakteristikama.

Tabela 12.3: Kanali Sentinel 2 misije.
Sentinel-2 Kanali	Talasna dužina (µm)	Rezolucija (m)	Opseg(nm)
Band 1 – Coastal aerosol	0.443	60	20
Band 2 – Blue	0.490	10	65
Band 3 – Green	0.560	10	35
Band 4 – Red	0.665	10	30
Band 5 – Vegetation Red Edge	0.705	20	15
Band 6 – Vegetation Red Edge	0.740	20	15
Band 7 – Vegetation Red Edge	0.783	20	20
Band 8 – NIR	0.842	10	115
Band 8A – Narrow NIR	0.865	20	20
Band 9 – Water vapour	0.945	60	20
Band 10 – SWIR – Cirrus	1.375	60	20
Band 11 – SWIR	1.610	20	90
Band 12 – SWIR	2.190	20	180

Podaci dela satelitskog snimka su učitani koristeći brick metod i kreiran je objekat s tipa RasterBrick, potom su preimenovani atributi objekta s, tako da su imena asocijativna.

s <- brick('S2A_USER_MSI_L2A_TL_SGS__20160427T133230_A004423_T34TDQ_fields.tiff')
names(s)<- c('B02','B03','B04','B05','B06', 'B07', 'B08', 'B08A', 'B11', 'B12')

Objekat s se može ispitati koristeći metode kao što su: summary, nlayers, crs, res.

summary(s)
# B02 B03 B04 B05 B06 B07 B08 B08A B11 B12
# Min. 12 138.0 1 313 758 825 872 959.00 444 1101
# 1st Qu. 207 410.0 197 717 1228 1335 1421 1757.75 1119 1586
# Median 494 643.5 871 1037 1492 1664 1658 2517.50 2291 1919
# 3rd Qu. 568 737.0 1007 1171 2484 3366 3479 2738.00 2547 3547
# Max. 1395 1771.0 2141 2280 4379 6270 6349 3461.00 3199 6486
# NA's 0 0.0 0 0 0 0 0 0.00 0 0
nlayers(s)
# [1] 10
crs(s)
# CRS arguments:
# +proj=utm +zone=34 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
ncell(s)
# [1] 90000
dim(s)
# [1] 300 300 10
res(s)
# [1] 10 10

Koristeći metod plotRGB može se kreirati karta koristeći neku od kolor kompozicija, ne nužno crveni, zeleni i plavi kanal, Slika 12.26.

plotRGB(s, r = 3, g = 2, b = 1, axes = TRUE, 
        stretch = "lin", main = "Senitnel R-G-B 27-04-2016")

Slika 12.26: Metoda raster::plot Senitnel R-G-B 27-04-2016.

Slično se može kreirati kompozicija kanala: bliskoinfracrveni, crveni i zeleni, Slika 12.27.

plotRGB(s, r = 7, g = 3, b = 2, axes = TRUE, 
        stretch = "lin", main = "Sentinel NIR-R-G 27-04-2016")

Slika 12.27: Metoda raster::plot Sentinel NIR-R-G 27-04-2016.

U sledećem primeru napravljena je funkcija VI koja generiše vegetacioni indeks. Funkcija ima 3 argumenta img koji je tipa RasterBrick i argumenata i i k koji predstavljaju kanale koji će biti korišćeni. Potom je primenjena funkcija VI i kreiran rasterski podatak ndvi (normalizovani vegetacioni indeks), a potom je ndvi i kartiran koristeći plot metod, , Slika 12.28.

# i,k redni broj atributa u objektu s
VI <- function(img, i, k) {
  bi <- img[[i]]
  bk <- img[[k]]
  vi <- (bk - bi) / (bk + bi)
  return(vi)
}
# Sentinel NIR = B08  -  s[[7]] , red = B04  - 3.
ndvi <- VI(s, 3, 7)
plot(ndvi, col = rev(terrain.colors(30)), 
     main = 'NDVI - Sentinel 27-04-2016')

Slika 12.28: Metoda raster::plot Sentinel NDVI.

12.5.3.3 RasterStack

RasterStack je lista rasterskih lejera, vrlo slične strukture kao RasterBrick, ali ne postoji ograničinje da se referiše na jedan GIS fajl, već može da bude kompozicija više fajlova koji su sačuvani na disku ili koji su u RAM memoriji računara, ili kombinacija istih. Praktično se može reći da je to lista objekata tipa RasterLayer. Podrazumeva se da je prostorni obuhvat i rezolucija ista za svaki RasterLayer koji je deo RasterStack objekta.

12.6 Prikaz reljefa

Koristeći paket raster i rasterVis (Perpinan Lamigueiro and Hijmans 2018) mogu se kreirati različite vizualizacije reljefa, naredni primer kreira perspektivni prikaz reljefa, koji koristi predefinisanu paletu boja terrain.colors (Slika 12.29), a umesto koje bi i korisnik sam mogao da definiše bilo koju paletu ili da je napravi sam.

library(RColorBrewer)
library(rgl)
require(raster)
library(rasterVis)
alt = getData('alt', country='SRB') # 1km rezolucija
plot3D(alt, drape=NULL, zfac=0.5,
       col=colorRampPalette(terrain.colors(12) ))

Slika 12.29: Metoda plot3D perspektivni prikaz reljefa.

Prikaz reljefa koristeći metodu senke i hipsometrijsku skalu (metodu boja) može se napraviti na sledeći način. Karta je plotovana direktno u jpg fajl, Slika 12.30.

slope = terrain(alt, opt='slope') # nagib
aspect = terrain(alt, opt='aspect') # orjentacija
hill = hillShade(slope, aspect, 45, 270) # senka deklinacija 45 , azimut 270
# plot u jpeg
jpeg(file = "dem.jpeg", width =2000, height = 2000, res = 300, pointsize = 14)
plot(hill, col=grey(0:100/100), legend=FALSE, 
     main='Srbija prikaz reljefa - senke i hipsometrija')
plot(alt, col=terrain.colors(12, alpha=0.35), add=TRUE) # metoda boja
dev.off()

Slika 12.30: Srbija prikaz reljefa - senke i hipsometrija.

Zanimljiva opcija vizualizacije reljefa koja daje marginalne profile koji uzimaju prosek visina u pravcu zapad-istok i u pravcu sever-jug data je u narednom primeru, Slika 12.31.

p0 <- levelplot(hill, margin=FALSE, par.settings = GrTheme)
p1 = levelplot(alt, layers = 1, colorkey=FALSE,
  col.regions= terrain.colors(12,alpha=0.35) ,
  margin = list(FUN = 'mean'), contour=F)
p1 + as.layer(p0, under = TRUE)

Slika 12.31: Srbija prikaz reljefa sa prosečnim profilima.

Metoda izohipsi i metoda boja za područje u okolini Zlatibora je ilustrovana sledećim primerom, Slika 12.32.

zlatibor = crop(alt, extent (19.65,19.75, 43.60, 43.75))
plot(zlatibor, col= terrain.colors(160,alpha=0.5) )
contour(zlatibor, add=T)

Slika 12.32: Srbija prikaz reljefa - izohipse.

12.7 Izrada tematskih karata u paketu ggplot2

Paket ggplot2 (Wickham and Chang 2016) uvodi nov sistem za izradu grafika i karata u odnosu na tradicionalne metode u R-u. Paket ggplot2 je razvio i održava ga Hadley Wickham. Paket je nastao na bazi knjige “The Grammar of Graphics” čiji je autor Leland Wilkinson (Wilkinson 2006). Izrada grafika i karata koristeći ggplot sintaksu je takva da se koristi lejerska struktura u izradi vizualizacije.

Lejerski pristup može biti dekoponovati na sledeće delove:

Podaci koji se plotuju ili kartiraju,
Stilizacija (estetika),
Geometrija,
Statistika (ako se podaci sumiraju ili grupišu i sl.).

U sledećem primeru je prikazano kako se kreira karta koristeći proporcionalne simbole nad meuse podacima (Slika 12.33), veličina simbola je proporcionalna atributnom podatku zinc (koncentracija cinka u ppm). Metoda ggplot očekuje da i prostorni podaci budu konvertovani u data.frame objekat.

library(sp)
data(meuse)
coordinates(meuse) = \~x+y
library(ggplot2)
ggplot(as.data.frame(meuse)) +
aes(x, y, size = zinc) +
geom_point() +
coord_equal() + labs(x = "E", y="N")

Slika 12.33: Metod ggplot: proporcionalni simboli.

Izrada horoplet karte koristeći ggplot je data u sledećem primeru. Gemetrija je iz sp objekta konvertovana u data.frame koristeći metod fortify, koji je kreira tabelu sa atributima, a koordinate prelomnih tačaka su grupisane prema atributu ID_Municip (id opštine). Potom su atributi iz tabele koja sadrži geometriju spojeni sa atributnim podacima koristeći funkciju merge, gde je definisan i zajednički atribut. Koristeći metod ggplot je napravljena horoplet karta broja stanovnika po opštinama u Boneovoj (Bonne) pseudocilindričnoj ekvivalentnoj projekciji, Slika 12.34. Paket ggplot koristi notaciju R paketa mapproj za zadavanje kartografskih projekcija i računanje koordinata u projekciji.

library(rgdal)
url=paste('http://osgl.grf.bg.ac.rs/geoserver/osgl_4/wfs?',
          'service=WFS&version=1.0.0&request=GetFeature&typeName=osgl_4:POP_age',
          sep="")
download.file(url, 'poage.gml')
pa <- readOGR(dsn ='poage.gml', layer = 'POP_age', stringsAsFactors=FALSE )
proj4string(pa) <- CRS("+init=epsg:4326")
pa$TOTAL_POP <-as.numeric(pa$TOTAL_POP)
f = fortify(pa, region = "ID_Municip")
ff = merge(f, pa@data, by.x = "id", by.y = "ID_Municip")
library(mapproj)
ggplot(ff) +
aes(x=long, y=lat, group = group, fill = TOTAL_POP) +
geom_polygon() +
coord_map("bonne", lat0 = 44)

Slika 12.34: Metod ggplot: horoplet karta.

12.8 Izrada tematskih karata koristeći paket tmap

Paket tmap (Tennekes 2018) je napravljen za izradu tematskih karata koristeći isti pristup kao i ggplot2 paket. Pored lejerske organizacije za izradu karata može se koristiti i standardan pristup za brzu izradu tematskih karata, metoda qtm. Za razliku od ggplot2, tmap koristi prostorne klase (sp, sf i raster klase) kao ulazne podatke, namenjen je isključivo izradi karata, i nudi veliku fleksibilnost kod izrade elemenata karte.

Lejerska struktura kod paketa tmap je sledeća:

Prostorni podaci;
- Tačke, linije, poligon ili raster,
- Prostorni okvir,
- Koordinatni referentni sistem.
Lejeri (kartografsko oblikovanje);
- Estetika,
- Statistika,
- Razmera.
Templejt karte.

Sledeći primer prikazuje izradu karte broja stanovnika po opštinama koristeći qtm metod, karta je u Lambertovoj azimutnoj ekvivalentnoj projekciji, Slika 12.35.

library(tmap)
library(tmaptools)
# transformacija u Lambertovu azimutnu ekvivalentnu proj
pa <- spTransform(pa, 
        CRS('+proj=laea +lat_0=44 +lon_0=21 +x_0=200000 +y_0=200000 +datum=WGS84'))
# izrade tematske horoplet karte
qtm(pa, fill="TOTAL_POP", style="gray", 
    fill.title="Broj stanovnika", fill.textNA="Nema podataka", 
    legend.format.text.separator="do")

Kartografsko oblikovanje u primeru izrade “brze tematske karte” definisano je u okviru jedne funkcije kao kod tradicionalnih plot metoda. Zadata je pozadina sive boje, vizualna promenljiva varira u zavisnosti od ukupnog broja stanovnika po poligonu opštine, pored naslova zadat je i separator za numeričke podatke u legendi, kao i kako se definiše oznaka za entitete koji imaju nedostajuće podatke.

Slika 12.35: Metod qtm: horoplet karta.

Sa istim podacima može se napraviti horoplet karta koristeći lejersku organizaciju za izradu tematske karte. Detaljno je prikazan postupak izrade uključujući: klasifikaciju atributa, iscrtavanje koordinatne mreže, orjentacije, definisanja legende i njene pozicije; i na kraju teme koja sadrži templejt za okvirnu i vanokvirnu stilizaciju, Slika 12.36.

Paket tmap nudi još mnogo mogućnosti za izradu tematskih karti, pored ostalog i interaktivnih tematskih web karata u kombinaciji sa OpenStreetMaps sadržajem koristeći Leaflet JavaScript biblioteku u klijentskom delu, takav primer je sadržan u vinjeti paketa.

tm_shape(pa ) +
tm_polygons(c("TOTAL_POP"),
textNA="Nema podataka",
style=c("pretty"),
palette="YlOrRd",
auto.palette.mapping=FALSE,
title="Broj stanovnika") +
tm_grid(projection="longlat", labels.size = .5) +
tm_compass(type= "arrow", color.light = "grey90") +
tm_style_classic(inner.margins=c(.04,.06, .04, .04), 
                 legend.position = c("right", "top"),
                 legend.frame = TRUE, bg.color="lightblue", 
                 legend.bg.color="lightblue",
                 earth.boundary = TRUE, space.color="grey90",
                 legend.format=list(text.separator="do",
                    format=c(format="f", big.mark = " ", digits=0)
))

Slika 12.36: Horoplet karta, paket tmap.

12.9 Izrada tematskih karata koristeći paket cartography

Paket cartogrphy (Giraud and Lambert 2017) koristi vrlo sličan pristup kao tmap paket.

library(cartography)
# Konvertovanje iz sp u sf tip
pasf <- st_as_sf(pa)
# Udeo odraslih u ukupnoj populaciji
pasf$ADULTS <- as.numeric(pasf$ADULTS)
pasf$deoStarijih <- 100 * pasf$ADULTS/pasf$TOTAL_POP
# Selektovanje svih poligona koji u imenu sadrže reč Beograd
library(stringr)
pasf <- pasf[str_detect(pasf$NameASCII, "Beograd"), ]
# Izbor palete boja
cols <- carto.pal(pal1 = "wine.pal", n1 = 5)
# Karta će biti štampana u jpeg fajlu
jpeg(file = "cart.jpeg", width =2000, height = 2000, res = 300, pointsize = 14)
opar <- par(mar = c(0,0,1.2,0))
# preuzimanje OSM tajlova
BGosm <- getTiles(pasf, type = "osmgrayscale", crop = TRUE)
########### TEMATSKA KARTA ########################
# Iscrtavanje osnovne karte
tilesLayer(BGosm)
# Iscrtavanje poligona, samo granice
plot(st_geometry(pasf), col =NA, border = "black",
bg = NA, lwd =3, add = TRUE)
# Proporcionalni simboli, veličina proporcionalna broju stanovnika, boja je su
# raspodeljene u odnosu na udeo odraslih u ukupnoj populaciji.
propSymbolsChoroLayer(x = pasf, # sf objekat
var = "TOTAL_POP", # atribut za proporcionalne simbole
var2 = "deoStarijih", # atr. za izbor boja
col = cols, # paleta boja
inches = 0.2, # radijus najvećeg simbola
method = "quantile", # klasifikacija atributa
border = "grey50", # boja granične linije prop. sim.
lwd = 1, # debljina granične linije prop. sim.
legend.var.pos = "topleft", # pozicija prve legende
legend.var2.pos = "left", # pozicija druge legende
legend.var2.title.txt ="Procenat \n odraslih \n stanovnika",
legend.var.title.txt = "Broj stanovnika",
legend.var.style = "c") # stil legende
# Dodavanje layouta
layoutLayer(title="Udeo odraslih u Beogradu, popis 2011", # Naslov karte
scale = 20, # veličina razmernika
north = TRUE, # pravac severa
col = "white",
coltitle = "black",
author = "Milan Kilibarda",
sources = "OSM, RSZ Popis, 2011, OSGL",
frame = TRUE)
# predhodno podešavanja grafičkih parametara
par(opar)
# zatvaranje aktuelnog plot prozora
dev.off()

Predhodni primer kombinuje granice opština, pozadinsku kartu koja je preuzeta sa OpenStreetMap servisa, proporcionalne simbole čija veličina odgovara ukupnom broju stanovnika po opštini, a boja proporcionalnog simbola varira u zavisnosti od udela starijih stanovnika u ukupnom broju stanovnika. U narednom primeru, objašnjenje koda je dato kroz komentare, Slika 12.37.

Slika 12.37: Vizualizacija sa više vizualnih promenljivih.

Očigledno je da u starom delu Beograda ima nekoliko opština koje su izrazito male po površini. Alternativni pristup u tematskoj kartografiji je da naprvimo heksagone istih površina, koji će dati naglasak samo na kartirane promenljive. Dakle, površine će biti zanemarene jer je svaki heksagon iste površine. Kod ovakvih karata nastoji se da prostorna orjentacija bude što sličnija izvornim podacima, Slika 12.38. Sledeći kod kreira heksagon kartu, a objašnjenje koda je dato u komentarima.

# devtools::install_github("jbaileyh/geogrid")
# paket sa funkcionalnostima da poligone kovertuje u heksagone
library(geogrid)
library(sf)
# kreiranje heksagona za opštine u Beogradu
hex_grid <- calculate_grid(shape = as(pasf, 'Spatial'), grid_type = "hexagonal", seed = 1)
# kreiranje sf poligona sa atributima iz obj. pasf sa heksagon geometrijom hex_grid obj. 
pol = st_as_sf  ( assign_polygons(as(pasf, 'Spatial'), hex_grid ) )
# kreiranje fajla u kom će biti smeštena karta
jpeg(file = "cart2.jpeg", width =2000, height = 2000, res = 300, pointsize = 14)
# granice heksagona
plot(st_geometry(pol), col =NA, border   = "black", 
     bg = NA, lwd =3 )
# proporcionalni simboli
propSymbolsChoroLayer(x =pol, # sf obj.
                      var = "TOTAL_POP", # atr. za veličinu simbola
                      var2 = "deoStarijih", #  atr. za izbor boja
                      col = cols, # boje
                      inches = 0.2, # maksimalni radijus prop. sim.
                      method = "quantile", # intervali 
                      border = "grey50", # boja granične linije za prop. sim.
                      lwd = 1, # debljina gran. lin. prop. sim.
                      legend.var.pos = "topleft", # pozicija prve legende
                      legend.var2.pos = "left", # pozicija druge legende
                      legend.var2.title.txt =  
                        "Procenat \n odraslih \n stanovnika", 
                      legend.var.title.txt = "Broj stanovnika",
                      legend.var.style = "c") # stil legende
# skraćivanje imena opš. u BGD za karakter "Beograd-"
pol$Naziv  = gsub("Beograd-", "", pol$NameASCII)
# kreiranje toponima
labelLayer(spdf = as(pol,'Spatial'), # SpatialPolygonsDataFrame za toponime
           txt = "Naziv", # atribut koji sadrži toponime
           col = "black", #  boja toponima
           cex = 0.5, # veličina za toponime
           font = 1) # font
# Ostali sadržaj
layoutLayer(title="Udeo odraslih u Beogradu, popis 2011", # Naslov
            scale = 20, # veličina razmernika u km
            north = TRUE, # Orjentacija
            col = "white",
            coltitle = "black",
            author = "Milan Kilibarda",  
            sources = "RSZ Popis, 2011, OSGL Beograd",
            frame = TRUE)  
dev.off()

Slika 12.38: Karta sa heksagon geometrijom.

12.10 Izrada interaktivnih tematskih karata koristeći paket plotly

R paket plotly je namenjen izradi interaktivnih grafika na bazi JavaScript grafičke biblioteke plotly.js. Veoma obimna dokumentacija za R paket je dostupna na sledećim adresama:

U sledećem primeru je kreirana karta opština u Beogradu, gde je svaka opština prikazana različitom bojom, a dobija se karta koja je interaktivna u okviru RStudio okruženja, Slika 12.39.

library(plotly)
cols <- c("#89C5DA", "#DA5724", "#74D944", "#CE50CA", "#3F4921",
          "#C0717C", "#CBD588", "#5F7FC7", "#673770", "#D3D93E",
          "#38333E", "#508578", "#D7C1B1", "#689030", "#AD6F3B",
          "#CD9BCD", "#D14285")
plot_ly(pasf, color = ~NameASCII, colors = cols, alpha = 1)

Slika 12.39: Paket plotly: interaktivna karta.

U narednom primeru se koristi “pipe” sistem sa operatorom %>%, kojim se skraćuje kod, jer se međurezultati direktno prosleđuju, a da predhodno nisu sačuvani u nekoj promenljivoj. Ovaj sistem je prvi put uveden sa paketom magrittr, a sada se koristi uglavnom kao podrazumevani paket u sklopu mnogih drugih kao što je plotly, tidyverse i drugih.

Jednostavan primer pokazuje najprostiju primenu ovog sistema.

 library(magrittr)
  x<- 60 *pi/180
  # prvi način sa ugnježdenim funkcijama
  round(exp(sin(tan(x))), 1)
  ## [1] 2.7
  # drugi način, korak po korak sa pomoćnim prom.
  x1 <- tan(x) 
  x2 <- sin(x1)
  x3 <- exp(x2)
  round(x3,1)
  ## [1] 2.7
  # pipe način x u tan onda sinus onda exp onda zaokruži
  x %>% tan() %>% sin() %>% exp() %>% round(1)  
  ## [1] 2.7

Dakle, svaki od tri načina računanja datog izraza daje iste rezultate, samo je pitanje pristupa u čitljivosti koda koji princip će se koristiti. Detaljan opis pipe koncepta dostupan je knjizi “R for Data Science” (Wickham and Grolemund 2016).

Paket nudi mogućnost za izradu izradu interaktivnih web karata baziranih na mapbox servisima. Sledeći primer kreira interaktivnu kartu opština u Beogradu, sa satelitskim snimcima koje obezbeđuje mapbox servis u pozadini (Slika 12.40), a korišćen je i pipe operator. Da bi primer radio neophodno je obezbediti mapbox ključ, koji je za ovaj primer snimljen kao sistemska promenljiva u R okruženju.

Sys.setenv('MAPBOX_TOKEN' = 'tajni ključ')
pasf <-  st_transform(pasf, crs = 4326)
# izrada karte
plot_mapbox(pasf, fillcolor = 'rgba(7, 164, 181, 0.2)' ) %>%
    layout(
      title = "Opštine u Beogradu",
      mapbox = list(style = "satellite-streets")
    )

Slika 12.40: Paket plotly: interaktivna web karta.

Sledeći primer pokazuje kreiranje interaktivne karte i tabele sa atributima, Slika 12.41. Interakcija sa kartom koja je u levom delu sinhronizovana je sa prikazom iz tabele. Tako da ako je selektovana jedna opština na karti onda su u tabeli prikazani atributi samo te opštine, selekcija radi i za selekciju više opština istovremeno u tabeli ili na karti.

library(crosstalk)
  pasf <-  st_transform(pasf, crs = 4326)
  bgd <- SharedData$new(pasf)
  # jedino radi sa WGS84
  bscols(
    plot_mapbox(bgd, text = ~NameASCII, hoverinfo = "text"),
    DT::datatable(bgd)
  )

Slika 12.41: Interaktivna karta i tabela sa atributima.

Sledeći primer je sličan prethodnom i radi sa sinhronizacijom histograma i karte. Nudi i mogućnosti više selekcija koje se prikazuju različitim bojama i na karti i na histogramu, Slika 12.42.

bscols(
    plot_mapbox(bgd) %>%
    highlight(dynamic = TRUE, persistent = T),
    plot_ly(bgd, x = ~TOTAL_POP) %>% 
    add_histogram(xbins = list(start = 20000, 
                    end = 215000, size = 30000)) %>%
    layout(barmode = "overlay") %>% 
    highlight("plotly_selected", persistent = T)
  )

Slika 12.42: Interaktivna karta i histogram.

12.11 Izrada web karata koristeći paket plotGoogleMaps

Jedan od prvih R paketa za izradu interaktivnih web karata je plotGoogleMaps (Kilibarda 2017) koji je nastao 2010 godine, paket je baziran na Google Maps API JavaScript biblioteci. Osnovne karte se mogu birati koristeći globalne karte kompanije Google. Paket je baziran na standardnom plot pristupu i radi sa sp i raster klasama.

Sledeći primer kreira web kartu nad meuse podacima (Slika 12.43), karta je interaktivna i klikom na entitet dobijaju se atributni podaci, a može se birati osnovna karta i slično.

Slika 12.43: Interaktivna karta meuse podataka, plotGoogleMaps

Jednostavno je kreirati kartu koja koristi proporcionalne simbole za vizualizaciju prostornih fenomena, Slika 12.44.

m2 <- bubbleGoogleMaps(meuse, 'zinc', filename = 'm.html')

Slika 12.44: Interaktivna karta meuse podataka, bubbleGoogleMaps

R sadrži odličan paket spacetime (Pebesma 2016) koji obezbeđuje klase i metode za prostorno vremenske podatke. Više o paketu videti u vinjeti paketa (Pebesma 2012, Bivand, Pebesma, and Gomez-Rubio (2013)). Naredni primer prikazuje jedan od klasičnih prikaza prostorno vremenskih podataka, a to je panel prikaz. Prostorno vremenski podaci prikazani kao niz prostornih podataka, gde svaki element odgovara jednom momentu u vremenu, Slika 12.45. Alternativa je animirani prikaz ili kombinacija karte i grafika. Sledeći kod dat je bez pojašnjenja, a na osnovu prostorno vremenskih podataka o broju nezaposlenih u SAD kroz vreme kreira objekat koji sadrži prostorno vremenske informacije.

# Data preparation
  # STFDF data from spacetime vignette spacetime: Spatio-Temporal Data in R
  library("maps")
  states.m = map('state', plot=FALSE, fill=TRUE)
  IDs <- sapply(strsplit(states.m$names, ":"), function(x) x[1])
  library("maptools")
  states = map2SpatialPolygons(states.m, IDs=IDs)
  yrs = 1970:1986
  time = as.POSIXct(paste(yrs, "-01-01", sep=""), tz = "GMT")
  library("plm")
  data("Produc")
# no data for Colorado
  Produc.st = STFDF(states[-8], time, Produc[order(Produc[2], Produc[1]),])
  Produc.st@sp@proj4string=CRS('+proj=longlat +datum=WGS84')

Slika 12.45: Prostorno vremenski prikaz nezaposlenosti u SAD, multi-panel plot

Za ovakvu strukturu podataka u plotGoogleMaps paketu moguće je napraviti prostorno vremenski prikaz u kome je karta kreirana po prosečnoj vrednosti nezapolsenosti u celom vremenskom intervalu (metoda boja za vizualnu promenljivu), a nakon klika na entitet dobija se grafik koji prikazuje vremensku seriju kartiranog atributa, Slika 12.46.

# Procenat nezaposlenih po državama u SAD od 1970-1986
m <- stfdfGoogleMaps(Produc.st, zcol= 'unemp', filename = 'mst.html')

Slika 12.46: Prostorno vremenski prikaz nezaposlenosti u SAD, stfdfGoogleMaps

12.12 Izrada web karata koristeći paket mapview

Paket mapview (Appelhans et al. 2018) je jedan od niza paketa koji se razvija u R u poslednjih nekoliko godina koji je namenjen izradi web karata, koje su interaktivne i bazirane na Leafleat JavaScript bibloteci. Paket je jednostavan za korićenje i ima klasičan plot pristup. Paket nudi odlične funkcionalnosti za sve standardne tipove prostornih podataka u R-u.

Sledeći primer prikazuje izradu Web karte nad podacima meuse, ako drugačije nije podešeno OpenStreetMap je osnovna karta. U ovom primeru ulazni podaci su sp objekat. Karta je interaktvina, mogu se dobitit atributni podaci za svaki entitet. Takođe, može se menjati osnovna karta i drugo.

library(mapview)
## Loading required package: leaflet
library(sp)
data(meuse)
coordinates(meuse) = ~x+y
proj4string(meuse) <- CRS("+init=epsg:28992")
mapview(meuse, zcol = "zinc", legend = TRUE)

Slika 12.47: Interaktivna karta meuse podataka, mapview

Sedeći primer kreira interaktivnu horoplet kartu po opštinama, a kartira se atribut koji predstavlja broj stanovnika po opštinama, podaci su korišćeni u ranijim primerima, a objekat je tipa sf. Za ovu kartu podešena je providnost na 20%.

mapview(pasf, zcol = c("TOTAL_POP"), alpha.regions = 0.2, legend = TRUE)

Slika 12.48: Interaktivna horoplet karta, mapview

Sledeći primer prikazuje vizualizaciju prostorno vremenskih podataka, koji predstavljaju NDVI indeks za područje planine Kilimandžaro. Napravljena je web karta koja ima četiri panela, koji su sinhronizovani, tako da je ovaj prikaz vrlo pogodan za vizualnu analizu promena u vremenu.

library(raster)
kili_data <- system.file("extdata", "kiliNDVI.tif", package = "mapview")
kiliNDVI <- stack(kili_data) # 23 lejera u seriji
library(RColorBrewer)
clr <- colorRampPalette(brewer.pal(9, "YlGn")) # paleta boja
m1 <- mapview(kiliNDVI[[1]] , col.regions = clr, 
              legend = TRUE, map.types = "Esri.WorldImagery")
m2 <- mapview(kiliNDVI[[5]], col.regions = clr ,
              legend = TRUE, map.types = "Esri.WorldImagery")
m3 <- mapview(kiliNDVI[[10]], col.regions = clr, 
              legend = TRUE, map.types = "Esri.WorldImagery")
m4 <- mapview(kiliNDVI[[20]],col.regions = clr, 
              legend = TRUE, map.types = "Esri.WorldImagery")
sync(m1, m2, m3, m4) # 4 panela sinhronizovana

Slika 12.49: Sinhronizoavne NDVI karte, mapview

Paket mapview nudi jednostavan prikaz prostorno vremenskih fenomena u 3D prostoru, gde je Z vremenska koordinata. Primer je dat nad istim podacima kao za sinhronizovani prikaz.

cubeView(kiliNDVI[[c(1,5,10,20)]], col.regions = clr)

Slika 12.50: Prikaz prostorno vremenskih fenomena u 3D prostoru

12.13 Paket plotKML za vizualizaciju na virtuelnim globusima

Paket plotKML (Hengl 2017) je paket koji obezbeđuje izradu tematskih karata na virtuelnim globusima, prvenstveno na Google Earth virtuelnom globusu, ali i na drugim kao što je CesiumJS. Praktično sva vizualizacija je smeštena u samom KML fajlu.

Sledeći primer kreira proporcionalne simbole za meuse podatke, gde veličina simbola varira po atributu cink.

library(plotKML)
shape = "http://maps.google.com/mapfiles/kml/pal2/icon18.png"
kml(meuse, colour = zinc, size = zinc,
alpha = 0.75, file = "meuse.kml", shape=shape, labels =zinc)

Slika 12.51: Proporcionalni simbole za meuse podatke na virtuelnom globusu.

Mogućnosti trodimenzionalnog prikaza omogućavaju da se kreira interaktivna karta gde je varijacija prikazana po visini simbola u odnosu na koncentraciju cinka, kao i u predhodnom simbolu kartografska promenljiva varira i u boji.

kml(meuse, colour = zinc,
alpha = 0.75, file = "meuseH.kml", shape=shape, altitude =zinc, labels=zinc)

Slika 12.52: Proporcionalni simboli za meuse podatke na virtuelnom globusu, kartografska preomenljiva varira po visini.

Ako vizualizujemo opštine u Beogradu tako da su prikazane kao prizme gde je visina prizme proporcionalna broju stanovnika, onda takav fajl možemo otvoriti u nekom drugom virtuelnom globusu, na primer CesiumJS. CesiumJS nudi veliki izbor osnovnih karata koje su prikazane na globusu.

library(stringr)
BG <- pa[str_detect(pa$NameASCII, "Beograd"), ]
kml(BG, colour = TOTAL_POP,
alpha = 0.75, file = "BGPop.kml", altitude =TOTAL_POP, labels=TOTAL_POP)

Slika 12.53: Horoplet karta sa varijacijama po visini.

Sledeći primer kreira animirani prikaz prostorno vremenskih podataka za meteorološku mrežu stanica u Hrvatskoj, prikazani atribut je temperatura. Google Earth nudi kontrolu animiranog prikaza na način da je moguće:

promena unapred,
promena unazad,
brzo,
sporo,
pauza,
selekcija vremenskog opsega.

Slično nudi i CesiumJS, osim selekcije vremenskog intervala.

library(plotKML)
  data(HRtemp08)
  # format the time column:
  HRtemp08$ctime <- as.POSIXct(HRtemp08$DATE, format="%Y-%m-%dT%H:%M:%SZ")
  # create a STIDF object:
  library(spacetime)
  sp <- SpatialPoints(HRtemp08[,c("Lon","Lat")])
  proj4string(sp) <- CRS("+proj=longlat +datum=WGS84")
  HRtemp08.st <- STIDF(sp, time = HRtemp08$ctime, data = HRtemp08[,c("NAME","TEMP")])
  # write to a KML file:
  HRtemp08_jan <- HRtemp08.st[1:500]
  shape <- "http://maps.google.com/mapfiles/kml/pal2/icon18.png"
  kml(HRtemp08_jan, dtime = 24*3600, colour = TEMP, shape = shape, labels = "", kmz=TRUE)

Slika 12.54: Animirani prikaz prostorno vremenskih podataka.

Literatura

R Core Team. 2018. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Pebesma, Edzer, and Roger Bivand. 2018. Sp: Classes and Methods for Spatial Data. https://CRAN.R-project.org/package=sp.

Hijmans, Robert J. 2017a. Geosphere: Spherical Trigonometry. https://CRAN.R-project.org/package=geosphere.

Bivand, Roger, Tim Keitt, and Barry Rowlingson. 2018. Rgdal: Bindings for the ’Geospatial’ Data Abstraction Library. https://CRAN.R-project.org/package=rgdal.

Bivand, Roger, and Colin Rundel. 2017. Rgeos: Interface to Geometry Engine - Open Source (’Geos’). https://CRAN.R-project.org/package=rgeos.

Bivand, Roger. 2017. ClassInt: Choose Univariate Class Intervals. https://CRAN.R-project.org/package=classInt.

Neuwirth, Erich. 2014. RColorBrewer: ColorBrewer Palettes. https://CRAN.R-project.org/package=RColorBrewer.

Brewer, Cynthia A., Geoffrey W. Hatchard, and Mark A. Harrower. 2003. “ColorBrewer in Print: A Catalog of Color Schemes for Maps.” Cartography and Geographic Information Science 30 (1). Taylor & Francis: 5–32. doi:10.1559/152304003100010929.

Sumner, Michael D. 2016. Graticule: Meridional and Parallel Lines for Maps. https://CRAN.R-project.org/package=graticule.

Pebesma, Edzer. 2018. Sf: Simple Features for R. https://github.com/r-spatial/sf/.

Lovelace, Robin, Jakub Nowosad, and Jannes Muenchow. 2018. “Geocomputation with R.”

Bivand, Roger, Jakub Nowosad, and Robin Lovelace. 2018. SpData: Datasets for Spatial Analysis. https://CRAN.R-project.org/package=spData.

Hijmans, Robert J. 2017b. Raster: Geographic Data Analysis and Modeling. https://CRAN.R-project.org/package=raster.

Perpinan Lamigueiro, Oscar, and Robert Hijmans. 2018. RasterVis: Visualization Methods for Raster Data. https://CRAN.R-project.org/package=rasterVis.

Wickham, Hadley, and Winston Chang. 2016. Ggplot2: Create Elegant Data Visualisations Using the Grammar of Graphics. https://CRAN.R-project.org/package=ggplot2.

Wilkinson, Leland. 2006. The Grammar of Graphics. Springer Science & Business Media.

Tennekes, Martijn. 2018. Tmap: Thematic Maps. https://CRAN.R-project.org/package=tmap.

Giraud, Timothée, and Nicolas Lambert. 2017. Cartography: Thematic Cartography. https://CRAN.R-project.org/package=cartography.

Wickham, Hadley, and Garrett Grolemund. 2016. R for Data Science: Import, Tidy, Transform, Visualize, and Model Data. “ O’Reilly Media, Inc.”

Kilibarda, Milan. 2017. PlotGoogleMaps: Plot Spatial or Spatio-Temporal Data over Google Maps. https://R-Forge.R-project.org/projects/

Pebesma, Edzer. 2016. Spacetime: Classes and Methods for Spatio-Temporal Data. https://CRAN.R-project.org/package=spacetime.

Pebesma, Edzer. 2012. “spacetime: Spatio-Temporal Data in R.” Journal of Statistical Software 51 (7): 1–30. http://www.jstatsoft.org/v51/i07/.

Bivand, Roger S., Edzer Pebesma, and Virgilio Gomez-Rubio. 2013. Applied Spatial Data Analysis with R, Second Edition. Springer, NY. http://www.asdar-book.org/.

Appelhans, Tim, Florian Detsch, Christoph Reudenbach, and Stefan Woellauer. 2018. Mapview: Interactive Viewing of Spatial Data in R. https://CRAN.R-project.org/package=mapview.

Hengl, Tomislav. 2017. PlotKML: Visualization of Spatial and Spatio- Temporal Objects in Google Earth. https://CRAN.R-project.org/package=plotKML.