Práce s vektorovými daty v R

Řada předchozích příkladů pro práci s multivariačními daty nám umožňuje výrazně zmenšit prostor (počet ukazatelů) a v ideálním případě pak aplikovat na výsledky těchto analýz prostorové analýzy monovariačních dat. Existují sice také monovariační verze, např. pro prostorovou autokorelace baliček ade4 a funkce multispati(), ale v tomto posledním cvičení se zaměříme na práci s vektorovými daty v prostředí R. Ukážeme si základy práce – načtení dat, možnosti vizualizace a některé základní operace s prostorovými daty. Již jsme se v úvodu setkali s vybranými balíčky pro práci s grafy, ale nabídka balíčků je opravdu široká a níže jsou uvedeny některé důležité pro práci s prostorovými daty. Na stránkách R projektu je uvedeno také celé téma (task) zaměřeno na prostorové analýzy dat, jsou zde podrobně rozepsány podoblasti prostorových analýz a v rámci každé podoblasti jsou vypsané do nich spadající balíčky

názvy balíčků	účel
ggplot2, ggmap, animation, mapplots	grafické výstupy
sp, rgdal, rgeos, maptools	práce s prostorovými daty
spatstat, spdep, splancs	prostorové analýzy
plotKML	export do prostředí Google Earth
RColorBrewer	práce s barvami

• Třídy pro prostorová data: Protože vstupem nebo výstupem analýz může být jiný, byť podobný, objekt, balíky spadající do této kategorie zabezpečují konverze mezi jednotlivými třídami prostorových dat. Základním balíkem definujícím prostorová data je balík sp.
• Manipulace s prostorovými daty: balíky s volně šiřitelnými geodatasety, pro výpočty vzdáleností a ploch, atd.
• Čtení a zápis prostorových dat: Vstup a výstup rastrových i vektorových formátů, konverze mezi nimi, spolupráce s webovými poskytovateli dat (Google, Open Street Maps – OSM), transformace souřadnic, atd. Základním balíkem této podoblasti je rgdal.
• Analýzy prostorových vzorů: hlavním balíčkem této části je spatial nebo spatstat, který má velmi kvalitní dokumentaci.
• Geostatistika: základní operace obsahuje balíček gstat a dále pak např. geoR a geoRglm
• Analýzy areálových dat: v této oblasti je důležitým balíčkem spdep, který umožňuje definovat prostorové sousedství, konstrukci vah a dále pak metody jako je prostorová autokorelace.
• Prostorová regrese: pro bodová data je vhodný balíček nlme a pro areálová data opět balíček spdep
• Ekologické analýzy: již výše uvedený balíček ade4 je hlavním zástupcem této kategorie

Základní zobrazovací metody v R (base graphic) mají za cíl funkčně, rychle a přehledně zobrazit objekty, což splňují. Mohou však budit rozpačitý estetický dojem, díky čemuž nemusí jít vždy o vhodný prezentační prostředek. Pomocí dodatečných úprav se z výchozích jednoduchých vizualizací dá vytvořit trochu pohlednější výsledek, nicméně možnosti úprav jsou velmi omezené. „Základní grafika funguje podle modelu pero na papíře. Můžete kreslit vždy navrch kresby, ale nelze modifikovat nebo vymazat existující obsah. Současně zde není pro uživatele dosažitelná grafická reprezentace oddělená od samotného zobrazení na obrazovce.“ (Wickham, 2009). To znamená, že v jediném příkazu dochází k nastavení i samotné prezentaci, nelze to provádět zvlášť.

V roce 1999 vydal Leland Wilkinson publikaci The Grammar of Graphics, která se stala inspirací pro vznik grafického balíku ggplot, následně ggplot2. Rozdíl mezi těmito balíky a výchozími zobrazovacími metodami R je v přístupu ke grafice. Konvenční přístup funguje na principu „zobrazit – co”, kdežto ggplot2 bere výsledný grafický výstup jako abstraktní objekt, který rozděluje na jednotlivé, konkrétně modifikovatelné a navzájem nezávislé komponenty. Zjednodušeně řečeno pracuje ve stylu zobrazit(co-jak + co-jak +…).

Wickham (2009) ve své knize ggplot2: Elegant Graphics for Data Analysis rozlišuje několik komponentů, jejichž seznam s ukázkami je možné nalézt v tomto dokumentu.

1. data (datový soubor a základní systém estetických pravidel, např. přiřazení proměnných k osám nebo rozlišení kategorií)
2. geoms (geometrické objekty definující to, co je ve výsledku vidět; způsob zobrazení dat, např. body, linie, polygony, histogram, atd.)
3. stats (statistické transformace, například count pro histogram, nebo bining – plošné součty)
4. scales (měřítka os a definování legendy)
5. coord (systém souřadnic + zobrazení pomocných liníí vedlejších os, výchozí systém bývá kartézský, k dispozici jsou možnosti například polárních nebo mapových souřadnicových systémů)
6. facet (aspekty, nabízí možnosti rozložení datového souboru do podmnožin, například zobrazení grafů pro různé roky zvlášť)
7. další vedlejší komponenty jako theme (témata, obstarává úpravy barev a velikostí všech popisků nebo pozadí), annotation, fortify nebo position.

Import vektorových dat

Spolupráci s formátem shapefile obstarává například balík maptools, který v sobě obsahuje kopii knihovny OGR pro čtení shapefile nebo balík rgdal, komunikující s množstvím dalších vektorových i rastrových formátů. Pokud však další formáty nepotřebujeme řešit, funkce readShapeSpatial() z balíku maptools je pro uživatele jednodušší a pro velkou část operací dostačující. Ke správnému načtení shapefile do prostředí R, konkrétně jako prostorová třída SpatialPolygonsDataFrame balíku sp, pak stačí jako argument zadat název souboru (příp. s cestou). Této funkci není třeba ani deklarovat, zda se jedná o bodovou, liniovou nebo polygonovou vrstvu. Pro export do shapefile pak slouží funkce writeSpatialShape(), kde jsou argumentem objekt z R a požadovaný název souboru na disku. maptools však neřeší čtení nebo zápis souboru s příponou .prj, který u shapefile uchovává informaci o systému souřadnic a projekci, což je pro některé operace nezbytné. Je však možnost tuto informaci k načtenému objektu přiřadit nebo při exportu tento textový soubor manuálně zkopírovat z jiných shapefile stejných souřadnicových systémů. V tomto má výhodu řešení importu a exportu přes balík rgdal, který soubor .prj bere v úvahu. Import provádí funkce readOGR(), kde argumenty jsou složka a soubor na disku bez přípony, a export funkce writeOGR() s argumenty objekt, složka, název exportu a ovladač požadovaného formátu.

Načteme a zobrazíme si vektorovou vrstvu obcí krajů České republiky.

data<-readShapeSpatial(“C:/rFolder/nuts3_010113.shp”) # knihovna maptools
data1 <- readOGR(dsn=”C:/rFolder”, layer=”nuts3_010113″) # knihovna rgdal
plot(data)

Geokódování

Pokud jsou data uložena v tabulce, kde každý záznam zahrnuje atribut adresa, existuje možnost, jak na jejím základě přiřadit k záznamům zeměpisné souřadnice. Zjednodušeně je možné jej popsat následovně: řetězec obsahující adresu v určité syntaxi je poslán na server poskytovatele geokódování, který vlastní databázi adres, a jako odpověď se vrátí souřadnice, opět v určité syntaxi. Existuje i opačný postup, na základě souřadnic vypisování adresy, jedná se o tzv. reverzní geokódování. V rámci R se lze setkat se 3 možnostmi geokódovacích procesů – Google Maps, OSM Nominatinm a mapy.cz.

Google Maps Geocoding API

Nejjednodušší způsob geokódování v R nabízí balík ggmap skrze funkci geocode(), kam stačí zadat jako argument adresu. Volitelným argumentem se nastavuje požadovaný výstup, zda chceme pouze souřadnice x a y, či i s adresou nebo i víc. Poskytovatelem tohoto způsobu geokódování je společnost Google přes svoji aplikaci Google Maps Geocoding API. Provede se zadaný dotaz a do další proměnné se podchytí celý výsledek. Z této proměnné se pak vypreparují požadované výstupy, v tomto případě zeměpisná šířka a délka, (ale lze například vytáhnout i informace o celé adrese nebo o administrativních regionech), které se stanou výstupem celé funkce.

geocode(location, output = c(“latlon”, “latlona”, “more”, “all”), messaging = FALSE, …)

kde location obsahuje údaje o adrese a outpup specifikuje objem výsledků.

souradnice <- geocode(“Ostrava”, output = “more”, messaging = TRUE)

geocodeQueryCheck() # pro kontrolu, kolik dotazů zbývá dotazů na adresu
write.csv(souradnice, “C:/rFolder/souradnice_export.csv”)

lon	lat	type	loctype	address	postal_ code	country	administrative_area_level_2	administrative_area_level_1	locality
18,263	49,821	locality	approximate	ostrava, czech republic	NA	czech republic	ostrava-city district	moravian-silesian region	ostrava

Nevýhodou geokódování přes Google Maps Geocoding API je limit, který povoluje maximálně 2 500 dotazů za den. Další vlastností tohoto způsobu geokódování, která může být zároveň výhodou i nevýhodou, je aproximace zadané adresy. Nicméně v jednom z atributů výsledku je i informace o typu určení souřadnic, při správném zadání ulice nabývá hodnoty GEOMETRIC_CENTER (souřadnice jsou geometrickým středem ulice), při správném zadání kompletní adresy na úrovni domu nabývá hodnoty ROOFTOP a právě při chybně zadané ulici nebo pouze zadání obce či větší administrativní jednotky vypíše APPROXIMATE.

S využitím Google Maps’ Geocoding API je možné provést geokódování následovně:

getGeoCode <- function(gcStr) {

library(“RJSONIO”) #Load Library

gcStr <- gsub(‘ ‘,’%20’,gcStr) #Encode URL Parameters, mezera se nahrazuje za %20

#Open Connection

connectStr <- paste(‘http://maps.google.com/maps/api/geocode/json?sensor=false&address=’,gcStr, sep=””)

con <- url(connectStr)

data.json <- fromJSON(paste(readLines(con), collapse=””))

close(con)

#Flatten the received JSON

data.json <- unlist(data.json)

if(data.json[“status”]==”OK”) {

lat <- data.json[“results.geometry.location.lat”]

lng <- data.json[“results.geometry.location.lng”]

gcodes <- c(lat, lng)

names(gcodes) <- c(“Lat”, “Lng”)

return (gcodes)

}

}

geoCodes <- getGeoCode(“Ostrava, 28. rijna”)

geoCodes

Problémem je také diakritika. Při zápisu adres je potřeba háčky a čárky odstranit. Problémem může být i přesnost lokalizace.

Prostorové operace

V R je k dispozici funkčně obsáhlý balík rgeos, představující rozhraní ke knihovně GEOS (Geometry Engine – Open Source, projekt OSGeo). V rámci něj jsou definovány různé prostorové metody, geometrie a topologické vztahy. Jednou z často využívaných operací je průnik. rgeos nabízí funkce gIntersetion(vrstva1,vrstva2), výsledkem které však je nová vrstva třídy spatial, která sice obsahuje prvky z jedné původní vrstvy, které polohově pronikají s (požadovanými nebo se všemi) prvky z druhé vrstvy, nicméně při tomto procesu dojde ke ztrátě všech atributů.

Balík sp nabízí alternativy pod funkcemi over() a overlay(). Výsledkem over(vrstva1, vrstva2) je matice o stejném počtu entit jako první vrstva a obsahující atributy prvků z druhé vrstvy prostorově se protínajících. To znamená například při over(body, polygon) bude výsledkem matice se stejným počtem entit jako původní bodová vrstva obsahující atributy polygonu, ve kterém se daný bod nachází. Tato tabulka se dá k původní bodové vrstvě připojit, identifikátorem je číslo řádku. Je možné libovolně kombinovat typy geometrie, pokud je větší množství entit z druhé vrstvy v průniku s první vrstvou, ve výsledku budou pouze atributy první z nich. Funkce overlay(vrstva1,vrstva2) je potom jen na průnik bodové a polygonové vrstvy s podobným výsledkem.

vrstva_1 <-readShapeSpatial(“C:/rFolder/nuts3_010113.shp”)

vrstva_2 <-readShapeSpatial(“C:/rFolder/nuts3_010113.shp”)

library(rgeos)

prunik <- gIntersection(vrstva_1,vrstva_2)

par(mfrow=c(2,2)) # funkce par, pro zobrazení několika obrázku o stejne velikosti vedle sebe

plot(vrstva_1)

plot(vrstva_2)

plot(prunik)

Balík rgeos nabízí širokou škálu prostorových operací, merge, buffer, contains. Obsahuje prakticky všechny základní funkce. Pro správné použití dané funkce je vždy vhodné se podívat do nápovědy (např. příklazem ?rgeos).

Užitečnou funkci nabízí i balík sp, jedná se o spsample(), která vytváří nové body v rámci zadaného polygonu, a to buď náhodně, nebo i pravidelně. Dalšími argumenty může být buďto počet nově vytvořených bodů, nebo při pravidelném způsobu i jejich vzdálenosti od sebe. Této funkce lze využít například při vytváření tzv. fishnet, pravidelné polygonové vrstvy ve formě sítě.

Vizualizace dat

Ukážeme si možnosti vizualizace na konkrétním příkladě:

Příklad využití funkce plot() na shp datech. Kombinace liniove, bodové a polygonove vrstvy.

Může být vytvořen také kartogram:

library(maptools)

library(RColorBrewer)

library(classInt)

colours <- brewer.pal(5, “Reds”) # vytvoření proměnné s barvami

kraje <- readShapePoly(“nuts3_010113.shp”)

brks<-classIntervals(kraje$PRUMERNA_M, n=5, style=”jenks”, dataPrecision=1) # vytvoření 5 intervalů z proměnné Jenksovou metodou se zaokrouhlením. style: “fixed”, “sd”, “equal”, “pretty”, “quantile”, “kmeans”, “hclust”, “bclust”, “fisher” nebo “jenks”

brks_n<- brks$brks # uložení hranic intervalů do proměnné

plot(kraje, col=colours[findInterval(hranice$PRUMERNA_M, brks_n, all.inside=TRUE)], axes=T)

title(paste (“Průměrná mzda v krajích ČRn v roce 2013”))

legend(-950000,-1155000,legend=leglabs(brks_n), fill=colours, bty=”n”, cex = 0.7,text.col=”black”)

legend(-950000,-1215000, c(“Hranice kraje”), lty=c(1), lwd=c(2.5), col=c(“black”), bty=”n”, cex = 0.7)

SpatialPolygonsRescale(layout.north.arrow(2), offset= c(-420000,-1200000), scale = 40000, plot.grid=F)

text(-465000,-1220000, “Souřadnicový systém S-JTSKnInstitut geoinformatiky; VŠB-TUO”, cex= 0.5)

text(-905000,-1155000, “Průměrná mzda:”, cex= 0.8)

Vizualizace pomocí balíku ggplot2

Balík ggplot2 přináší nový systém vizualizace. Hlavními výhodami jsou propracované řešení skládání komponent, kontrola detailů, možnost uložení části i celku vizualizačního postupu do proměnné, pohledné výstupy, automatické postupy a spousta dostupných návodů na diskuzních fórech apod. Naopak jeho nevýhodou je nepřímá spolupráce s objekty prostorových tříd, které je nutno převést do některé ze tříd běžného datasetu.

Bodové vrstvy (mající souřadnice v atributech) stačí konvertovat pomocí funkce as.data.frame(). Pro převod liniových a polygonálních vrstev slouží funkce fortify(), kde se jako parametr zadává kromě názvu prostorového objektu také název sloupce obsahující ID instancí. Výsledkem je rozsáhlá matice zlomových bodů se souřadnicemi x, y a se společným identifikátorem pro region nebo linii. Dojde však ke ztrátě případných dalších atributů objektu, které je tak nutno připojit zpět.

Řešením je využití metody s názvem poly_coords() má výhodu zachování všech původních atributů. Potřeba je však, aby prostorový objekt obsahoval sloupec s číselnými identifikátory přesně pojmenovaný ID.

Vhodné je seřadit nově vzniklé body podle atributu order (funkcí order()), občas je toto pořadí proházené a pro správné vykreslování hranic je pochopitelně potřeba běžného vzestupného pořadí.

Nejjednodušší poloautomatické funkci qplot() stačí k výsledku zadat proměnné pro osy x a y a typ geometrie. Funkce qplot() zpracovává některé volitelné atributy jako například titul, názvy os, limity os, rozdělení na kategorie či vykreslení spojité proměnné a další, nicméně pro práci jdoucí více do hloubky grafiky ggplot2 je vhodnější komplexnější funkce ggplot().

Důležitým stavebním kamenem ggplot() je funkce aes(), což je zkratka anglického slova aesthetics značící estetiku. Do ní se vkládají zobrazující informace jako název datového souboru, proměnné pro osy x a y a proměnné tvořící legendu, kategoriální rozdělení zobrazovaného objektu do různých barev (argument colour), velikostí (argument size) a tvaru (argument shape). aes() se dá specifikovat pro celý zobrazovací postup (v kódu se vloží samostatně jako prvek oddělený od ostatních znaménkem +) nebo pro konkrétní geom. Lze použít i obojí, v tomto případě jsou hodnoty aes pro geom nadřazené v rámci onoho geomu, pro zbytek prvků v ggplot() platí oddělené aes.

Práci s bodovými daty

library(ggplot2)

okresy_body <- readShapePoints(“okresy_body.shp”)

okresy_df<-as.data.frame(okresy_body) # konverze z prostorových tříd do běžných dataframe. U bodové vrstvy, ve které jsou obsazeny sloupce souřadnic stačí funkce as.data.frame(). Při jejich absenci (např. import z shp) je nutný obdobný postup jako s polygony pomocí fortify() (vyžaduje mít nainstalovaný balík rgeos).

names(okresy_df)<- c(“LAU1”, “NAZKRAJ”, “NAZOK”, “OBAKT”, “X”, “Y”, “x1”, “Y1”) # pomocí této funkce přejmenujeme atributy tak, aby je následně bylo možné pomocí funkce qplot() využít.

qplot(data = okresy_df, X, Y, geom=”point”, colour=NAZKRAJ, size=OBAKT, xlab=”X”, ylab=”Y”, main=”Počet obyvatel v okresech v roce 2013″) # argumenty barvy a velikosti bodu, názvy os a obrázku. Možno doplnit i o argument shape s 6 různými tvary.

Práci s areálovými daty

Při práci s polygonovými daty je možné data vizualizovat 4 způsoby podle hodnoty argumentu geom.

Zadání pro samostatnou práci

• načtěte si vektorová data okresů, napojte k němu data unemploy.xls a zobrazte sloupec míra nezaměstnanosti
• dokončete mapový výstup

Data:

Data pro vedené a samostatné cvičení

Cvičení je vytvořeno v rámci projektu Inovace bakalářských a magisterských studijních oborů na Hornicko-geologické fakultě VŠB-TUO pod číslem CZ.1.07/2.2.00/28.0308. Tento projekt je realizován za spoluúčasti EU.